CC BY
75
Вычислительные технологии
Том 17, № 6, 2012
Численное моделирование обратной связи p53 — Mdm2 в биологическом процессе апоптоза*
О.Ф. ВОРОПАЕВА, Ю.И. ШОКИН Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: vorop@ict.nsc.ru, dir@ict.nsc.ru
Выполнен численный анализ решений 2D нелинейной математической модели взаимосвязи белков p53 и Mdm2, играющих одну из ключевых ролей в биологическом процессе апоптоза.
Ключевые слова: р53-зависимый апоптоз, p53 — Mdm2-сеть, ОДУ с запаздывающим аргументом, численное моделирование, бифуркация Андронова — Хопфа.
Введение
Молекулярно-генетические механизмы старения связываются с накоплением дефектов в клеточной структуре, в процессах транскрипции и трансляции белков, с возникновением ферментов с дефектным функционированием, с которыми механизмы возмещения не могут справиться. Ответом здорового организма на накопленные дефекты является апоптоз — генетически запрограммированная гибель клеток, не только позволяющая удалить ненужные или дефектные клетки, но и выполняющая центральную роль в поддержании клеточного гомеостаза. Апоптотическая гибель клеток происходит без воспаления, что отличает данный процесс "разумной" самоликвидации клеток от губительного для организма некроза.
Масштабы апоптоза впечатляют: ежедневно в здоровом организме рождается несколько десятков миллиардов новых клеток (за год их общая масса оказывается сопоставима с массой тела) и такое же их количество гибнет в основном за счёт апоптоза. В этих условиях последствия сбоев в апоптозе, с которыми связываются такие тяжелейшие заболевания как нейродегенерация (болезни Альцгеймера, Паркинсона, рассеянный склероз, эпилепсия и др.), остеопороз, остановка роста, преждевременное старение внутренних органов и рак, зачастую оказываются катастрофическими. Изучению апоптоза посвящены работы [1-20], содержащие более полный обзор литературы по данному вопросу. Особое место в современных исследованиях занимает связь апоптоза, старения и рака [1].
Одну из ключевых ролей в апоптозе играет белок p53, кодируемый геном p53. Получая сигналы о клеточном повреждении, белок p53 либо останавливает клеточный цикл для устранения повреждения (восстановления генома), либо инициирует апопто-тическую гибель клетки. Одновременно белок р53 ("фактор некроза опухоли") широко известен как важный подавитель опухоли — согласно наблюдениям, его функция нарушена более чем в 50 % опухолей человека. Лабораторному изучению р53 посвящено значительное число работ, ссылки на которые можно найти, в частности, в [1-9].
* Работа выполнена при финансовой поддержке гранта Президента РФ НШ-6293.2012.9.
Считается, что в нормальных клетках содержание и активность р53 должны быть низкими, поскольку чрезмерная активация р53 способствует преждевременному старению или развитию хронических патологий. С другой стороны, снижение функции р53 должно быть своевременно компенсировано, чтобы обеспечить необходимый ответ р53 на возможный стресс. Таким образом, в организме должен действовать механизм регулирования как избыточного продуцирования р53, так и его чрезмерного подавления (последнее повышает риск развития рака).
Одно из центральных мест в этом механизме регулирования занимает взаимосвязь р53 с белком Мёш2, который способствует угнетению р53 [1-9]. В ответ на сигналы о повреждении ДНК белок Мёш2 быстро деградирует, что позволяет р53 накапливаться и полностью активироваться. Одновременно ген Мёш2 активируется при повышении уровня белка р53, способствуя его разрушению, что важно для восстановления низкого уровня р53 после завершения стресса. Поэтому Мёш2 имеет решающее значение для содержания р53 под контролем. Необходимо отметить, что описанный механизм представляет собой лишь весьма упрощённую схему апоптоза, за каждым элементом которой кроется сложнейшая система взаимодействий. Многие аспекты данного процесса остаются пока малоизученными или неизвестными и находятся в зоне пристального внимания исследователей.
Математическое моделирование взаимосвязи р53 — Мёш2, основанное на моделях разной степени сложности, осуществлялось в целом ряде работ [9-21]. Основной принцип построения этих моделей — получение качественного согласия с лабораторными исследованиями, свидетельствующими о колебательном характере процесса при воздействии на клетки УФ-излучением (см., например, [9, 12, 14]) и запаздывании реакции Мёш2 на изменение активности р53. В работах [16, 18, 19] математическое моделирование дополняет углублённый теоретический анализ используемых моделей.
Интересным примером нелинейной 2В-модели взаимосвязи р53 — Мёш2, включающей в себя ОДУ с запаздывающим аргументом, является модель, предложенная в [10] (см. также [13]). В работах [10, 14, 16] представлены некоторые теоретические оценки решений для данной модели. В [10] выполнены численные эксперименты, в ходе которых изучалось возможное поведение р53 и Мёш2 под воздействием стресса, моделируемого как скачок в значении одного из параметров математической модели при больших значениях времени, когда решение системы в дострессовом состоянии выходит на стационарный режим. Полученные результаты [10, 13] указывают на возможность таких ситуаций, когда существенное уменьшение одного из параметров системы — константы диссоциации комплекса р53 — Мёш2 — запускает механизм реагирования, характеризующийся возникновением периодических колебаний в системе р53 — Мёш2, которые могут позволить клеткам восстановить свою ДНК, не рискуя необратимыми последствиями продолжающейся чрезмерной активации р53. В работах [10, 15] обсуждается реакция системы на изменение сразу двух параметров модели [10] и приводятся результаты численных экспериментов, в ходе которых выделялись только режимы, иллюстрирующие колебательный характер взаимосвязи белков р53 — Мёш2. Неполнота и противоречивость ряда полученных данных [15] послужили поводом для настоящей работы, целью которой является детальное численное моделирование взаимосвязи р53 — Мёш2 путем исследования возможных решений системы уравнений, представленной в [10], и способов управления этой системой. Основное внимание уделено методическим расчётам при достаточно малом времени запаздывания в окрестности базальных значений коэффициентов системы.
1. Постановка задачи
Математическая модель взаимосвязи белков р53 и Mdm2 представляет собой следующую систему уравнений с запаздывающим аргументом [10]:
= s - af (yi(t),y2(t)) - byi(t), (1)
dd2 = cig(yi(t - T),y2(t - т)) - C2y2(t), (2)
где взаимодействие белков определяется функциями
f(УЪ y2) = 2 (yi + У2 + k/ - \J(yi + y2 + k/)2 - , (3)
/ ч yi - f(УЬЫ ...
£(УьУ2Н^-^-7. (4)
yi + kg - f (yi,y2)
Здесь yi и y2 — концентрации р53 и Mdm2 соответственно; s — скорость производства р53; a — скорость деградации р53 посредством убиквитинирования (присоединения к белку-мишени небольшого белка убиквитина для инициирования деградации, а также для регуляции локализации и функционирования этого белка); b — скорость самопроизвольного распада р53; ci — скорость производства белка Mdm2, в том числе за счёт взаимодействия с р53; c2 — скорость деградации белка Mdm2; k/ — константа диссоциации комплекса p53 — Mdm2; kg — константа диссоциации белка p53 и гена Mdm2 (более детальное биологическое описание параметров см. в [10, 13]).
Начальные данные для системы (1)-(4) задаются в виде функций "истории"
yk(0) = Pk(0), 0 е [-т, 0] , k = 1, 2. (5)
Значения параметров
s = 1 с-\ a = 3 ■ 10-2 с-1, b = 10-4 с-1,
ci = 1 с-1, c2 = 10-2 с-1, kf = 180, kg = 28 (6)
согласованы с данными различных лабораторных исследований и вслед за [10, 13, 15] принимаются как базальные, т. е. соответствующие нормальному состоянию здорового организма (обсуждение правомерности отнесения значений (6) к категории базальных выходит за рамки настоящего исследования). Параметр т определяет время запаздывания реакции Mdm2 на изменение состояния белка р53 и, согласно наблюдениям, может варьироваться от нескольких минут до нескольких часов в зависимости от фактора стресса системы (как правило, в лабораторных исследованиях факторами стресса, инициирующими апоптоз, являются ультрафиолетовое и другие виды облучения клеток).
2. Методы решения и тесты
Система (1)-(4) обладает рядом особенностей, которые необходимо учитывать при выборе численного алгоритма. К их числу относится наличие запаздывания, вследствие чего одной из ключевых является проблема задания функций "истории" (5). При этом
запаздывание т является варьируемым параметром, от значения которого решение зависит самым существенным образом. Отметим также проблему нарушения гладкости решения в точке £ = 0 и последующего переноса разрывов производных в точки, кратные величине запаздывания т. Кроме того, согласно оценкам [10, 15], при некоторых условиях, связанных со значительными отклонениями значений параметров от базаль-ных (6), система (1)-(4) может быть жёсткой. Однако подобные ситуации останутся за рамками настоящей работы, поскольку авторов в первую очередь интересуют решения, соответствующие достаточно близким к базальным значениям параметров, для которых жёсткость не характерна [10]. Тем не менее при проведении расчётов предпочтение отдавалось неявным численным методам и особое внимание уделялось контролю за шагом сетки и его влиянием на решение.
В основе процедуры решения задачи (1)-(5) лежит широко известный (в первую очередь как метод проведения теоретических исследований) метод шагов [22, 23]. Его суть сводится к следующему: вначале решение уравнений выполняется на интервале (0, т] с привлечением функций "истории" ^\(9), ^2(е), заданных на интервале [—т, 0]; далее расчёты проводятся последовательно на интервалах размером т с "историей", полученной как решение системы на предыдущем интервале. Для простоты численной реализации во всех расчётах величина запаздывания т полагалась кратной шагу расчётной сетки (при условии, что величина шага является достаточно малой величиной, это ограничение не оказывает существенного влияния на характер решения).
Для численного решения использовались следующие известные методы: метод предиктор-корректор, основанный на методах Эйлера и Рунге — Кутты, методы Адамса и Рунге — Кутты третьего и четвёртого порядка, одноитерационный метод типа метода Розенброка третьего порядка, методы Гира разных порядков. Дополнительно были выполнены расчёты модифицированным методом Адамса, специально предназначенным для решения ОДУ с запаздыванием (с корректировкой решения в окрестности разрывов производных) [24], а также методом Адамса — Бэшфорта — Моултона, в рамках которого решение, вычисленное по явному методу Адамса четвёртого порядка, используется затем в неявной реализации этого метода. На первом этапе анализ работоспособности методов проводился в численных экспериментах с простейшим уравнением с запаздыванием
Ит
= -т(£ - т), т(е) = 1, е е [-т, 0] , т =1,
имеющим точное решение
[ Т ]+1 т (£)= £ (-1)
(£ - (П - 1) т)п
п!
п=0
Получено, что наиболее предпочтительным в смысле точности и простоты численной реализации является метод Адамса. В частности, в расчётах на сетке с 960 узлами (при £ е [0, 12]) относительная погрешность решения в равномерной норме по методу Адамса составила 2.66683 • 10-7, по методу Рунге — Кутты — 2.42773 • 10-4, по методу типа метода Розенброка — 5.66895 • 10-3. Вычислительные эксперименты показывают, что нарушение гладкости решения в точке £ = 0 может компенсироваться при расчётах на достаточно больших интервалах времени (см. также [22, 23]). Это указывает на возможность успешного применения обычных численных методов решения систем ОДУ в условиях, когда запаздывание является одним из варьируемых параметров системы.
п
Основываясь на результатах тестовых расчётов, для решения задачи (1)-(5) использовался метод Адамса. Дополнительно (в методических целях) применялись методы предиктор-корректор, Рунге — Кутты и Адамса — Бэшфорта — Моултона. Отметим, что недостаток многошаговых методов, связанный с необходимостью особого подхода к вычислению решения в первых точках, при решении уравнений с запаздыванием перестаёт быть таковым в силу особенностей задания начальных условий. Во всех случаях для уточнения решения системы (1)-(5) выполнялись итерации по нелинейности. При этом результаты, полученные по указанным методам, оказались весьма близкими. Расчёты системы (1)-(4) с начальными данными (5) проводились на сетках с разрешением от 10 до 10 000 точек на 1 с, а основная рабочая сетка имела 100 точек на 1 с.
3. Численное решение при базальных значениях параметров
Прежде чем приступить к исследованию работоспособности принятой математической модели (1)-(5) под влиянием факторов, которые в литературе определяются как стрессовые по отношению к нормальному состоянию здорового организма, необходимо получить более детальное представление о решении, соответствующем нормальному состоянию (6).
Была проведена серия расчётов задачи (1)-(6) с нулевыми значениями функций "истории" ^к на интервале [—т, 0] (такой способ реализован, в частности, в работе [10]). В базальных условиях единственным ключевым параметром, способным повлиять на решение, является параметр запаздывания. Наши расчёты показали, что при т = 0 поведение решения можно охарактеризовать как весьма стремительный переход с нулевых значений концентраций белков на новый уровень — стационарные предельные значения у* и 154.6347 (р53), у* и 81.3105 (Мёш2).
При варьировании параметра запаздывания т возникают затухающие колебания концентраций, иллюстрирующие взаимодействие белков по принципу обратной связи. При этом было установлено, что решение системы (1)-(6) имеет состояние равновесия у*, у*, к которому сходятся все затухающие решения, соответствующие достаточно широкому диапазону значений параметра запаздывания. Таким образом, полученное при т = 0 решение является единственным устойчивым пределом на рассмотренном интервале значений т (эти данные соответствуют известным представлениям о поведении затухающих решений динамических систем, см. [18]). На рис. 1 представлены некоторые результаты расчётов (здесь и далее в тексте и на рисунках у1 — р53, у2 — Mdш2). Эти данные иллюстрируют начальный этап линейного роста на котором достигает-
ся пиковое значение концентрации р53, обусловленное первоначально низким уровнем Mdш2 и запаздыванием реакции Mdш2 на рост р53 (величина участка роста прямо пропорциональна величине запаздывания). Спустя время т начинается рост у2(£) до уровня, необходимого для запуска механизма обратной связи, который характеризуется затухающими колебаниями обеих компонент решения (практически в противофазе) и сходимостью этого процесса к предельным точкам у* ~ 154.6347, у* ~ 81.3105 (рис. 1, а). На рис. 1, б приведены максимальные значения концентраций укт = тах (у к(¿)) (к = 1, 2) в зависимости от т. Видно, что с увеличением т пиковые значения концентрации р53 также увеличиваются, а рост максимальных значений концентрации белка Mdш2 ограничен сверху при всех т > 1200 с. При достаточно больших значениях времени можно наблюдать весьма продолжительные по времени участки релаксации с близкими к постоянным значениями концентраций.
У1>У2
400
200
0
О
2000 4000 6000 8000 с
100
1000
Рис. 1. а — решение задачи (1)-(5) с базальными значениями параметров (6) в зависимости от времени и параметра запаздывания: линии 1-6 соответствуют т = 52, 120, 240, 480, 1480, 2480 с, пунктир — т = 0; б — максимальные значения компонент решения как функции параметра запаздывания т (при базальных значениях коэффициентов модели): 1 — ухт, 2 — у2т
Дополнительно была проведена серия численных экспериментов, в которых в качестве начальных данных использовались стационарные значения концентраций у1 ~ 154.6347, Уз ~ 81.3105 с внесением в них "стрессовых" возмущений — от 20% до двух порядков (в сторону убывания и возрастания поочередно для каждой компоненты решения и одновременно для обеих). При этом варьировался также параметр запаздывания. Расчёты показали, что при всех рассмотренных стрессовых ситуациях характер решения системы (1)-(6) меняется слабо: решение возвращается к установленным ранее предельным значениям уЗ, уЗ, т. е. имеет устойчивый фокус. Полученные в расчётах данные о стационарном решении задачи имеют важное значение для понимания моделируемого процесса и более тщательного моделирования функций "истории".
Описанные численные эксперименты проводились для т Е [0, 5000] (с). При дальнейшем увеличении параметра запаздывания не удалось получить других типов решений. Можно отметить, что период затухающих колебаний возрастал пропорционально значению т. При этом математическая модель всегда обеспечивала запуск биологического механизма отрицательной обратной связи в виде затухающих колебаний двух концентраций, происходящих в противофазе и характеризующихся наличием при больших значениях времени предельных значений экстремумов концентраций.
4. Зависимость решения от параметров системы
Проведено численное исследование поведения решений системы (1)-(4) при отклонении значений основных параметров системы в, а, Ь, с\, с2, kf, кд от базальных (6). Ряд результатов подобных исследований представлен в работах [10, 15, 25]. В [10, 15] основное внимание уделялось варьированию значений в, kf, кд. Было показано, что для некоторых комбинаций указанных параметров решение имеет предельные циклы (представленные в [15] данные требуют более тщательного рассмотрения из-за имеющихся
в этой работе несоответствий зависимостей у1(^) и у2(£) и фазовых портретов решений). В настоящем разделе будут представлены результаты детального численного исследования решений, характерных для фиксированного набора параметров системы в условиях отклонения значения одного из них от базального (6). Основные расчёты этой серии выполнялись с нулевыми значениями функций "истории" при параметре запаздывания т = 120 с (выбор данного значения обусловлен соображениями удобства проведения методических расчётов, на которые нацелена настоящая работа).
Для оценки влияния изменения скорости в производства белка р53 на решение задачи проводились расчёты, в которых значение в варьировалось как в сторону возрастания, так и в сторону убывания (дополнительно варьировался параметр запаздывания т). Отклонения параметра в от вЬа8а1 = 1 с-1 составляли от 20 % до четырёх порядков. Некоторые результаты расчётов представлены на рис. 2 (здесь и на рис. 3 точками помечены линии, соответствующие решению при базальных (6) значениях параметров). При незначительном завышении значения в получены либо монотонно возрастающие решения, либо решения в виде затухающих колебаний, обладающие собственным устойчивым фокусом. Видно также, что существенное завышение значения параметра я может приводить к резкому росту концентрации р53, при котором ингибитор (Mdш2) практически не оказывает должного регулирующего воздействия. В свою очередь снижение я может иметь негативные последствия в виде чрезмерно низкого (в сравнении с базальными стационарными значениями) содержания обоих белков, особенно р53. Рост времени запаздывания приводит к изменениям в решении, аналогичным описанным в предыдущем разделе.
При варьировании значений скорости деградации белка р53 под действием Mdш2 (параметра а) были получены следующие результаты. Уменьшение а приводит к монотонному возрастанию концентраций у1(^) и у2(£). При этом возможны ситуации чрезмерного роста концентрации белка р53, в которых регулирующие функции Mdш2 не реализуются. При увеличении данного параметра при всех т Е [0, 5000] (с) удалось получить только затухающие решения, обладающие собственным устойчивым фокусом. При значительном увеличении параметра а существенную роль могут играть релаксационные эффекты. Характерные для данной ситуации решения (в фазовой плоскости) представлены на рис. 2, б. Здесь и ниже расчётный интервал времени составлял от 32 000 до 100 000 с.
В ходе численных экспериментов варьировался также параметр Ь — скорость самопроизвольного распада р53. Представление полученных решений в фазовой плоскости позволяет наглядно продемонстрировать падение концентрации обоих белков при существенном (на два — четыре порядка) увеличении параметра Ь, что в итоге может привести к полному угасанию обратной связи р53 — Mdш2 (аналогичный эффект можно наблюдать и при варьировании других коэффициентов уравнения (1), см. рис. 2). Расчёты показали, что базальное значение параметра Ь представляет собой минимальный предел: при дальнейшем уменьшении Ь решение (при фиксированном т) практически не меняется. Так, на рис. 2, в линии 7-9, иллюстрирующие решение при уменьшении параметра Ь на один — три порядка, оказались весьма близкими друг к другу. Таким образом, при варьировании параметра Ь решения в виде неконтролируемого роста р53 обнаружены не были. При этом с ростом т более выраженными становятся релаксационные эффекты.
При существенном (на порядок и более) уменьшении базального значения скорости производства белка Mdш2 (параметра с1) возможны ситуации чрезмерного роста
уровня p53, когда уровень Mdm2 не поднимается выше некоторого достаточно низкого предельного значения, и, напротив, при таком же увеличении 0\ наблюдаются в основном затухающие решения с весьма низкими стационарными значениями p53, обусловленными резким ростом уровня Mdm2 (рис. 3, а). При фиксированном запаздывании (т = 120 е) можно указать значение 0\ ~ 0.64 с-1, разделяющее эти две ситуации. При больших значениях т и с1 проявляются релаксационные эффекты. На рис. 3, б представлены также решения, соответствующие различным значениям с2 (скорости деградации Mdm2). Расчёты показали, что характер изменения решений аналогичен наблюдаемому при варьировании параметра с1 .
При изменении значения параметра кд (константа диссоциации белка р53 и гена Mdm2) относительно базального при т = 120 е были получены только монотон-
Рис. 2. Решение при отклонении параметров системы от базальных значений (6): а — в = д§Ьа8а1 при ^ = 0.5, 1, 1.5, 2, 2.1, 2.2, 2.4 (линии 1-7, т = 120с); б — а = ^аЬа8а1 при q = 3, 1.5, 1, 0.833, 0.667, 0.588, 0.5 (линии 1-7, т = 120с); в — Ь = qЬbasa1 при q = 100, 50, 35, 10, 5, 1, 0.1, 0.01, 0.001 (линии 1-9, т = 120с); штрих — q = 100, т = 1200с
но возрастающие решения или затухающие колебания. Результаты расчётов представлены на рис. 3, в. Возрастающие решения, не допускающие возникновения обратной связи p53 — Mdm2, получены при значительном завышении значения параметра kg относительно базального. Этим решениям, в частности, могут соответствовать ситуации чрезмерного роста уровня p53. При уменьшении kg более чем на два порядка решение практически не меняется, в результате чего фиксируется уровень p53 (весьма близкий к базальному), что является препятствием для неконтролируемого роста Mdm2.
Варьирование значения параметра kf (константа диссоциации комплекса p53 — Mdm2) относительно базального значения kbasal = 180 показало, что в окрестности
Рис. 3. Решение при отклонении параметров системы от базальных значений (6): а — С1 = дс^ при д = 0.1, 0.5, 0.6, 0.64, 0.8, 1, 1.2, 2, 5 (линии 1-9, т =120 с) и д = 1, 2 (линии 10-11, т = 280с); б — С2 = дсЬа8а1 при д = 0.2, 0.833, 1, 1.2, 1.425, 1.5, 1.7, 2 (линии 1-8, т = 120с); в — кд = дк^1 при д = 0.001, 0.01, 0.2, 1, 5, 10, 15, 20, 30 (линии 1-9, т = 120с); г — к/ = дкЬа5а1 при д = 0, 0.01, 0.04, 0.2, 1, 2, 2.5, 3, 4, 5 (линии 1-10, т = 120 с)
этого значения имеются затухающие решения, обладающие фокусами (рис. 3, г). Как и для большинства других параметров системы (1)-(5), для kf имеется пороговое значение, с превышением которого запуск механизма обратной связи невозможен: концентрация р53 растет независимо от Мёш2. Для т = 120 с данная величина оказалась весьма близкой к значению kf ~ 2.5^^' = 450. При этом наряду с уже упомянутыми видами решений при уменьшении kf возможны также решения, имеющие предельный цикл. Для запаздывания т = 120 с значение к^ = ^^'/40 = 4.5 весьма близко к бифуркационному, при котором затухающие решения с предельной точкой переходят в периодические решения с предельным циклом (бифуркация Андронова — Хопфа). Все периодические решения при kf < Щ образуют вложенную последовательность предельных циклов с выраженными релаксационными явлениями; максимальный предельный цикл соответствует kf = 0 (линия на рис. 3 г, помеченная ромбами). Бифуркационные значения kf существенным образом зависят от величины запаздывания.
На рис. 4 приведены неподвижные точки численных решений, полученных при поочередном варьировании значений коэффициентов системы (1)-(4). Все рассчитанные координаты неподвижных точек в фазовой плоскости с достаточно высокой точностью согласуются с аналитическим стационарным решением задачи [10] — точкой пересечения кривых (при соответствующих значениях коэффициентов уравнений (1)-(2)):
(1) _ С1 (а + Ь)у1 — в _ а ( _akg
у21) = -Ч ^ , „, = -М 1 -
с2 (а + Ь)у1 — в + akg с2 \ (а + Ь)у1 — в + а^ (2)_ (в — Ьу1)((а + Ь)у1 + akf — в) _ (в—Ьу) Л , akí
у2 = „// , гл„ „^ = „ 1 +
а((а + Ь)у1 — в) а \ (а + Ь)у1 — в/
На рис. 4 представлены также кривые аналитического стационарного решения [10] при базальных значениях параметров (6). Отметим, что область значений у1 и у2, где стационарное решение [10] имеет особые точки, в настоящих расчётах не рассматривается в силу нефизичности значительной части этих решений (отрицательные значения
Рис. 4. Неподвижные точки решений при поочерёдном варьировании значений коэффициентов; маркированные линии — численные расчёты, штриховые линии — аналитическое стационарное решение [10] при базальных значениях параметров (6)
концентраций). Из рис. 4 видно, что геометрическое место неподвижных точек, полученных при варьировании коэффициентов уравнения (1), в фазовой плоскости представляет собой возрастающую кривую, которая при у1 > 100 практически совпадает с теоретической кривой у21)(у1). Неподвижные точки решений, полученных при варьировании коэффициентов уравнения (2), образуют в фазовой плоскости убывающую кривую, весьма близкую к теоретической у22)(у1). Точка пересечения двух семейств кривых соответствует неподвижной точке — стационарному решению при базальных значениях параметров (6): у* ~ 154.6347, у* ~ 81.3105. Примечательно, что аналитическая кривая оказалась близкой к весьма грубой (на этом участке) приближённой кривой (линия, маркированная полыми квадратами), образованной стационарными решениями при малых к/, которые обладают не фокусом, а предельными циклами.
Таким образом, в рамках принятой модели при отклонении значения только одного из параметров от базального могут возникать следующие ситуации: 1) значительное накопление р53 при весьма низком уровне его ингибитора Mdш2; 2) значительное накопление Mdш2, приводящее к заниженному уровню р53; 3) запуск периодического процесса; 4) быстрое затухание колебаний концентраций белков с выходом на постоянные стационарные значения. В первых двух случаях речь идёт о монотонно возрастающих решениях. Отметим, что чрезмерное накопление р53, как и нежелательное затухание его активности, может приводить к весьма серьёзным негативным последствиям для организма в виде преждевременного старения органов — в первом случае или возрастания риска онкологических заболеваний — во втором. Расчёты показали, что запаздывание является одним из ключевых параметров, влияющих на состояние системы р53 — Mdш2. Для большей части возникающих ситуаций характерны затухающие решения, что указывает на достаточно устойчивую реакцию системы на отклонение параметров от нормальных значений. Отметим, что полученные данные не противоречат известным представлениям о моделируемом явлении и позволяют сделать вывод об адекватности рассмотренной математической модели взаимосвязи р53 — Mdш2. В частности, наличие лишь нескольких пиков в концентрации р53 согласуется с известными данными лабораторных экспериментов (см., в частности, [6, 7]).
5. Периодические решения системы
При проведении следующей серии численных экспериментов основное внимание уделялось поиску условий, приводящих к периодическим решениям, поскольку именно такие состояния системы характеризуют нормальную реакцию организма на стрессы. Появление периодических решений в данной модели связано с бифуркациями решений из неподвижной предельной точки в предельный цикл. В работе [10] изучалось возможное поведение р53 и Mdш2 под воздействием стресса, моделируемого как скачок в значении одного из параметров математической модели при больших значениях времени, когда решение системы в дострессовом базальном состоянии выходит на стационарный режим. Было показано, что механизм обратной связи белков р53 и Mdш2 в виде периодических решений реализуется в рамках модели (1)-(5) только при снижении значения константы диссоциации к/ относительно базального, при этом запаздывание также играет роль бифуркационного параметра. Эти данные подтвердились результатами настоящих расчётов, причём стало очевидно, что реализованный в [10] вариант стресса фактически можно рассматривать как иллюстрацию устойчивости решения си-
стемы к возмущению начальных данных: использованные в настоящей работе нулевые значения функций "истории" и ненулевые стационарные значения ^ (0) = у^ в работе [10] дают весьма близкие по амплитуде и периоду решения. Вместе с тем до сих пор за рамками обсуждения остается вопрос о том, действительно ли малость константы диссоциации kf (при подходящем значении параметра запаздывания) гарантирует возникновение периодических решений.
С целью исследования этого вопроса были выполнены расчёты, в которых, как и в [10, 15], варьировались одновременно два параметра — 5 и kf — вокруг базальных значений. В ходе расчётов для каждого из 5 = 0.01, 0.1, 10 (с-1) параметр kf принимал значения от 0.18 до 1800. Значения остальных параметров, если не оговорено иное, оставались базальными (6). Подобные исследования проводились и в работах [10, 15], где были определены также собственные числа и показаны некоторые характерные периодические решения. В настоящей работе указанные ситуации будут исследованы более детально. Основные расчёты выполнены для т Е (0; 5000] (с).
Рассмотрим решения, соответствующие фиксированному 5 = 0.01 с-1. При kf = кЬава1/1000 = 0.18 и малых т получены затухающие колебательные решения с фокусом, а при всех т Е [150; 5000] (с) — только периодические решения с периодом Т > 2т (последнее согласуется с известными оценками характеристик периодических решений в динамических системах [18]). Полученная в наших расчётах последовательность предельных циклов представлена на рис. 5, а (в центре цифрой 1 показан фокус, соответствующий решению при т = 10с). Отметим, что предельные циклы, наблюдаемые при значениях т, близких к бифуркационному, с достаточно высокой степенью достоверности можно охарактеризовать как эллипсы. При больших значениях параметра запаздывания наблюдается трансформация предельных циклов, обусловленная релаксационными эффектами: фазовые траектории состоят из участков быстрого убывания у2 (¿) и чрезвычайно медленного изменения уровня концентрации у1 (¿) в окрестности
Рис. 5. а — решение при й = 0.01 с 1, kf = 0.18, линии 1-5 соответствуют значениям т =10, 150, 200, 250, 500 с; б — решение при й = 0.1 с-1, kf = 0.18, линии 1-6 соответствуют значениям т =10, 59, 65, 100, 150, 250 с
минимума данной функции. При этом во всех случаях решения у1(^) и у2(£) имеют ненулевые предельные нижние границы значений концентраций. При kf = ^^/100 = 1.8 бифуркация Андронова — Хопфа происходит при т ^^ 550 с, а при всех т Е [550; 5000] (с) получены только периодические решения с периодом Т > 2т. При kf = ^^'/10 = 18 бифуркационное значение параметра запаздывания возрастает, в результате чего периодические решения наблюдаются при всех рассмотренных т Е [3000, 10 000] (с).
Аналогичная зависимость решений от изменения параметров kf и т получена и при в = 0.1 с-1. Так, при kf = 1.8, как и в предыдущих примерах, наблюдается бифуркация Андронова — Хопфа, а бифуркационное значение параметра запаздывания составляет т ~ 210 с, при kf = 18 — т ~ 600 с и далее, вплоть до рассмотренного в расчётах значения т ~ 20 000 с, решение остается периодическим с периодом Т ^ 2т. На рис. 5, б для сравнения с рис. 5, а представлены предельные циклы (цифрой 1 показан фокус решения при т = 10с), соответствующие в = 0.1 с-1 и kf = 0.18 при различных значениях запаздывания т. При в = 0.1 с-1 размер предельных циклов оказался существенно больше, чем при в = 0.01 с-1 для одних и тех же значений т, а релаксационные эффекты стали проявляться раньше.
Проведённые численные эксперименты показали, что при достаточно малых значениях в с приближением к kf = kbasa1 = 180 бифуркационные значения параметра т достаточно быстро возрастают. Так, при в = 0.1 с-1 периодические решения наблюдались лишь при существенно больших значениях т Е [7000, 20 000] (с). Это обстоятельство указывает на то, что устойчивые при умеренном (от десятков секунд до одного часа) запаздывании решения в случае значительного возрастания т могут терять устойчивость. Вместе с тем для более детального анализа этих ситуаций необходимы дополнительные численные эксперименты и теоретический анализ системы (1)-(4).
В ходе экспериментов при в = вЬайа1 • 10 = 10 с-1 для всех рассмотренных kf (в том числе kf = 0) периодические решения не обнаружены: большую часть решений можно охарактеризовать как переход концентраций на новый постоянный стационарный уровень (даже при т ~ 20000с решение имеет только один пик), который определяет неподвижную предельную точку решения. В случае kf = kbasa1 • 10 = 1800 при всех рассмотренных в Е [0.01,10] (с-1) и т были получены только затухающие решения, имеющие неподвижную точку. Таким образом, существенно завышенные значения в и kf (при умеренных т, принятых в данных численных экспериментах) практически исключают появление периодических решений, если остальные параметры принимают базальные значения. Малость параметра kf не гарантирует существование периодических решений, а наблюдаемая в [10] реакция системы на уменьшение kf в 15 раз в виде бифуркации Андронова — Хопфа стала возможной лишь при условии сохранения базальных значений остальных коэффициентов системы.
На рис. 6 представлены некоторые дополнительные данные для случая kf = 0. Видно, что даже при больших значениях в возможны периодические решения. Так, в одном из "предельных" случаев — в = 20 с-1 — бифуркационным параметром является с1 (рис. 6, а), а при базальном значении с1 параметр в сам является бифуркационным (рис. 6, б). В последнем случае бифуркация решения от предельной точки к предельному циклу происходит при в ^^ 2.15 с 1 (периодические решения наблюдаются и при дальнейшем уменьшении в). Было установлено, что при в = 20с-1 и kf = 0 бифуркационными являются также коэффициенты а и с2.
При варьировании параметра в (для простоты полагалось kf = 0, т = 180 с) исследовались также характеристики решений, обладающих предельными циклами. Установле-
но, что с увеличением 5 период колебаний возрастает (ранее была найдена зависимость этой характеристики от величины запаздывания т). Амплитуда колебаний обеих компонент решения также зависит от изменения параметра 5. Особенностью полученных решений является отсутствие монотонной зависимости амплитуды колебаний у1(£) от параметра 5 — максимальная амплитуда достигается не при максимальном 5, как на начальном этапе развития процесса, а при среднем (из рассмотренных) значении 5 = 2 с-1 (на возможность таких ситуаций указывалось в [10, 13]). Следует отметить существование практически одинакового для всех рассмотренных 5 верхнего предельного значения
Рис. 6. а — решение при 8 = 20 с 1, kf = 0 и варьируемом значении С1 = 1,5, 6.52, 7, 10, 20, 30 с-1 (линии 1-7, т = 180с); б — решение при kf = 0 и варьируемом значении 8 =1, 2.15, 2.3, 2.8, 3 с-1 (линии 1-5, т = 180с)
Рис. 7. Решение при 8 = 2 с 1 и варьируемом значении kf = 0, 10, 100, 200 (линии 1-4, т = 180с); сплошные линии — Ыёш2, штриховые — р53
концентрации y2(t). В случае стационарных решений с неподвижной предельной точкой этот факт согласуется с результатами теоретического анализа системы (см. рис. 4). В случае периодических решений подобное поведение также наблюдается. В частности, в расчётах при s = 2 c-1, т = 180 c концентрация y2(t) имела практически один и тот же максимум y2m(t) для всех малых kf, при которых решение имеет предельные циклы. Примечательно, что это максимальное значение весьма близко к стационарным значениям ^2, соответствующим решениям с неподвижными точками для всех kf > kjasal (рис. 7). Таким образом, для параметров системы имеются пороговые значения, дальнейшее изменение которых уже не приводит к повышению уровня Mdm2.
Выполненные численные эксперименты имеют прежде всего методический характер. Более детальное изучение решений рассмотренной системы, в том числе сопоставление с данными лабораторных измерений, представляет собой предмет будущих исследований. Вместе с тем, несмотря на относительную простоту используемой математической модели, выявленные в ходе численных экспериментов взаимосвязи наглядно демонстрируют, насколько сложны и многогранны протекающие в живом организме медико-биологические процессы.
Заключение
В работе выполнен численный анализ решений системы уравнений, описывающей динамику концентраций белков p53 и Mdm2 в живом организме. Проведено детальное сопоставление с аналитическим стационарным решением задачи и расчётными данными [10, 13]. Результаты расчётов качественно согласуются с известными лабораторными измерениями. Проведено детальное численное исследование решения при базальных значениях параметров и отклонениях параметров от нормальных значений. В рамках этих численных экспериментов получены состояния системы, соответствующие как угрозе неуправляемой апоптотической гибели клеток (при чрезмерном росте уровня "фактора некроза опухоли"), ускоряющей процессы старения организма, так и ситуации чрезмерного подавления апоптоза, при которой увеличивается риск возникновения рака. Получены новые сведения о возможности описания в рамках используемой математической модели механизма обратной связи белков p53 и Mdm2, гарантирующего адекватную реакцию организма на серьёзные повреждения ДНК в виде регулируемой гибели клеток.
Авторы выражают благодарность д.м.н. С.Р. Сенчуковой за полезные обсуждения, а также студентам ММФ НГУ К. Бродту, А. Кудряшову, А. Мельниченко и И. Санкову, принимавшим участие в численных экспериментах на разных этапах их проведения.
Список литературы
[1] Lane D., Levine A. p53 research: The past thirty years and the next thirty years // Cold Spring Harb. Perspect. Biol. 2010; 2:a000893.
[2] Kastan M.B., Onyekwere O., Sidransky D. et al. Participation of p53 protein in the cellular response to DNA damage // Cancer Res. 1991. Vol. 51. P. 6304-6311.
[3] Lu X., Lane D.P. Different induction of transcriptionally active p53 following UV or ionizing radiation: Defects in chromosome instability syndromes? // Cell. 1993. Vol. 75. P. 765-778.
[4] Haupt Y., Maya R., Kazaz A., Oren M. Mdm2 promotes the rapid degradation of p53 // Nature. 1997. Vol. 387. P. 296-299.
[5] Moll U.M., Petrenko O. The Mdm2 —p53 Interaction // Mol. Cancer Res. 2003. No. 1, P. 1001-1008.
[6] Lahav G., Rosenfeld N., Sigal A. et al. Dynamics of the p53 —Mdm2 feedback loop in individual cells // Nature Genetics. 2004. Vol. 36, No. 2. P. 147-150.
[7] Абраменко И.В., Злвгородняя А.В., Балан В.И. и др. Индукция p53-зависимого апоптоза под действием ионизирующего излучения в лимфоидных клетках больных в-клеточным хроническим лимфолейкозом // Онкология. 2008. Т. 10, № 2. С. 225-229.
[8] Желтухин А.О., Чумаков П.М. Повседневные и индуцируемые функции гена p53 // Успехи биолог. химии. 2010. Т. 50. С. 447-516.
[9] Lev Bar-Or R., Maya R., Segel L.A. et al. Generation of oscillations by the p53 — Mdm2 feedback loop: A theoretical and experimental study // PNAS. 2000. Vol. 97, No. 21. P. 11250-11255.
[10] Tiana G., Jensen M.H., Sneppen K. Time delay as a key to apoptosis induction in the p53 network // Europ. Phys. J. B. 2002. Vol. 29. P. 135-140.
[11] Ciliberto A., Novak B., Tyson J.J. Steady states and oscillations in the p53 —Mdm2 network // Cell Cycle. 2005. Vol. 4, No. 3. P. 488-493.
[12] Geva-Zatorsky N., Rosenfeld N., Itzkovitz Sh. et al. Oscillations and variability in the p53 system // Molecular Systems Biology. 2006. No. 2. doi:10.1038/msb4100068.
[13] Tiana G., Krishna S., Pigolotti S. et al. Oscillations and temporal signalling in cells // Phys. Biol. 2007. Vol. 4. P. R1-R17.
[14] Chiokarmane V., Ray A., Sauro H. M., Nadim A. A. Model for p53 dynamics triggered by damage // SIAM J. Appl. Dynamical Systems. 2007. Vol. 6, No. 1. P. 61-78.
[15] Horhat R.F., Neamtu M., Miroea G. Mathematical models and numerical simulations for the p53 —Mdm2 network // Appl. Sci. 2008. Vol. 10. P. 94-106.
[16] Horhat R.F., Neamtu M., Opris D. The qualitative analysis for a differential system of the p53 —Mdm2 interaction with delay kernel // 1st WSEAS Intern. Conf. on Biomedical Electronics And Biomedical Informatics (BEBI '08). Rhodes, Greece, 2008.
[17] Batohelor E., Mook C.S., Bhan I. et al. Recurrent initiation: A mechanism for triggering p53 pulses in response to DNA damage // Molecular Cell. 2008. Vol. 30. P. 277-289.
[18] Лихошвай В.А., Колчанов Н.А., Фадеев С.И. и др. Теория генных сетей // Системная компьютерная биология. Новосибирск: Изд-во СО РАН, 2008. С. 395-480.
[19] Гайдов Ю.А., Голубятников В.П. Одна модель генной сети, исправляющей повреждение ДНК // Математика в приложениях. Всерос. конф., приуроченная к 80-летию академика С.К. Годунова: Тез. докл. Новосибирск: Ин-т математики СО РАН, 2009.
[20] Hamada H., Tashima Y., Kisaka Y. et al. Sophisticated framework between cell cycle arrest and apoptosis induction based on p53 dynamics // PlOS ONE. 2009. Vol. 4, No. 3. e4795.
[21] Abou-Jaoude W., Chaves M., Gouze J.-L. A Theoretical Exploration of Birhythmicity in the p53 —Mdm2 Network. InRIA Research Report Model. INRIA-00523270. Version 2. 2010. 34 p.
[22] Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежёсткие задачи. М.: Мир, 1990. 512 с.
[23] Эльсгольц Л.Э., Норкин С.Б. Введение в теорию дифференциальных уравнений с отклоняющимся аргументом. М.: Наука, 1971. 296 с.
[24] Зверкина Т.С. Модификация конечноразностных методов для интегрирования обыкновенных дифференциальных уравнений с негладкими решениями // Журн. вычисл. математики и матем. физики. 1964. Т. 4, дополнение к № 4. С. 149-160.
[25] Воропаева О.Ф., Шокин Ю.И. Численное моделирование в медицине: Некоторые постановки задач и результаты расчётов // Вычисл. технологии. 2012. Т. 17, № 4. С. 29-55.
Поступила в 'редакцию 8 ноября 2012 г.
