ФАЗОВЫЕ ПОРТРЕТЫ МОДЕЛЕЙ ДВУХ ГЕННЫХ СЕТЕЙ Текст научной статьи по специальности «Математика»

Математические заметки СВФУ Январь—март, 2021. Том 28, № 1

УДК 514.745.82

ФАЗОВЫЕ ПОРТРЕТЫ МОДЕЛЕЙ ДВУХ ГЕННЫХ СЕТЕЙ В. П. Голубятников, Н. Е. Кириллова

Аннотация. Построены математические модели функционирования двух малокомпонентных генных сетей, участвующих в регуляции суточных ритмов в живых организмах посредством комбинаций положительных и отрицательных обратных связей между их компонентами. Доказано, что фазовые портреты этих моделей, представленных в виде нелинейных динамических систем биохимической кинетики, содержат в точности по одной стационарной точке. Установлено, что в обоих случаях при всех значениях параметров указанных динамических систем матрицы их линеаризаций в окрестностях стационарных точек имеют либо отрицательные собственные значения, либо собственные значения с отрицательными вещественными частями. Это означает, что стационарные точки рассмотренных систем устойчивы. Построены инвариантные окрестности этих точек, дано описание поведения траекторий обеих систем и биологическая интерпретация полученных результатов.

Б01: 10.25587/8УРи.2021.68.70.001

Ключевые слова: модель генной сети, нелинейная динамическая система, фазовый портрет, стационарная точка, устойчивость, критерий Вышнеградского.

Введение

Рассматривается математическая модель одной трехкомпонентной генной сети, которая участвует в регуляции циркадных (суточных) ритмов активности процессов метаболизма многих видов живых организмов. В работах [1-3] приводится описание схем соответствующих биохимических процессов, которые регулируются комбинациями положительных и отрицательных обратных связей между компонентами аналогичных генных сетей.

На рис. 1 приведена схема одной такой генной сети, компоненты которой для краткости обозначены через X, YиZ. Функционирование этой сети происходит следующим образом: вещество X действует на вещества Y и Z посредством положительных обратных связей — прямые стрелки (чем выше концентрация вещества X, тем выше скорость синтеза обоих веществ Y и Z). Вещества Y и Z ингибируют вещество X (скорость синтеза вещества X уменьшается с ростом концентраций каждого из веществ Y и Z), это отрицательные обратные связи, обозначенные на схеме прерывистыми «кривыми линиями» с Т-образными окончаниями. Подобные обозначения используются в литературе при описании

Работа выполнена в рамках государственного задания ИМ СО РАН (проект № 0314— 2019-0011) и поддержана РФФИ (грант 20-31-90011).

© 2021 Голубятников В. П., Кириллова Н. Е.

Рис. 1. Схема моделируемой трехкомпонентной генной сети.

схем функционирования разнообразных природных генных сетей (см., например, [3-5]).

Следуя этим публикациям, в качестве математической модели рассматриваемой генной сети будем изучать нелинейную динамическую систему с гладкими правыми частями

^ = Li(y) + L2{z) - k]_x, ^ = Til» - k2y, ^ = T2(x) - k3z, (1)

в которой x(t), y(t), z(t) — концентрации указанных компонент X, Y и Z соответственно, гладкие класса C1 монотонно убывающие неотрицательные функции

L1(y) и L2(z) моделируют отрицательные обратные связи Y---|X и Z---IX,

как на рис. 1. Монотонно возрастающие ограниченные сверху неотрицательные гладкие функции Г1(ж), Г2(х) моделируют положительные обратные связи X ^ Y, X ^ Z.

