Математическое моделирование функционирования системы биомаркеров дегенеративных заболеваний Текст научной статьи по специальности «Математика»

АПВПМ-2019

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ФУНКЦИОНИРОВАНИЯ СИСТЕМЫ БИОМАРКЕРОВ ДЕГЕНЕРАТИВНЫХ ЗАБОЛЕВАНИЙ

О, Ф, Воропаева1, К, С, Гаврилова1,2, С, Д. Сенотруеова1

1 Институт, вычислительных технологий СО РАН, 630090, Новосибирск 2 Новосибирский государственный университет, 630090, Новосибирск

УДК 519-6; 51-76

Б01: 10.24411/9999-016А-2019-10014

Работа посвящена численному исследованию функционирования наиболее важных петель положительной и отрицательной обратной и прямой связи сигнального пути белка р53. Рассматриваются четыре известные математические модели, основанные на разных биологических идеализациях процесса. Выполнен сопоставительный анализ математических моделей, даны оценки области применимости моделей. Обнаружены бифуркации Андронова-Хопфа и Неймарка-Сакера в диапазонах фазовых состояний р53 и его ингибиторов, соответствующих нормальной реакции сигнального пути р53 на стресс. В рамках принятых моделей исследованы варианты гипотетических терапевтических противораковых стратегий.

Ключевые слова: численное моделирование, уравнение с запаздыванием, онкомаркер, микроРНК, р53.

Введение

Белок р53 (супрессор опухолей) играет одну из ключевых ролей в управлении ростом, старением и отмиранием клеток. Нарушение функции этого белка приводит к развитию дегенеративных заболеваний, характеризующихся чрезмерным накоплением в организме дефектных клеток или, наоборот, массовой преждевременной клеточной смертью. Установлено, что функции р53 нарушены в клетках большинства видов рака, что позволяет рассматривать р53 в качестве перспективного диагностического маркера и терапевтической мишени. В отсутствие стрессовых воздействий р53 сохраняется на достаточно низком уровне, что предотвращает ненужную в нормальных условиях экспрессию его генов-мишеней. В случае превышения физиологически допустимого уровня повреждений ДНК р53 активируется для избавления от генетически опасных дефектных клеток. Результат активации р53 - остановка клеточного цикла и репарация ДНК при незначительном повреждении или запуск программы апоптоза (генетически запрограммированной смерти клеток) при сильном стрессовом сигнале [1,2]. В регулировании динамической реакции белка р53 участвует целый ряд положительных и отрицательных циклов обратной и прямой связи, через которые с р53 взаимодействуют, в частности, его ингибиторы — белки Мс1т2 и \¥1р1, а также многочисленные семейства микроРНК (гшТ1). Исследование роли и особенностей функционирования этих циклов в р53-сети — это ключ к более глубокому пониманию процессов дегенерации, старения, механизмов развития онкологических заболеваний и их подавления. Учитывая чрезвычайную сложность структуры сети р53, механизмов взаимодействия и функций ее участников, планирование новых лабораторных экспериментов в этой области биомедицины может быть существенным образом упрощено предварительными оценками, основанными на результатах математического моделирования.

Математические модели, предназначенные для изучения особенностей кинетики, регуляторной функции и/или терапевтического потенциала ключевых участников сети р53 при различных заболеваниях, предложены в целом ряде работ (см. обзоры в [3,4] и цитируемой там литературе). Многообразие существующих математических моделей сети р53, как правило, отражает разнонаправленность этих исследований. Однако скрупулезный подбор экспериментальных фактов, определяющих границы применимости каждой разработанной модели, является пока нерешенной проблемой.

!ЯВ.\ 978-5-901548-42-4

Одной из главных целей настоящей работы является сопоставительный анализ четырех известных математических моделей динамики сигнального пути р53, основанный на использовании лабораторных данных [1]. Предложены модификации моделей, расширяющие область их применимости за счет уточнения описания механизма отрицательной обратной связи р53-ингибитор (Мс1т2 или \У1р 1). В работе представлены результаты численного анализа решений моделей в широком диапазоне значений параметров, позволившего выявить, в частности, бифуркации рождения предельных циклов и тороидальных аттракторов как свидетельство весьма тонкой организации ответа биологической системы на стрессовые воздействия.

1 Математические модели

