CC BY
9
Вычислительные технологии
Том 21, № 2, 2016
Численный анализ перехода от уравнения с запаздыванием к системе ОДУ в математической модели сети онкомаркеров
О.Ф. ВОРОПАЕВА1'2*, А. О. КОЗЛОВА2, С. Д. СЕНОТРУСОВА2
1Институт вычислительных технологий СО РАН, Новосибирск, Россия 2 Новосибирский государственный университет, Россия *Контактный e-mail: vorop@ict.nsc.ru
Выполнено численное исследование решений систем уравнений, описывающих динамику концентраций белков р53 (фактора некроза опухоли) и Мёш2 при их взаимодействии. Рассмотрены две взаимосвязанные математические модели сети р53-Мёш2 — модель на основе уравнения с запаздывающим аргументом и модель, включающая в себя систему ОДУ для описания гипотетических промежуточных стадий функционирования сети. Численно показано, что предельный переход к модели с достаточно большим числом стадий приводит к модели с запаздыванием. Исследованы конкретные условия численной реализации этого перехода.
Ключевые слова: численный анализ, уравнение с запаздыванием, онкомаркеры, р53-Мёш2-сеть.
Введение
Изучение онкомаркеров, направленное на выявление и прогнозирование исхода заболевания, является бурно развивающейся отраслью медицины и биологии. Число показателей, претендующих на роль потенциальных маркеров злокачественных опухолей, постоянно растет. Так, в клинической практике в качестве факторов онкопрогноза предлагается рассматривать белки р53 и Мёш2 и кодирующие их гены. Роль Мёш2 считается ключевой в сложной системе отрицательных регуляторов белка р53, который экспрессируется во всех клетках организма и активируется при повреждении ДНК для остановки клеточного цикла, репликации ДНК или запуска программы клеточной смерти [1]. Данные о прогностической значимости р53 и М^ш2 как онкомаркеров пока недостаточно обоснованы и поэтому требуют глубокого всестороннего изучения, в том числе средствами математического моделирования.
При построении актуальных математических моделей основное внимание уделяется учету всех существенных биологических факторов, влияющих на функционирование сети р53-М^ш2. Одним из наиболее значимых из них представляется запаздывание реакции М^ш2 на изменение состояния р53 (см., например, [1-3]). Поэтому введение
© ИВТ СО РАН, 2016
времени запаздывания т в число параметров математической модели является в достаточной мере стандартным инструментом, позволяющим достичь качественного согласия с данными лабораторных исследований (уравнения с запаздывающим аргументом включены, в частности, в известные модели работ [2, 4-7] и их модификации).
Следует отметить, что численный анализ решений нелинейных моделей, включающих уравнения с запаздывающим аргументом, представляет определенные сложности. В то же время хорошо известны результаты теоретических и численных исследований о существовании предельных переходов от системы ОДУ большой размерности к уравнению с запаздыванием и наоборот (см., например, работы [8-11] и цитируемую в них литературу). Речь идет о том, что для уравнения с запаздывающим аргументом с кусочно-непрерывной функцией начальных условий ф(^) (заданной на интервале, равном величине запаздывания, с конечным числом точек разрыва, причем все разрывы первого рода) можно сформулировать последовательность задач Коши для системы ОДУ размерности n с компонентами вектора решения (хг (t) , ..., xn (t)), такую что последовательность последних компонент этого вектора {xn (t)} при n ^ <х сходится равномерно к обобщенному решению исходного уравнения с запаздыванием. Связь указанных постановок задач можно интерпретировать как моделирование одного и того же многостадийного процесса на микро- и макроуровнях соответственно, причем под микроуровнем следует понимать, как правило, некоторые относительно быстро протекающие промежуточные стадии, обеспечивающие адекватную функциональную взаимосвязь начального и последнего звеньев рассматриваемого процесса [8]. Тогда т может иметь смысл суммарного времени протекания процесса из начального состояния (г =1) в конечное (г = п).
В настоящей работе мы будем опираться на результаты работы [8], где для автономной системы ОДУ вида
dx\ п — 1
-¡г = J Ы--Хг,
dt т
~ГТ = --1 (Хг-1 — Xi) , i = 2,...,П — 1, (1)
at т
d/^Сга 1 ,
7~ Xn— 1 UXn
at т
с нулевыми начальными условиями для компонент вектора x = (хг,х2, ...,хп) строго доказано существование пределов (1 < i < п)
lim Xi(t) ^ 0, lim xn(t) = y(t), (2)
если y(t) — решение дифференциального уравнения с запаздывающим аргументом т следующего вида:
^ = f (y(t — г)) — ву, t>r, y(t) = 0 при t<r. (3)
В уравнениях (1), (3) 9 > 0, f (хп) — достаточно гладкая функция, удовлетворяющая условию Липшица. Этот результат [8] послужил поводом для рассмотрения в настоящей работе деталей перехода от уравнения с запаздыванием к системе ОДУ без запаздывания при численном моделировании функционирования сети р53-М^т2 с применением математической модели работы [5]. При этом основной задачей исследования является не столько демонстрация основанных на результатах численных экспериментов
предельных свойств полученных решений, сколько определение оптимальной стратегии проведения больших серий расчетов с использованием вычислительной технологии замены уравнений с запаздыванием системой ОДУ при "разумном" числе уравнений в этой системе.
1. Математическая модель сети р53—М^ш2 с запаздыванием
Интересным примером математической модели взаимосвязи р53 и является нели-
нейная многопараметрическая система уравнений с запаздывающим аргументом, предложенная в [5] и наиболее детально изученная в [5, 7, 12]:
^ = ^ - а/ЫЪШ) - ЬУ1(г), (4)
^ = схд^ухЦ - т),у2Ц - т)) - С2У2СО, (5)
где взаимодействие белков определяется функциями
/(УъУ2) = 2 (у! + У2 + к/ - (У1 + У2 + к/)2 , (6)
/ ч У1 - I(У1,У2) (Г7,
9(У1,У2) = -1-V • (7)
У1 + кд - /[У1 ,У2)
Здесь у1 и у2 — концентрации белков р53 и М^ш2 соответственно; в — скорость генерации р53; а — скорость деградации р53 посредством убиквитинирования; Ь — скорость самопроизвольного распада р53; с1 — скорость производства М^ш2, в том числе за счет взаимодействия с р53; с2 — скорость деградации Mdm2; к/ — константа диссоциации комплекса p53-Mdm2; кд — константа диссоциации белка р53 и гена Mdm2. Начальные данные для системы (4)-(7) задаются в виде функций "истории":
Ук(0) = Щ(в), 0 е [-Т, 0] , к = 1, 2• (8)
Значения параметров
5 = 1 с-1, а = 3 • 10-2 с-1, Ь =10-4 с-1, с1 = 1 с-1, с2 = 10-2 с-1,
к1 = 180, кд = 28 (9)
согласованы с данными лабораторных исследований и вслед за [5] принимаются в настоящей работе как базальные, т. е. соответствующие нормальному состоянию здорового организма. Параметр т определяет время запаздывания реакции Mdm2 на изменение концентрации белка р53 и согласно наблюдениям может варьироваться от нескольких минут до десятков часов в зависимости от фактора стресса системы (как правило, в лабораторных исследованиях факторами стресса являются ультрафиолетовое и другие виды облучения клеток).
Аналитическое стационарное решение задачи (4)-(8) получено в [5] в виде системы
у(1) = С\ (а + Ъ)У1 - з у(2) = (з - Ьу1)((а + Ь)у1 + ак? - в) (10)
с2 (а + Ь)у1 - в + акд , ^2 а((а + Ь)у1 - 8)
Подробные сведения о результатах численного решения задачи (4)-(8), включая методические расчеты при достаточно малых значениях т < 120 с, исследование вопросов функционирования системы в условиях стресса при ультрадианных колебаниях в сети р53-Мёш2 и ее регуляции (перехода от состояния стресса к состоянию нормы) обсуждались в работах [13-17]. В связи с этим ниже будут представлены лишь некоторые наиболее важные для настоящего исследования факты.
В основе алгоритма решения задачи лежит широко известный метод шагов. Из соображений простоты численной реализации величина запаздывания т полагалась кратной шагу расчетной сетки (при условии, что шаг сетки является достаточно малой величиной, это ограничение не оказывает существенного влияния на характер решения). На предварительном этапе исследований [14, 17] для численного решения задачи привлекался ряд численных методов из семейств Рунге — Кутты, Адамса, Розенброка и Гира до 4-го теоретического порядка точности включительно; оценивалась роль итераций по нелинейности. При этом реальный порядок сходимости (точности) решений, полученных с применением указанных методов на последовательности сеток, составлял от 1 до 2.83. В расчетах на последовательности сеток было установлено, что наиболее предпочтительными по скорости сходимости являются явный метод Адамса и Рунге — Кутты или метод типа предиктор-корректор Адамса — Бэшфорта — Моултона (все 4-го порядка). В [17] показано также, что при условии достаточно малого шага интегрирования выбор численного метода имеет второстепенное значение. При этом простота численной реализации метода Эйлера вполне может компенсировать относительно небольшую потерю порядка точности в больших сериях расчетов. Методические расчеты проводились на последовательности сеток с разрешением от 1 до 10 000 точек на 1 с. На основании анализа результатов численных экспериментов оптимальной признана сетка с разрешением не более 100 точек на 1 с.
Расчеты выполнялись при нулевых начальных данных, что согласуется с условиями предельной теоремы [8] для задач (1)-(3). Отметим, что для задачи (4)-(8) предварительно была проведена серия численных экспериментов, в которых в качестве функций "истории" фь (к — 1, 2) использовались либо постоянные, в том числе стационарные значения концентраций, варьируемые в диапазоне от 10 % до четырех порядков, либо периодические решения данной задачи. Расчеты показали, что при достаточно больших значениях времени характер решения системы (4)-(8) полностью определяется выбором значений параметров системы, а влияние начальных условий можно считать кратковременным (более детально эти эксперименты описаны в [17]).
На рис. 1 приведены некоторые интересные для клинической практики результаты численного исследования реакции системы (4)-(8) на отклонение параметров модели от нормальных значений (более подробное описание этих данных представлено в [12, 17]). Сплошные линии образованы предельными точками или предельными циклами решений, полученными при варьировании поочередно каждого параметра модели относительно значений (9); штриховые линии соответствуют аналитическому стационарному решению (10) в базальных условиях (9) и при кд — 7.3к„ .
Расчеты проводились в достаточно большом диапазоне значений запаздывания — от нескольких секунд до полутора десятков часов, что характерно для данной биологической системы (см., например, [1, 2, 18, 19]). Получено, что неподвижные точки (фокусы) приближенных решений с достаточно высокой точностью согласуются с аналитическим решением (10). В частности, при значениях параметров (9) точка пересечения аналитических кривых у1^^ (у1) и у^^(у1) практически совпадает с неподвижной
точкой у\ав « 154.6347, у^8 ~ 81.3105 приближенного решения задачи (4)-(9). В условиях (9) система имеет решения в виде монотонных функций времени (при т = 0), обеспечивающих переход концентраций р53 и М^ш2 с нулевого на новый стационарный уровень значений, или в виде затухающих колебаний (при т > 0) с той же стационарной точкой (у \а8 ,у\а8). При отклонении параметров от значений (9) возможны бифуркации Андронова —Хопфа. В частности, в ходе численных экспериментов было получено, что к бифуркациям приводит варьирование одной из констант диссоциации — kf или кд относительно базальных значений. При этом ключевую роль играет также параметр запаздывания.
Все полученные решения модели (4)-(8) имеют достаточно ясный медико-биологический смысл. Речь идет об описании в рамках принятой модели как функционирования системы р53-М^ш2 в норме, так и о предсказании возможных опасных для организма ситуаций. В частности, наличие лишь нескольких пиков в затухающих колебаниях концентрации р53 отмечалось в известных лабораторных экспериментах по изучению поведения р53 в нормальных условиях [18]. Подобное поведение может соответствовать также реакции р53 на стресс в виде ультрафиолетового облучения живой клетки [18]. В то же время возникновение периодических колебаний в системе р53-М^ш2 считается благоприятной реакцией на воздействие гамма-излучения, позволяющей клеткам восстановить свою ДНК [1-3, 18]. Данные на рис. 1 указывают также на возможность чрезмерного накопления р53 или М^ш2. Согласно многочисленным клиническим исследованиям, эти ситуации дисбаланса возникают на фоне затухания взаимодействия в сети р53-М^ш2 и могут приводить к серьезным негативным последствиям в виде дегенерации органов в первом случае или возрастания риска онкологических заболеваний — во втором.
Рис. 1. Неподвижные точки и характерные предельные циклы решений системы (4)-(8). Сплошные линии соответствуют варьированию параметров а, Ь, 5 (линия 1), kf (линия 1'), С1, С2, кд (линия 2), к/ = 0.2Ща8 при т = 120 с (линия 3), кд = 7.3къда8 при т = 25 000 с (линия 4); пунктиры — у^ (линии 1, 5) и у)^2 (линия 2) аналитического стационарного решения при базовых параметрах (9) (линии 1, 2) и при кд = 7.3кьдаз ( линии 5, 2). Пересечение линий 1 и 2 — стационарное решение (у\а8,у2^,а8) в нормальных условиях (9), пересечение пунктирных линий 5 и 2 — стационарное решение при кд = 7.3кьда8
С учетом вышеизложенного при численном анализе технологии предельного перехода в настоящей работе будут рассмотрены два наиболее характерных варианта решений — с неподвижной точкой, соответствующее базальным значениям параметров (9), и периодическое решение (в этом случае принималось kf — 1.8, а остальные значения параметров полагались равными (9)). Из соображений простоты численной реализации в ходе расчетов полагалось т — 120 с.
2. Математическая модель сети р53—М^ш2 без запаздывания
Пусть в согласии с (1), (2) исходная система уравнений (4), (5) может быть заменена следующей системой ОДУ:
£ —.-. д.,,.) -
¿х1 п — 1
—¡Г = С19 (У1,хп)--х1,
аъ т
^ГТ = --1 (Хг-1 —Хг), 1 — 2,...,П — 1, (11)
йхп п — 1
7~ хп— 1 С-2ХП1
где х1 (¿) , ..., хп (¿) — дополнительные переменные, имеющие смысл гипотетических промежуточных стадий в функционировании сети р53-Мёш2. Предполагается, что первая и последняя компоненты решения задачи Коши для системы ОДУ (11) при начальных условиях у1(0) — 0, хД0) — 0 (г — 1,..., п) сходятся к решению у1(Ь) и у2(Ь) исходной системы уравнений (4), (5) при фк(в) — 0: |у1(1) — у1(^)| ^ 0 и |хп(£) — у2(Ь)1 ^ 0 равномерно при п ^ то. Функции /ид задаются формулами (6), (7). В ходе численных экспериментов решение задачи (11) сопоставляется с известным решением исходной системы с запаздыванием, полученным на достаточно подробной расчетной сетке и принятым в качестве "точного". Решение задачи Коши (11) проводилось методом Эйлера, дополнительно привлекались методы Рунге — Кутты четвертого порядка и алгоритм первого порядка, основанный на идее метода Зейделя, в котором при вычислении правых частей уравнений участвуют уже полученные на новом временном слое значения компонент решений, а остальные значения берутся с предыдущего временного слоя (ниже — метод Эйлера — Зейделя).
На рис. 2 показана сходимость (с ростом размерности системы (11)) компонент у1 и хп численного решения задачи (11) к соответствующему решению задачи (4)-(7) при нулевых начальных условиях и базальных значениях параметров (9). Согласно расчетам уже при п — 4 система (11) имеет то же стационарное решение (у\аз, у2аз), что и система с запаздыванием. При этом основные погрешности сосредоточены на временном интервале (0, г], где по условию задачи наиболее наглядно проявляется запаздывание реакции Мёш2 на изменение р53. При достаточно больших п компоненты решения системы (11) весьма близки на всем рассмотренном интервале времени к решению задачи в исходной постановке с запаздыванием.
На рис. 3 продемонстрирована сходимость последовательности компонент у1 (¿) и хп(к) решений системы (11), начиная с п — 4, к решению у1(^, у2(¿) системы с запаздыванием для одного из характерных вариантов периодического решения. Как и в случае
Рис. 2. Изменение во времени компонент решения системы (11) при п = 4 (пунктирные линии) и решения системы с запаздыванием (сплошные линии) в базальных условиях: 1 — р53, 2 — М^ш2 (а); сходимость последовательности решений системы (11) при п = 4 (пунктирная линия), 8, 16, 32 к решению системы с запаздыванием (линия 1) в условиях (9) (б)
Рис. 3. Сходимость последовательности решений ух (£), хп(£) системы (11) без запаздывания при п = 4 (пунктирные линии 2), 8, 16, 32 к решению системы (4)-(8) с запаздыванием у\(£), у2(£) (линии 1) при kf = 1.8 для белков р53 (а) и М^ш2 (б)
решения с предельной точкой, на временном интервале (0, г] с ростом п наблюдается постепенное "выполаживание" хп (¿), соответствующее начальному условию для функции у2 исходной задачи с запаздыванием. При этом период колебаний достаточно хорошо согласуется с величиной Т ^ 2т, полученной для решения задачи с запаздыванием, так что с ростом п наиболее существенной представляется коррекция решения по амплитуде колебаний. Фазовые портреты этих решений представлены на рис. 4. Можно видеть, что в случае kf — 1.8 все численные решения системы (11) при п — 4, 8, 16, ... образуют вложенную последовательность предельных циклов, сходящихся при достаточно больших п к предельному циклу решения системы с запаздыванием (линия 1 на рис. 4).
40 60 80 р53
Рис. 4. Предельные циклы периодических решений при kf = 1.8: сходимость последовательности решений системы без запаздывания при п = 4 (пунктирная линия 2), 8, 16, 32 к решению системы с запаздыванием (линия 1)
Для более детального численного анализа процесса сходимости исследовались данные об изменении погрешности вычисления компонент вектора решений и хп(Ь) по отношению к решению у1{Ъ), у2(¿), принятому в качестве "точного", в зависимости от размерности системы (11). Численные эксперименты проводились для 4 < п < 220 000. Погрешность еп определялась в матричной чебышевской норме и отдельно для у^) и хп(Ъ) в согласованной векторной норме [20]. Согласно теоретическим оценкам [9], скорость сходимости решений в предельном переходе от задачи (1) при нулевых начальных условиях к задаче (3) определяется асимптотическим законом п-0'5. В настоящих расчетах на рассмотренном интервале значений п получено несколько более быстрое убывание еп: от п-0'7 до п-0'9 в зависимости от принятой стратегии проведения численных экспериментов (при этом следует отметить более быструю равномерную сходимость компоненты решения у^) к у\^), для которой в расчетах получен степенной закон п-1).
На рис. 5 представлены некоторые данные об изменении еп, полученные методами Эйлера — Зейделя при значениях параметров (9); в каждой серии расчетов шаг по времени полагался фиксированным. Видно, что при достаточно умеренных п < 103, представляющих наибольший интерес с точки зрения практической реализации технологии перехода, поведение еп соответствует указанным степенным закономерностям, а при больших п имеют место существенные отклонения, если шаг расчетной сетки по времени не корректируется при увеличении п (ниже эти вычислительные эффекты будут обсуждаться более подробно). В случае периодического решения процесс сходимости характеризуется аналогичными закономерностями. Важно отметить, что все степенные зависимости, приближающие (с точностью до константы) функцию изменения погрешности решения системы ОДУ по отношению к решению исходной системы с запаздыванием в зависимости от п, определены численно и, разумеется, не претендуют на роль точных асимптотических оценок.
102
8,
\
ХО"3 11111N_11м 11м!_11м 11м!_11м 11м!_11м 11111
104 п ю:
10
10'
10
Рис. 5. Изменение нормы еп погрешности предельного перехода в базальных условиях (9) в зависимости от размерности системы ОДУ: 1 — метод Эйлера — Зейделя, к = 0.01 с; 2 — метод Эйлера — Зейделя, к = 0.001 с; 3-5 — метод Эйлера, к = 0.01 с; сплошные линии — погрешность в матричной равномерной норме; штрихпунктирные 4 и 5 — соответственно погрешности вычисления ух (£) и хп (£); пунктирные — приближенные степенные зависимости
В ходе расчетов анализировалась также погрешность численных решений на последовательности конечно-разностных сеток (при фиксированном п). Было установлено, в частности, что при п = 1200 и к = 0.01 с она достигает значения 0.0196, что составляет не более 0.0125 % по отношению к характерному значению решения при параметрах (9), а фактический порядок численных методов на этих решениях весьма близок к единице. Эти данные являются свидетельством того, что для больших серий численных экспериментов такие значения параметров перехода весьма близки к оптимальным.
Функции Хг (£), описывающие гипотетические промежуточные стадии процесса, показаны на рис. 6 при различных значениях п для решения с предельной точкой и для периодического решения соответственно. Из анализа данных о поведении Х{ (£) следует lxi(¿)| — 0 (2 < г < п — 1), что согласуется с (2) (см. также [8]). Сами эти функции имеют тот же характер зависимости от времени, что и соответствующее решение исходной задачи (4)-(8): выходят на постоянное значение в случае базальных значений параметров модели или на режим периодических колебаний при kf = 1.8. При этом для каждого фиксированного п с возрастанием номера функции г можно наблюдать коррекцию Х{ (£) в виде "фазового" сдвига по временной шкале на величину, соответствующую значению параметра запаздывания г. В случае периодического решения можно отметить также незначительную коррекцию х^(£) по амплитуде, наблюдаемую с возрастанием номера функции г только при весьма малых п.
Дополнительно изучался вопрос о зависимости линий нейтральности, полученных в результате численного анализа решений для двух рассмотренных постановок задачи, от размерности системы (11). В частности, на рис. 7 представлены линии нейтральности, разделяющие плоскость параметров (т,kf) на область, в которой стационарное решение асимптотически устойчиво (над соответствующей линией нейтральности), и область периодических решений (под линией).
На рис. 8 линии нейтральности в плоскости параметров (т, кд) отражают более сложную зависимость решений от указанных параметров, когда при достаточно больших т решение характеризуется не одной, а двумя бифуркациями (область периодических решений располагается "внутри" линий нейтральности). Данные численного анализа демонстрируют сходимость последовательности линий нейтральности, полученных для решений системы (11), к линиям нейтральности, соответствующим задаче (4)-(8), при нулевых начальных условиях.
Рис. 6. Характерные функции (£) при различных значениях п: а — в случае базальных значений параметров (пунктиры — г = 2, сплошные жирные линии — г = п — 1); б —в случае kf = 1.8 (линии 1-3 соответствуют расчетам при п = 8, 16, 64; пунктиры — г = 2, сплошные линии — г = п — 1
Рис. 7. Линии нейтральности в плоскости параметров (т,kf) для последовательности решений у1 (£), хп(£) при п = 4 (пунктир), 8, 16, 32, 64, 128, 256 и для решения системы с запаздыванием у1 (£), у2 (£) (линия 1)
Рис. 8. Линии нейтральности в плоскости параметров (г, кд) для последовательности решений у1 (£), хп(£) при п = 8 (пунктир), 16, 32, 64, 128 и для решения системы с запаздыванием у1 (£), у2 (£) (линия 1)
Анализировалась зависимость погрешности перехода от уравнения с запаздыванием к системе ОДУ от используемого численного метода. Выше отмечалась существенная роль метода при п > 103. При этом следует отметить достаточно слабый разброс значений погрешностей, полученных разными методами при фиксированных п для умеренных значений этого параметра. Так, в случае решения с предельной точкой (для параметров (9)) величина погрешности еп в расчетах методом Эйлера составляла 0.31167 • 101 (п = 32) и 0.18169 • 101 (п = 64), а методом Рунге — Кутты получены значения еп = 0.31347 • 101 (п = 32) и еп = 0.18416 • 101 (п = 64).
Для оценки влияния параметров конечно-разностной сетки на погрешность предельного перехода вновь обратимся к данным, представленным на рис. 5. Отметим, что в рамках каждого метода при умеренных значениях п получена достаточно слабая зависимость величины £п от шага по времени к (для принятого на основании анализа результатов расчетов исходной задачи диапазона значений к). Например, в случае ба-зального решения при п = 200 метод Эйлера — Зейделя дает погрешность еп = 0.878 346 на основной сетке с шагом к = 0.01 с и еп = 0.849 642 на сетке с шагом к = 0.001с; в случае периодического решения (kf = 1.8) для явного метода Эйлера при к = 0.01 с получено £п= 0.55 151 952 • 102, а при к = 0.00125с — еп = 0.54 510 549 • 102 для системы (11) при п = 32. Вместе с тем в расчетах с шагом к = 0.01 с (отметим, что эта сетка обеспечивает высокое разрешение — 24 000 узлов на 1 характерный период колебаний концентраций у1 и у2) при п ^ 104 могут возникать вычислительные трудности, которые выражаются в существенном замедлении убывания еп и свидетельствуют о доминировании погрешностей численного алгоритма (метода и сетки).
Так, в расчетах методом Эйлера — Зейделя при переходе от п = 54 000 к п = 220 000 погрешность уменьшилась менее чем на 8 %. Аналогичные трудности наблюдались и в расчетах другими методами. Эти ситуации указывают на необходимость корректировки расчетной сетки для исключения дисбаланса в погрешностях предельного перехода и численного алгоритма (подобные меры предпринимались, например, в расчетах авторов [10]). Эта коррекция, согласно результатам наших численных экспериментов, предполагает выполнение некоторого соотношения для п, т и к. В частности, из анализа результатов расчетов явным методом Эйлера следует, что при весьма больших п эффективный счет возможен, если п < т/к (в данном случае это условие является также условием устойчивости); подобное соотношение возникает и в расчетах методом Эйлера — Зейделя. Очевидно при этом, что существенное изменение шага сетки может усилить и требования к решению исходной задачи, принятому в качестве "точного". Численный анализ этих и других результатов расчетов показывает, что для эффективной реализации рассматриваемой в настоящей работе вычислительной технологии главную роль в погрешности решения системы (11) по отношению к решению системы (4)-(8), по-видимому, должна играть размерность системы ОДУ. Более детальное изучение данного вопроса является предметом дальнейших исследований.
Заключение
В ходе численных экспериментов показана возможность перехода от уравнения с запаздывающим аргументом к системе ОДУ в нелинейной математической модели функционирования системы белков р53 и М^ш2, участвующих в определении "судьбы" клетки. Изучены детали основанного на результатах численных экспериментов предельного пе-
рехода для двух характерных ситуаций — в норме и в случае возникновения периодических колебаний в системе p53-Mdm2 в ответ на стресс.
Получена равномерная сходимость первой и последней компонент решения задачи Коши для системы ОДУ достаточно большой размерности, описывающей гипотетические промежуточные стадии функционирования сети p53-Mdm2, к решению соответствующей системы нелинейных уравнений с запаздывающим аргументом при нулевых начальных условиях.
Показано, что для рассмотренных математических моделей скорость сходимости превосходит известные теоретические оценки, однако для достижения погрешности предельного перехода, сопоставимой с погрешностью используемых численных методов решения исходной задачи с запаздыванием, размерность системы ОДУ без запаздывания должна быть достаточно велика. Дальнейшее ускорение сходимости, по-видимому, возможно на основе коррекции вида уравнений для промежуточных стадий. Обсуждается вопрос о выборе оптимальной стратегии проведения численных экспериментов, обеспечивающих достаточно высокую точность решений задачи при умеренных значениях размерности системы ОДУ и шага по времени.
Благодарности. Авторы благодарят д. ф.-м. н. профессора С.И. Фадеева за постановку задачи (11) и полезные обсуждения.
Список литературы / References
[1] Lane, D., Levine, A. p53 research: The past thirty years and the next thirty years // Cold Spring Harb Perspect Biol. 2010. Vol. 2, No. 12. a000893.
[2] Geva-Zatorsky, N., Rosenfeld, N., Itzkovitz, S., Milo, R., Sigal, A., Dekel, E., Yarnitzky, T., Liron, Y., Polak, P., Lahav, G., Alon, U. Oscillations and variability in the p53 system // Molecular Systems Biology. 2006. Vol. 2, No. 1. P. 1-13.
[3] Toettcher, J.E., Mock, C., Batchelor, E., Loewer, A., Lahav, G. A synthetic-natural hybrid oscillator in human cells // Proc. of the Nat. Acad. of Sci. of the United States of America. 2010. Vol. 107, No. 39. P. 17047-17052.
[4] Mihalas, G.I., Simon, Z., Balea, G., Popa, E. Possible oscillatory behavior in p53-Mdm2 interaction computer simulation // J. of Biological Systems. 2000. Vol. 8, No. 1. P. 21-29.
[5] Tiana, G., Jensen, M.H., Sneppen, K. Time delay as a key to apoptosis induction in the p53 network // The European Physical Journal B-Condensed Matter and Complex Systems. 2002. Vol. 29. P. 135-140.
[6] Ma, L., Wagner, J., Rice, J., Hu, W., Levine, A., Stolovitzky, G. A plausible model for the digital response of p53 to DNA damage // Proc. of the Nat. Acad. of Sci. of the United States of America. 2005. Vol. 102, No. 4. P. 14266-14271.
[7] Horhat, R.F., Neamtu, M., Mircea, G. Mathematical models and numerical simulations for the p53 —Mdm2 network // Appl. Sci. 2008. Vol. 10. P. 94-106.
[8] Лихошвай В.А., Фадеев С.И., Демиденко Г.В., Матушкин Ю.Г. Моделирование уравнением с запаздывающим аргументом многостадийного синтеза без ветвления // Сиб. журн. индустр. матем. 2004. Т. 7, № 1. С. 73-94.
Likhoshvai, V.A., Fadeev, S.I., Demidenko, G.V., Matushkin, Yu.G. Modeling nonbranching multistage synthesis by an equation with retarded argument // Sib. Zh. Ind. Mat. 2004. Vol. 7, No. 1. P. 73-94. (In Russ.)
[9] Мельник И.А. Об одной нелинейной системе дифференциальных уравнений, моделирующей многостадийный синтез вещества // Вест. ТГУ. Естеств. и техн. науки. 2011. Т. 16, № 5. С. 1254-1259.
Melnik, I.A. On the nonlinear system of differential equations modelling of multiphase synthesis of a substance // Vestnik TSU. Natural and Engineering Sciences. 2011. Vol. 16, No. 5. P. 1254-1259. (In Russ.)
[10] Лихошвай В.А., Фадеев С.И., Ш^токало Д.Н. Об исследовании нелинейных моделей многостадийного синтеза вещества. Новосибирск, 2010. 37 с. (Препр. / Ин-т математики им. С.Л. Соболева СО РАН; № 246).
Likhoshvay, V.A., Fadeev, S.I., Shtokalo, D.N. On the study of nonlinear models for multiphase synthesis of a substance. Novosibirsk, 2010. 37 p. (Preprint / Institut Matematiki im. Soboleva; No. 246). (In Russ.)
[11] Демиденко Г.В. Системы дифференциальных уравнений высокой размерности и уравнения с запаздывающим аргументом // Сиб. матем. журн. 2012. Т. 53, № 6. С. 1274-1282. Demidenko, G.V. Systems of differential equations of higher dimension and delay equations // Sibirskii Matematicheskii Zhurnal. 2012. Vol. 53, No. 6. P. 1021-1028.
[12] Воропаева О.Ф., Шокин Ю.И., Непомнящих Л.М., Сенчукова С.Р. Математическое моделирование функционирования и регуляции биологической системы p53-Mdm2. М.: Изд-во РАМН, 2014.
Voropaeva, O.F., Shokin, Yu.I., Nepomnyashchikh, L.M., Senchukova, S.R.
Mathematical model of function and regulation of the p53-Mdm2 biological system. Moscow: Publishing House RAMS, 2014. (In Russ.)
[13] Воропаева О.Ф., Шокин Ю.И. Численное моделирование в медицине: Некоторые постановки задач и результаты расчетов // Вычисл. технологии. 2012. Т. 17, № 4. С. 29-55. Voropaeva, O.F., Shokin, Y.I. Numerical simulation in medicine: Formulations of the problems and some results of calculations // Computat. Technologies. 2012. Vol. 17, No. 4. P. 29-55. (In Russ.)
[14] Воропаева О.Ф., Шокин Ю.И. Численное моделирование обратной связи p53-Mdm2 в биологическом процессе апоптоза // Вычисл. технологии. 2012. Т. 17, № 6. С. 47-63. Voropaeva, O.F., Shokin, Y.I. Numerical simulation of feedback p53-Mdm2 in biological process of apoptosis // Computat. Technologies. 2012. Vol. 17, No. 6. P. 47-63. (In Russ.)
[15] Воропаева О.Ф., Шокин Ю.И., Непомнящих Л.М., Сенчукова С.Р. Математическое моделирование функционирования системы белков p53-Mdm2 // Бюл. экспер. биологии и медицины. 2014. Т. 157, № 2. С. 261-264.
Voropaeva, O.F., Shokin, Yu.I., Nepomnyashchikh, L.M., Senchukova, S.R.
Mathematical modeling of functioning of the p53-Mdm2 protein system // Bulletin of Experimental Biology and Medicine. 2014. Vol. 157, No. 2. P. 291-294.
[16] Воропаева О.Ф., Шокин Ю.И., Непомнящих Л.М., Сенчукова С.Р. Математическое моделирование регуляции биологической системы p53-Mdm2 // Бюл. экспер. биологии и медицины. 2014. Т. 157, № 4. С. 539-542.
Voropaeva, O.F., Shokin, Yu.I., Nepomnyashchikh, L.M., Senchukova, S.R.
Mathematical modeling of p53-Mdm2 protein biological system regulation // Bulletin of Experimental Biology and Medicine. 2014. Vol. 157, No. 4. P. 535-538.
[17] Воропаева О.Ф., Сенчукова С.Р., Бродт К.В., Гарбузов К.Е., Мельниченко А.В., Старикова А.А. Численное моделирование ультрадианных колебаний в биологической системе p53-Mdm2 в условиях стресса // Матем. моделирование. 2014. Т. 26, № 11. С. 105-122.
Voropaeva, O.F., Senchukova, S.R., Brodt, K.V., Garbuzov, K.E., Mel'nichen-ko, A.V., Starikova, A.A. Numerical simulation of ultradian oscillations in p53-Mdm2 network under stress conditions // Mathematical Models and Computer Simulations. 2015. Vol. 7, No. 3. P. 281-293.
[18] Loewer, A., Batchelor, E., Gaglia, G., Lahav, G. Basal dynamics of p53 reveals transcriptionally attenuated pulses in cycling cells // Cell. 2010. Vol. 142, No. 1. P. 89-100.
[19] Batchelor, E., Loewer, A., Mock, C., Lahav, G. Stimulus-dependent dynamics of p53 in single cells // Molecular Systems Biology. 2011. Vol. 7, Article number 488. 8 p.
[20] Бабенко К.И. Основы численного анализа. М.: Наука, 1986. 744 с.
Babenko, K.I. Fundamentals of numerical analysis. Moscow: Nauka, 1986. 744 p. (In Russ.)
Поступила в редакцию 5 октября 2015 г., с доработки — 3 декабря 2015 г.
Numerical analysis of the transition from the equation with retarded argument to the ODE system in a mathematical model of the tumor markers network
VOROPAEVA, OLGA F.1,2*, KOZLOVA, ALINA O.2, SENOTRUSOVA, SOFYA D.2
institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia 2 Novosibirsk State University, Novosibirsk, 630090, Russia * Corresponding author: Voropaeva, Olga F., e-mail: vorop@ict.nsc.ru
The p53 protein (tumor necrosis factor) involved in many life and death processes is expressed in all cells of the organism including the formation of tumors and aging. Mdm2 protein is considered to be the key negative p53 regulator. The investigation of the mechanism of p53 and Mdm2 interaction is paramount for developing new approaches to cancer treatment and determining the prevention strategy for many diseases.
This paper is devoted to a development of the effective numerical technology for solution of a system of equations describing the dynamics of the p53-Mdm2 network. We consider two interrelated mathematical models of the p53-Mdm2 network. The first model of the proteins concentrations dynamics includes the system of two nonlinear differential equations with the retarded argument. The second model describes hypothetical stages of process and uses the simplest ODE system of sufficiently higher dimension. We numerically show that in the transition to the limit in which the second model has sufficiently many stages, we obtain model based equation with retarded argument. We examined the specific conditions of the numerical realization of this transition.
Keywords: numerical analysis, equation with retarded argument, tumor marker, p53-Mdm2 network.
Received 5 October 2015 Received in revised form 3 December 2015
© ICT SB RAS, 2016