Как ив [1,4], положительные коэффициенты k1, k2, k3 характеризуют скорости естественного разложения компонент моделируемой генной сети. Будем называть эту генную сеть «первой», далее рассмотрим еще одну подобную конструкцию. Обозначим через A максимальное значение функции . поскольку слагаемые монотонно убывают, А = ^ xjepe3 ^^ и (jбудем обозначать точные верхние грани значений монотонно возрастающих функций Г1(ж) и Г2(х) соответственно.

Во всех уравнениях системы (1) и системы (5), рассмотренной в заключительном разделе этой работы, все переменные предполагаются неотрицательными. Фазовые портреты подобных динамических систем и качественные методы анализа поведения их траекторий изучались ранее при моделировании биохимических процессов в ряде других природных генных сетей [6-8].

1. Стационарные точки системы (1)

Положения равновесия системы (1), т. е. точки, в которых все производные Ж' if' if Равны нулю, находятся из уравнений

dx

— = Li(y) + L2(z) — k\x = 0, или k\x = L\(y) + L2(z),

dy r i \ i n Г1(х)

— = Щж) - k2y = 0, или y=—-, (2)

dt k2

йг -г I \ ; п Г2(х) /оч — = 12\х) — к^г = I), или 2 = —--. (Л)

аЬ к3

Отсюда следует, что

Поскольку композиция монотонно возрастающей и монотонно убывающей функции является монотонно убывающей функцией, в уравнении (4) обе композиции ¿1 (Г1Р) и монотонно убывают с ростом ж, как и их сумма.

Следовательно, графики монотонно возрастающей функции / = кхх и монотонно убывающей функции (Щ^) +¿2 (н^) имеют в точности одну точку пересечения; пусть х = х0 — ее абсцисса.

По этой координате х0 находим единственным образом из уравнений (2), (3)

/1\ 0 Г1 (х0) 0 Г2(ж0)

остальные координаты стационарной точки системы (1): у = ^ ', г = ^ '. В точке (х0,у0,г0), лежащей во внутренности положительного октанта, левые части всех уравнений системы (1) обращаются в нуль. Значит, имеет место

Лемма 1. Система (1) имеет в точности одну стационарную точку Б0 = (х0, у0, г0) при любых параметрах системы.