Модель 1. Для описания взаимосвязи р53-ингибитор-микроРНК привлекается нелинейная система дифференциальных уравнений с запаздывающими аргументами [4]. Модель составлена из балансных соотношений, отражающих вклад в совместную динамику системы р53-белок-ингибитор—микроРНК следующих трех ключевых механизмов: генерации, деградации (в рамках модели эти термины подразумевают как конститутивные спонтанные процессы, так и обусловленные влиянием каких-либо не описанных явно молекул), а также особо выделяемых процессов взаимодействия р53-ингибитор-гш11 (р53-зависимая активация генерации белка-ингибитора р53 и микроРНК, инактивация/деградация р53 под влиянием ингибитора и микроРНК-зависимое подавление генерации ингибитора):

с1Р

— = 13Р - ар1 ^(Р(г), I(г), к}) - «рр(г), ^ = ¡3!РС(Р(г - тр), I(г - тр)) - агI(г) - а1т^(I(г - п), т(г - тТ), кт),

= Рт + ртр^(Р(£ - Тт),т(г - Тт), Кр) - атт(г).

Для определения функций взаимодействия использовался подход, основанный на аппроксимациях типа Гольдбетера-Кошланда и Хилла:

Р(Р,1, К) = 2(Р + I + К - у/(Р + I + К)2 - 4Р1),

с(Р1Ю= р - ¥ (р,1,к') С(Р,Т,Кд ) = Р + кд - ^[Р,1,к{) ■

Здесь Р, I, т — уровни р53, белка-ингибитора р53 и микроРНК, соответственно; и — константы скорости генерации белка р53 и микроРНК; ар[ — константа скорости деградации р53 под влиянием ингибитора; ар, а.1, ат — константы скорости деградации р53, ингибитора и микроРНК; Др и @тр — константы скорости р53-зависимой генерации ингибитора и микроРНК; — константа скорости микроРНК-зависимого подавления белка-ингибитора; параметры К^ Кд, Кр, Кт имеют смысл констант диссоциации, которым в рамках принятой модели отводится роль регуляторов, обеспечивающих «тонкую настройку» уровней взаимосвязи образующих комплексы участников сети. Параметры тр, т^ и тт определяют время запаздывания реакции белков и микроРНК на сигналы. При анализе адекватности модели 1 в качестве белков-ингибиторов р53 рассматриваются Мс1т2 и \Vipl.

Модель 2. Для моделирования совместного функционирования двух петель отрицательной обратной связи р53-Мс1т2 и р53-\¥1р1 привлекается предложенная в [2] система нелинейных уравнений с запаздывающими аргументами. Поскольку данная модель рассматривалась в [2], в первую очередь, как инструмент для описания серии лабораторных экспериментов, в уравнения дополнительно включены слагаемые, позволяющие рассматривать воздействие нутлина — вещества, используемого для ингибирования взаимодействия р53 и Мс1т2 в терапевтических целях:

¿р. в"*

—т~ = 0р — амр-^-тт—;-г МРг — вяр-тт Р. — ар. Р.,

л ' мр кп + мп (г - тк) г ' в"* + т"* г

¿Ра _ я"* кп

= 13зр + т"'Рг - амра кп + мп (г - т„)

d¡M = Pmí + РмРа (t - тм) - aSMSM - амМ, at

dW = ftwPa (t - tw) - awW., dS = WS S

~ж = PS- aww^ttWfS-asS'

Здесь Pi, Pa — уровни неактивированного и активного белка р53; М, W — уровни Mdm2 и Wipl; S — уровень активного сигнала о повреждении; т^, тм, tw- параметры запаздывапия; N — уровень нутлина.

Модель 3. Последние данные свидетельствуют о бимодальном переключении р53, наблюдаемом при решении «судьбы» клетки в условиях стресса — при малом повреждении в системе возникают периодические колебания уровней белков, а при сильном повреждении наблюдается монотонное увеличение р53, ассоциированное с запуском программы клеточной смерти (апоптоза). Для описания этого эффекта в [5] предложена математическая модель, представляющая собой систему нелинейных дифференциальных уравнений с запаздывающими аргументами. Модель содержит 10 уравнений динамики ключевых белков сети р53 с учетом 28 биохимических реакций, которые считаются наиболее важными (в рамках принятой биологической модели) в механизме, обеспечивающем ответ сигнального пути р53 на стресс. В модели 3 рассматриваются две основные петли отрицательной обратной связи — p53-Mdm2 и рГ>:{ Wipl ATM. Модель учитывает не только конститутивные спонтанные процессы генерации и деградации р53 и двух его ключевых ингибиторов Mdm2 и Wipl, но и активацию белков через механизм фосфорилирования. Уровень этопозида рассматривается как входной параметр, инициирующий стресс. Предполагается, что этопозид как источник стресса может активировать белок ATM, что согласуется с известными экспериментальными наблюдениями. Поскольку известно, что общий уровень белка ATM практически не изменяется в течение 72 часов, но существенно изменяется его активность, в рамках принятой модели центральное место отводится именно активной форме ATM. Активированный белок ATM может фосфорилировать р53 и Mdm2, что способствует изменению активности р53 и Mdm2. Фосфорилированный р53 может активировать Wipl, который, в свою очередь, обеспечивает (в рамках принятой модели) регуляцию механизма фосфорилирования.

Модель 4. В работе [6] предложена математическая модель, которая предназначалась ее авторами для изучения функционирования механизма запуска программы р53-зависимого апоптоза в ответ на сигналы о повреждении ДНК. Модель представляет собой нелинейную систему 14 обыкновенных дифференциальных и четырех алгебраических уравнений. Она разделена на два блока. Такое разделение обусловлено принятой биологической моделью, в рамках которой ключевую роль, как и в моделях 1-3, играет петля отрицательной обратной связи p53-Mdm2 (блок 1). Во втором блоке, где непосредственно решается «судьба» клетки, р53 функционально разделен на три группы. Группа р53-киллеров активирует апоптотические гены, такие как PUMA, p53DINPl и р53А1Р1. PUMA-белок индуцирует выделение из митохондрий цитохрома, который стимулирует белковый комплекс АРА К1 для активации апоптосом, непосредственно (в рамках принятой модели) запускающих программу апоптоза. Другая форма р53 — р53-хелперы — ответственна за индукцию белка р21. Последний инициирует остановку клеточного цикла, давая время клетке для восстановления поврежденной ДНК, или участвует в запуске программы необратимого прекращения деления (старения). Третья, вспомогательная, форма р53 индуцирует Wipl, который, согласно принятой модели, сдерживает образование р53-киллеров. Отметим, что конкретная реализация биологической связи между блоками и внутри каждого из них представляет собой предмет исследований как в [6], так и в рамках настоящей работы.

В связи с тем, что уравнения моделей 1-4 представляют собой нелинейные системы дифференциальных уравнений, часть из которых содержит функции с запаздывающими аргументами, возникает проблема нарушения гладкости решения в точке t = 0 и последующего переноса разрывов производных в точки, кратные величине запаздывания. Численное решение представленных систем уравнений основано на применении идеи метода шагов, который позволяет рассматривать возникающие задачи с начальными условиями как задачи Коши для нелинейных систем ОДУ первого порядка, решаемые последовательно на интервалах времени, равных величине запаздывания. На каждом из этих интервалов функции с запаздывающими аргументами заданы в узлах расчетной сетки — они определяются из начальных условий или из численного решения, полученного на предыдущем этапе реализации метода шагов. Такой подход позволяет привлечь к решению задач широко распространенные численные методы решения задачи Коши для системы ОДУ — метод предиктор-корректор второго порядка, методы Адамса-Бэшфорта-Моултона или Гира (фактический порядок точности, как правило, не превосходил 3). Из соображений простоты численной реализации алгоритма решения систем уравнений шаг по времени h выбирался таким образом, чтобы величины запаздывания

были кратными Л. Расчеты показали, что данное упрощение не оказывает существенного влияния на характер и точность решения, если величина шага сетки является достаточно малой величиной. Окончательный выбор численного метода осуществлялся на основе анализа результатов серий методических расчетов. При решении каждой нелинейной системы дифференциальных уравнений использовалась идея метода Зейделя.

2 Результаты численного анализа

Прежде всего, представляет интерес оценка области применимости моделей 1-4, адекватность каждой из которых в [2,4-6] обосновывалась сопоставлением с существенно разными лабораторными экспериментами. В настоящей работе в качестве единого критерия для оценки всех моделей используются лабораторные данные [1], представленные в виде совокупности характерных фазовых состояний системы p53-Wipl, наблюдаемых in vitro при облучении раковых клеток с р53 дикого типа. Данные [1] показывают, что значительная часть клеток характеризуется относительно низким уровнем р53 и его отрицательного регулятора Wipl. Это следует рассматривать, по-видимому, как нормальную реакцию системы на стрессовое воздействие, приводящее к остановке клеточного цикла и репарации поврежденной ДНК. Из анализа экспериментальных данных [1] следует также, что клеткам с относительно низким сигналом р53 могут соответствовать весьма высокие уровни сигнала Wipl. Такие состояния системы p53-Wipl характеризуют, в частности, ситуацию затухания р53-зависимого апоптоза и риска рака. Точно так же наблюдаются состояния с высоким уровнем р53 при низком уровне Wipl, которые связываются с запуском программы клеточной смерти (р53-зависимого апоптоза). Известно, что именно два последних состояния позволяют рассматривать р53, его ингибиторы и р53-зависимые микроРНК как биомаркеры наиболее опасных дегенеративных заболеваний, включая рак, ишемическую болезнь сердца, болезни Альцгеймера и Паркинсона. Результаты численных экспериментов показали, что совокупность рассчитанных неподвижных точек решений модели 1 (стационарных значений уровней р53 и Wipl) располагается в фазовой плоскости в согласии с данными [1]. При этом модель 1 воспроизводит также наблюдаемые in vitro [2] колебательные режимы функционирования системы (характерные, как правило, для низкого уровня повреждения ДНК), которые регулируются петлей отрицательной обратной связи p53-Wipl. Важно, что все характерные значения численных решений модели 1 располагаются на фазовой плоскости (р53, Wipl) в том же диапазоне величин, что и экспериментальные данные. Подобная степень согласия с экспериментальными данными получена и для модели 3. В то же время, численные эксперименты показали, что в рамках моделей 2 и 4 отрицательная обратная связь р53 и его ингибиторов выражена слабо — изменение р53 вызывает сколько-нибудь значимые изменения уровня Wipl (и наоборот) лишь в относительно узком диапазоне значений уровней р53 и Wipl. Это существенно ограничивает область применимости моделей 2 и 4, фактически не позволяя их использовать как модели функционирования системы онкомаркеров p53-Wipl и/или как инструмент анализа эффективности противоопухолевых терапевтических стратегий.

Модель 1 нацелена, в первую очередь, на исследование динамики р53-зависимых микроРНК на основе минимального подхода к моделированию. Ранее авторами разработана иерархия минимальных математических моделей функционирования сети р53-ингибитор-гшТ1 с положительным прямым и обратным вариантами связи p53-miR и выполнен комплекс численных экспериментов, иллюстрирующих адекватность принятого подхода. Модель 1 — единственная модель этой иерархии, учитывающая не только прямую положительную связь р53-микроРНК, но и позволяющая «различать» разные микроРНК по уровню связи с р53 (речь идет о конкретных р53-зависимых микроРНК, для которых обратное воздействие микроРНК на р53 является несущественным либо, напротив, значимость этого воздействия чрезвычайно высока для функционирования системы в целом). Показано, в частности, что численные решения модели 1 качественно и количественно согласуются с экспериментальными данными [4] о динамике р53, а также адекватно воспроизводят изменение уровня микроРНК в ответ на активацию р53. Отдельно анализировалась адекватность математической модели взаимосвязи белка-ингибитора и микроРНК. С этой целью использовались полученные в культивируемых клетках карциномы экспериментальные данные для мРНК Mdm2 и miR-143. Численные решения модели 1 удовлетворительно согласуются с результатами лабораторных исследований, в которых было показано, что экспрессия Mdm2 отрицательно коррелирует с экспрессией miR-143. На этой основе в рамках модели 1 и ее более простых модификаций даны обобщенные оценки диагностического потенциала р53-зависимых микроРНК при онкологических и нейродегенеративных заболеваниях. Рассмотрены возможные стратегии восстановления нормального уровня р53 и связанных с ним микроРНК в целях профилактики угрозы рака. Изучены варианты противораковой терапии, связанные с одновременной гиперактивацией двух регуляторов апоптоза - р53 и микроРНК. В рамках модели 1 показана потенциально

высокая эффективность противораковой терапии, мишенью которой является белок-ингибитор р53 как основное звено петли положительной обратной связи р53—микроРНК. Выполненный комплекс исследований дает основания для того, чтобы рассматривать модель 1 как базовую при разработке более полных математических моделей сигнального пути р53.

Анализ численного решения модели 2 показал, что, в соответствии с известными представлениями, при отсутствии повреждения р53 находится в неактивном состоянии. Процесс активации р53 сопровождается стремительным ростом уровня активного р53 в ответ на сигнал и последующей регуляцией этой активности белками-ингибиторами Mdm2 и Wipl, при этом следует отметить определенную цикличность данного затухающего процесса. Результаты расчетов демонстрируют доминирующую роль Mdm2 как отрицательного регулятора р53, что, как представляется, соответствует принятой биологической модели [2]. Предложена модификация модели 2, корректирующая описание отрицательной обратной связи р53-ингибитор для согласования с наблюдениями [1]. Анализ адекватности модели 2 и ее модификации основан на сопоставлении полученных численных решений с экспериментальными данными о динамике р53 в условиях стресса, порожденного ингибированием взаимодействия р53—Mdm2 малыми молекулами нутлина. Лабораторные измерения показывают, что реакция р53 и системы р53—Mdm2 на введение нутлина является достаточно сложной и зависит от уровня нутлина, продолжительности и режима его введения в клетку. В лабораторном эксперименте постоянное введение нутлина (при постепенном наращивании его дозы на первом этапе эксперимента) достаточно быстро приводит к росту и стабилизации р53 на уровнях, значительно превышающих «безнутли-новый» в те же самые моменты времени. Аналогичный результат демонстрируют и приближенные решения модели 2 и ее модификации. При этом особенность принятой математической модели 2 состоит в том, что

N

решения эта зависимость выражена крайне слабо. Выполнен численный анализ решений исходной модели и ее модификации в широком диапазоне параметров, который дает представление о качественных свойствах решений. Обнаружены бифуркации Андронова-Хопфа и Неймарка-Сакера. Для модифицированной модели обнаружено существование двух решений - с неподвижной предельной точкой и с предельным циклом -при одинаковых значениях параметров модели.

В численных экспериментах с привлечением моделей 3 и 4 главное внимание уделялось анализу области применимости, о котором шла речь в начале данного раздела. Дополнительно с использованием модели 3 проведено численное моделирование динамики сети р53 под влиянием стрессового сигнала и терапевтического воздействия малыми молекулами этопозида. Показано, что результаты настоящих расчетов согласуются с известными представлениями (см., в частности, [5]) о бимодальном переключении пути р53 в зависимости от уровня повреждения: при малом повреждении возникают периодические колебания уровней белков, а при сильном повреждении наблюдается монотонное увеличение р53, ассоциированное с массовой гибелью клеток. Выполнено численное моделирование механизма запуска р53-зависимого апоптоза на основе математической модели 4. Рассмотрены модификации модели 4, корректирующие описание отрицательной обратной связи p53-Wipl. Выполнено численное исследование функционирования сигнального пути р53 как ключевого звена в программе клеточного ответа на кратковременное повреждение ДНК в зависимости от уровня, количества и длительности сигналов о повреждении. Показано, в частности, что при нескольких последовательных повреждениях ДНК может наблюдаться усиление процессов гибели клеток, связанное исключительно с взаимным влиянием разных, даже относительно слабых, сигналов о повреждении.

Список литературы

fl] Batchelor Е., Mock С.S., Bhan I., Loewer A., Lahav G. Recurrent initiation: A mechanism for triggering p53 pulses in response to DNA damage // Molecular Cell. 2008. V. 30, iss. 3. P. 277-289.

[2] Purvis J.E., Karhohs K.W., Mock C., Batchelor E., Loewer A., Lahav G. p53 dynamics control cell fate // Science. 2012. V. 336, iss. 6087. P. 1440-1444.

[3] Воропаева О.Ф., Сенотрусова С.Д., Шокин К).II. Дерегуляция р53-зависимых микроРНК: результаты математического моделирования // Математическая биология и биоинформатика. 2017. Т. 12, iss. 1. С. 151-175.

[4] Сенотрусова С.Д., Воропаева О.Ф. Математическое моделирование функционирования положительной связи в системе онкомаркеров р53-микроРНК // Сиб. жури, вычисл. математики. 2019. Т. 22, iss. 3.

[5] Sun Т. , Cui J. A plausible model for bimodal p53 switch in DNA damage response // FEBS Letters. 2014. V. 588. P. 815-821.

[6] Zhang Т., Brazhnik P., Tyson J.J. Exploring Mechanisms of the DNA-Damage Response: p53 Pulses and their Possible Relevance to Apoptosis // Cell Cycle. 2007. V. 6, iss. 1. P. 85-94.

Воропаева Ольга Фалалеевна — д.ф.-м.н., вед. науч. сотр.

Института вычислительных технологий;

e-mail: vorop@ict.nsc.ru; Сенотрусова Софья Дмитриевна — мл. науч. сотр.

Института вычислительных технологий;

e-mail: senotrusova.s@mail.ru; Гаврилова Ксения Сергеевна — студентка Новосибирского государственного университета;

e-mail: ksu483@yandex.ru. Дата поступления — 30 апреля 2019 г.