Рассмотрим параллелепипед ( = [0,А] х [0, В] х [0, С]. Так же, как и в [4,8,9], где рассматривался ряд аналогичных моделей генных сетей, проверяется, что стационарная точка Б0 лежит в его внутренности и что ( является положительно инвариантной областью системы (1), т. е. с ростом Ь траектории точек этого параллелепипеда не выходят за его границы.

Доказательство этих фактов состоит в проверке знаков производных йх/йЬ, йу/йЬ, йг/йЬ на гранях параллелепипеда (. Траектории всех точек его границы с ростом Ь входят в его внутренность и остаются в ней.

Хорошо известно, что в окрестностях гиперболических стационарных точек поведение траекторий нелинейных динамических систем полностью определяется собственными числами матриц линеаризации этих систем в таких точках.

Отметим, что модели генных сетей с более сложными, чем на рис. 1, 2, схемами взаимодействия компонент могут содержать несколько стационарных точек, в том числе и устойчивых, и поведение траекторий соответствующих динамических систем становится намного сложнее, чем в рассматриваемых здесь случаях. В частности, в таких моделях могут наблюдаться явления мультиста-бильности — различные траектории могут притягиваться к различным устойчивым точкам (см., например, [4, 7, 8] и цитированную там литературу).

2. Устойчивость стационарной точки £0

Матрица линеаризации системы (1) в окрестности точки Б0 имеет вид

-кх —рх —Р2 ' М0 = | дх -к2 0

42 0 -кз

здесь —р 1 = -щ^, ~Р2 = все производные этих функций вычисляются в стационарной точке, коэффициенты рх, р2 положительны, так как функции Ьх, Ь2 монотонно убывают.

Коэффициенты дх = , 92 = тоже положительны, поскольку функции Гх, Г2 монотонно возрастают.

Характеристический многочлен матрицы М0 имеет вид

Р(А) = -(А + кх)(А + к2)(А + кз) - (А + к2>242 - (А + кз)рх4х-

Подстановкой А = -к2 и А = -к3 можно убедиться в том, что при к2 = к3 выполняется равенство Р(-к2) = Р(-к3) = 0. Если же к2 = к3, то Р(-к2) и Р(-к3) имеют противоположные знаки, и, значит, отрезок между -к2 и -к3 содержит корень характеристического многочлена. Обозначим минимальный отрицательный корень многочлена Р(А) через Ах. Для проверки знаков вещественных частей остальных корней этого многочлена воспользуемся следующим критерием Вышнеградского [10] (см. также [11]).

У кубического уравнения ах3 + Ьх2 + сх + й = 0 с положительными коэффициентами все корни имеют отрицательные вещественные части или сами отрицательны тогда и только тогда, когда ай < Ьс.

Рассмотрим многочлен - Р(А) = А3 + А2(кх + к2 + к3) + А(кхк2 + кхк3 + к2к3 + Р242 + Рх4х)

+ кх к2к3 + к2Р242 + к3Рх4х = 0. Неравенство ай < Ьс в нашем случае принимает вид

кхк2к3 + к2Р242 + к3Рх4х < (кх + к2 + к3)(кхк2 + кхк3 + к2к3 + Р242 + Рх4х) и очевидным образом выполняется при любых положительных параметрах динамической системы (1) и ее линеаризации в окрестности стационарной точки. Значит, условие Вышнеградского для матрицы М0 этой линеаризации выполнено, все ее собственные числа либо отрицательны, либо имеют отрицательные вещественные части, и стационарная точка $0 является гиперболической. Тем самым установлена

Теорема. Стационарная точка 50 системы (1) устойчива.

Итак, траектории динамической системы (1) в инвариантной области ( ведут себя следующим образом.

(a) Если матрица М0 имеет всего один вещественный корень Ах, то траектории всех точек области ( , за исключением двух — приближающихся к стационарной точке $0 в направлении, соответствующем собственному значению Ах,— притягиваются стационарной точкой 50 «по спиралям».

(b) Если же все корни многочлена Р(А) вещественны, например, в случае, когда произведения Р242 и Рх4х достаточно малы, то траектории всех точек инвариантной области ( притягиваются к стационарной точке $0 по соответствующим направлениям.

Заключение

Таким образом, рассмотренная трехкомпонентная модель (1) генной сети в изоляции от других генных сетей не имеет осциллирующих траекторий. Для моделирования циркадных ритмов в живом организме соответствующую динамическую систему следует расширить. Подобные явления возникновения ос-цилляций в динамических системах, моделирующих комплексы всего из одной пары малокомпонентных генных сетей (двухкомпонентных и трехкомпонент-ных), хорошо известны (см. [12,13] соответственно, см. также [8]). Отметим, что для подобных динамических систем размерности 5 и выше аналогичные параллелепипеду Q инвариантные области могут содержать циклы, в том числе и по несколько циклов, а также инвариантные поверхности (см., например, [9,14]).

Рис. 2. Схема второй трехкомпонентной генной сети.

Дословно такими же рассуждениями, как и выше, устанавливаются единственность и устойчивость стационарной точки в модели генной сети, описываемой схемой, изображенной на рис. 2. Будем называть эту модель «второй». Здесь, как и для первой модели, положительные обратные связи обозначены стрелками, а отрицательные — прерывистыми линиями с Т-образными окончаниями. Нелинейная динамическая система, моделирующая вторую генную сеть, «противоположную» первой, представлена в следующей форме:

-^ = ц{у) + 12^)-к1х, — = ¿1{х) - к2у, = £2{х)-к3г. (5)

Матрица линеаризации этой динамической системы в окрестности ее стационарной точки принимает вид

(6)

Как и для первой модели, положительные коэффициенты а± и а2 обозначают вычисленные в стационарной точке производные монотонно возрастающих положительных ограниченных функций 71 и 72, описывающих положительные обратные связи; отрицательные коэффициенты (—Ь1) и (—Ь2) — вычисленные в той же точке производные монотонно убывающих положительных функций 11 и 12, которые описывают отрицательные обратные связи. Аналогичным образом

-к1 а1 а2

-Ь1 -к2 0

-Ь2 0 -кз

в положительном октанте пространства переменных (ж, y, z) для динамической системы (5) конструируется инвариантный параллелепипед Q и описывается поведение ее траекторий.

Так же, как в разд. 2, из критерия Вышнеградского следует, что собственные числа матрицы (6) либо отрицательны, либо имеют отрицательные вещественные части, значит, и траектории модели второй генной сети, рассматриваемой в изоляции от других генных сетей организма, притягиваются к устойчивой стационарной точке системы (5).

Авторы выражают благодарность О. А. Подколодной, Н. А. Колчанову и В. А. Чуркину за полезные советы и обсуждения.

ЛИТЕРАТУРА

1. Banks H. T., Mahaffy J. M. Stability of cyclic gene models for systems involving repression // J. Theor. Biology. 1978. V. 74. P. 323-334.

2. Bass J. Circadian topology of metabolism // Nature. 2012. V. 491, N 7424. P. 348-356.

3. Подколодная О. А. Молекулярно-генетические аспекты взаимодействия циркадных часов и метаболизма энергетических субстратов млекопитающих // Генетика. 2014. Т. 50, № 2. С. 1-13.

4. Лихошвай В. А., Голубятников В. П., Демиденко Г. В., Евдокимов А. А., Фадеев С. И. Теория генных сетей // Системная компьютерная биология. Новосибирск: Изд-во СО РАН, 2008. C. 395-480.

5. Bukharina T. A., Golubyatnikov V. P., Kazantsev M. V., Kirillova N. E., Furman D. P. Mathematical and numerical models of two asymmetric gene networks // Sib. Electron. Math. Rep. 2018. V. 15. P. 1271-1283.

6. Golubyatnikov V. P., Mjolsness E., Gaidov Yu. A. Topological index of a model of p53 — Mdm2 circuit // Информ. Вестн. Вавилов. о-ва генетиков и селекционеров. 2009. Т. 13, № 1. С. 160-162.

7. Чумаков Г. А., Чумакова Н. А. Гомоклинические циклы в одной модели генной сети // Мат. заметки СВФУ. 2014. Т. 21, № 4. C. 97-106.

8. Аюпова Н. Б., Голубятников В. П. Трехклеточная модель ранней стадии развития одного пронейрального кластера // Сиб. журн. индустр. математики. 2017. Т. 20, № 2. С. 15-20.

9. Голубятников В. П., Градов В. С. О неединственности циклов в некоторых кусочно-линейных моделях генных сетей // Мат. тр. 2020. V. 23, № 1. С. 107-122.

10. Вышнеградский И. А. О регуляторах прямого действия // Изв. Технолог. ин-та. СПб: Императ. Акад. наук., 1877. C. 21-62.

11. Постников М. М. Устойчивые многочлены. М.: УРСС, 2004.

12. Smale S. A mathematical model of two cells via Turing's equation // Lecture Notes in Applied Mathematics. Providence, RI: Amer. Math. Soc., 1974. V. 6. P. 15-26.

13. Акиньшин А. А., Бухарина Т. А., Голубятников В. П., Фурман Д. П. Математическое моделирование взаимодействия двух клеток в пронейральном кластере крылового има-гинального диска D. melanogaster // Сиб. журн. чистой и прикл. математики. 2014. Т. 14, № 4. С. 3-10.

14. Аюпова Н. Б., Голубятников В. П. Строение фазового портрета одной кусочно-линейной

динамической системы // Сиб. журн. индустр. математики. 2019. Т. 22, № 4. С. 19—25.

Поступила в редакцию 1 февраля 2021 г. После доработки 1 февраля 2021 г. Принята к публикации 26 февраля 2021 г.

Голубятников Владимир Петрович Институт математики им. С. Л. Соболева СО РАН, пр. Академика Коптюга, 4, Новосибирск 630090; Новосибирский государственный университет, ул. Пирогова, 1, Новосибирск 630090 Vladimir.golubyatnikov1@fulbrightmail.org

Кириллова Наталья Евгеньевна

Институт математики им. С. Л. Соболева СО РАН, пр. Академика Коптюга, 4, Новосибирск 630090 kne@math.nsc.ru

Математические заметки СВФУ Январь—март, 2021. Том 28, № 1

UDC 514.745.82

PHASE PORTRAITS OF TWO GENE NETWORKS MODELS V. P. Golubyatnikov and N. E. Kirillova

Abstract: We construct mathematical models of functioning of two few-components gene networks which regulate circadian rhythmes in organisms by means of combinations of positive and negative feedbacks between components of these networks. Both models are represented in the form of non-linear dynamical systems of biochemical kinetics. It is shown that the phase portraits of both models contain exactly one equilibrium point each and in both cases for all values of parameters of these dynamical systems, eigenvalues of their linearization matrices at their equilibrium points are either negative or have negative real parts. Thus, these equilibrium points are stable. We construct their invariant neighborhoods and describe the behavior of trajectories of these systems. Biological interpretations of these results are given as well.

DOI: 10.25587/SVFU.2021.68.70.001 Keywords: gene networks models, non-linear dynamical systems, phase portrait, equilibrium point, stability, Vyshnegradskii criterion.

REFERENCES

1. Banks H. T. and Mahaffy J. M., "Stability of cyclic gene models for systems involving repression," J. Theor. Biology, 74, 323-334 (1978).

2. Bass J., "Circadian topology of metabolism," Nature, 491, No. 7424, 348-356 (2012).

3. Podkolodnaya O. A., "Molecular genetic aspects of the interaction of the circadian clocks and metabolism of energy substrates of mammals [in Russian]," Genetics, 50, No. 2, 1-13 (2014).

4. Likhoshvai V., Golubyatnikov V., Demidenko G., Evdokimov A., and Fadeev S., "Gene networks theory," in: Computational Systems Biology [in Russian], pp. 395-480, Izdat. SO RAN, Novosibirsk (2008).

5. Bukharina T. A., Golubyatnikov V. P., Kazantsev M. V., Kirillova N. E., and Furman D. P., "Mathematical and numerical models of two asymmetric gene networks," Sib. Electron. Math. Rep., 15, 1271-1283 (2018).

6. Golubyatnikov V. P., Mjolsness E., and Gaidov Yu. A., "Topological index of a model of p53 — Mdm2 circuit," Inform. Vestn. Vavilov. Obshch. Genetikov i Seleksts., 13, No. 1, 160162 (2009).

7. Chumakov G. A. and Chumakova N. A., "Homoclinic cycles in one gene network model [in Russian]," Mat. Zamet. SVFU, 21, No. 4, 97-106 (2014).

8. Ayupova N. B. and Golubyatnikov V. P., "A three-cells model of the initial stage of development of one proneural cluster," J. Appl. Ind. Math., 11, No. 2, 1-7 (2017).

9. Golubyatnikov V. P. and Gradov V. S., "Non-uniqueness of cycles in piecewise-linear models of circular gene networks," Sib. Adv. Math., 31, No. 1, 1-12 (2021).

10. Vyshnegradskiy I. A. " On regulators of direct action [in Russian]," Izv. Tekhnolog. Inst., Imper. Akad. Nauk, Saint-Petersburg, 21-62 (1877).

11. Postnikov M. M., Stable Polynomials [in Russian], URSS, Moscow (2004).

12. Smale S., "A mathematical model of two cells via Turing's equation," in: Lecture Notes in Applied Mathematics, vol. 16, pp. 15-26, Amer. Math. Soc., Providence, RI (1974).

© 2021 V. P. Golubyatnikov and N. E. Kirillova

13. Akinshin A. A., Bukharina T. A., Golubyatnikov V. P., and Furman D. P., "Mathematical modeling of the interaction of two cells in the proneural cluster of the wing imaginal disc D. melanogaster [in Russian]," Sib. Zhurn. Chistoy Prikl. Mat., 14, No. 4, 3-10 (2014).

14. Ayupova N. B. and Golubyatnikov V. P., "The structure of the phase portrait of one piecewise linear dynamic system," J. Appl. Ind. Math., 13, No. 4, 1-8 (2019).

Submitted February 1, 2021 Revised February 1, 2012 Accepted 26 February, 2021

Vladimir P. Golubyatnikov Sobolev Institute of Mathematics, 4 Koptyug Avenue, Novosibirsk 630090, Russia; Novosibirsk State University, 1 Pirogov Street, Novosibirsk 630090, Russia vladimir.golubyatnikov1@fulbrightmail.org

Nataliya E. Kirillova Sobolev Institute of Mathematics, 4 Koptyug Avenue, Novosibirsk 630090, Russia kne@math.nsc.ru