CC BY
10
структуры и моделирование 2017. №2(42). С. 59-67
УДК 514.745.82
ВОПРОСЫ СУЩЕСТВОВАНИЯ УСТОЙЧИВОГО ЦИКЛА В ОДНОЙ МОДЕЛИ МОЛЕКУЛЯРНОГО РЕПРЕССИЛЯТОРА
В.П. Голубятников
профессор, д.ф.-м.н., г.н.с., e-mail: glbtn@math.nsc.ru
Институт математики им. С.Л. Соболева СО РАН Новосибирский государственный университет
Аннотация. Рассматривается нелинейная шестимерная динамическая система, моделирующая функционирование простейшего молекулярного ре-прессилятора. Установлены достаточные условия существования устойчивого цикла в её фазовом портрете.
Ключевые слова: нелинейная динамическая система, модели генных сетей, гиперболические стационарные точки, циклы, теорема Брауэра о неподвижной точке, устойчивость.
1. При моделировании широкого класса генных сетей используются нелинейные динамические системы вида:
dm\ Ф1
— = -kimi + Fi(pa); -dt
^imi - eipi;
dm2 = -k2m2 + F2(pi); = ^2m2 - A2P2;
dt dt
dm3 = -k3m3 + F3(P2); dpT3 = ^ama - взРз.
dt dt
(1)
Здесь нелинейные слагаемые в уравнениях ^(р) — гладкие положительные монотонно убывающие функции неотрицательного аргумента, описывающие отрицательные обратные связи; к], ^, вз — положительные параметры, характеризующие динамику биохимических процессов в моделируемой генной сети. Всюду в дальнейшем будем предполагать, что 3 = 1,2, 3 и что если 3 = 1, то 3 -1 = 3.
В наших предыдущих работах [1,3] эта система изучалась в случае ^ = вз, для которого были получены условия существования по крайней мере одного цикла в её фазовом портрете.
В так называемом симметричном безразмерном случае, в котором к! = к2 = кз = 1, = ^2 = ^э, Д? = в] и ^\(р) = ^2(р) = ВД = ^(р) = а(1 + р7)-! + ао , система (1) была предложена в [13] для описания динамики синтетической кольцевой генной сети, связывающей три белка с концентрациями р!(^), р2(£),
рз(Ъ), репрессирующими друг друга по циклу, и соответствующие мРНК с концентрациями т1 (Ъ), т2(Ъ), m3(í). Здесь и далее все параметры предполагаются положительными. Такая динамическая система инвариантна относительно циклических замен пар переменных (т1,р1) ^ (т2,р2) ^ (тз,рз) ^ (т1,р1). В недавних работах [7,8,14] исследования этого симметричного случая были продолжены с целью доказательства существования и устойчивости бегущих волн у таких динамических систем. В работе [9] подобные симметричные системы изучались с целью рассчёта кольцевых электрических сетей, состоящих из идентичных однонаправленно связанных генераторов.
В несимметричном случае подобные системы изучались и в [16], где в качестве нелинейных слагаемых рассматривались также и монотонно возрастающие функции. Следуя этой работе, совершим линейную замену переменных ^зРз := Рз, !з(Рз-1) := (рз-1), после чего система (1) принимает более симметричный вид:
¿О^ = -к1т1 + /1(рз); ^ = т1 - Р1;
¿т2 , , ¿р2
-к2т2 + /2(Р1); ~тТ = т2 - в2Р2; (2)
аъ
-кзтз + /з(Р2); аР3 = тз - взРз.
аъ 2 2 2 1 аъ атз ¿'Рз
-аъ- = -кзтз +/з(Р2); -ъ
Основной целью настоящей работы является установление достаточных условий существования по крайней мере одного устойчивого цикла системы (2) в её инвариантной области, которая будет описана ниже.
Пусть Аз := /3(0)/кз, Б3 := А,/вз и д := [0, А1] х [0, Б1] х [0, А2] х [0,Б2] х [0, Аз] х [0,Бз].
Также, как и в [4, 16], где изучались фазовые портреты подобных динамических систем, в том числе и других размерностей, можно проверить, что параллелепипед д, лежащий в положительном октанте является инвариантной областью системы (2); все траектории, начинающиеся в д при Ъ = 0, остаются в нём при всех Ъ > 0.
Координаты любой стационарной точки системы (2) удовлетворяют соотношениям тз = взРз и находятся из уравнения
к1т1 = /1 ( /з( ** ) | , (3)
в котором правая часть — композиция трёх монотонно убывающих гладких функций — монотонно убывает с ростом т1, а левая возрастает. Следовательно, уравнение (3) имеет единственное решение, и потому стационарная точка Б0 системы (2) существует и единственна. Пусть £0 = (т°;р1; т°;р2; тз;р°°) — её координаты; через (-) будем обозначать производную О/(р)/ф, вычисленную при р = р0-1, здесь q3• > 0, так как производная монотонно убывающей функции отрицательна.
Для описания фазового портрета системы (2), следуя [6,10,18], где подобные построения проводились для аналогичных динамических систем, разобьём
область Q плоскостями, шу
т0, ру
проходящими через точку 50. Обо-
значим полученные параллелепипеды (блоки) разбиения бинарными индексами:
Е = {£1^2£з£4£5£6} =
{х е Q 1 ш1 ш1; р1 р?; ш2 ш°; Р2 р0; шз ^ ш0; Рз Р0, },
(4)
где X = (ш?,р?,ш2,р2,ш3,р3), £2, £3, £4, £5, £б е {0,1}, и отношения порядка задаются следующим образом: символ означает а символ означает
Так же, как в [3,6,18], проверяется, что при такой дискретизации фазового портрета динамической системы (2) для любых двух соседних блоков Е? и Е2 разбиения (4), имеющих общую пятимерную грань Е? П Е2, траектории всех точек этой грани переходят только в один из этих двух блоков — либо из Е? в Е2: Е? ^ Е2, либо наоборот: Е2 ^ Е?.
Матрица линеаризации системы (2) в точке 50 имеет вид:
Мо
-к? 0 0 0 0
1 -в? 0 0 0 0
0 -Я2 -к2 0 0 0
0 0 1 -в2 0 0
0 0 0 -дз -кз 0
0 0 0 0 1 -вз
В аналитическом выражении определителя этой матрицы всего два ненулевых слагаемых — произведение диагональных элементов и произведение а6 := д^дз недиагональных элементов, поэтому характеристический многочлен матрицы М0 может быть записан в форме
Р (Л) = (к? + Л)(^2 + Л)(кз + Л)(в? + Л)(в2 + Л)(вз + А) + а6.
Напомним, что стационарная точка динамической системы называется гиперболической, если собственные числа соответствующей матрицы линеаризации имеют как положительные, так и отрицательные вещественные части, но не имеют мнимых вещественных частей. Упорядочим пары собственных чисел матрицы М0 в порядке убывания их вещественных частей: КеЛ?,2 > И.еЛ2,з > И,еЛ5,6. В дальнейшем мы будем рассматривать именно такую комбинацию знаков этих вещественных частей, при которой стационарная точка 50 является гиперболической.
Следующая диаграмма показывает, по каким двенадцати блокам разбиения (4) может проходить цикл динамической системы (2).
{110011} -► {010011} -► {000011} -► {001011}
{110010}
{001111}
(5)
{110000}
{001101}
{110100} ц- {111100} ц- {101100} ц- {001100}
Для каждого блока, указанного в диаграмме, траектории системы (2) могут переходить из него в другие блоки разбиения (4) только в направлении, указанном в диаграмме. Обозначим через Ш объединение всех этих двенадцати блоков.
Пусть и ~ Б2 х Б4 — достаточно малая открытая окрестность стационарной точки 50; здесь двумерный диск Б2 параллелен плоскости и построен по собственным числам Л1,Л2, а четырёхмерный диск Б4 аналогичным образом построен по остальным собственным числам матрицы М0, имеющим отрицательные вещественные части.
Следующие два утверждения доказываются дословно так же, как их аналоги, установленные в [3] для системы (1) в случае вз = ^з.
Лемма. При достаточно больших значениях параметра а характеристический многочлен Р(Л) имеет в точности два комплексных корня с положительными вещественными частями и четыре комплексных корня с отрицательными вещественными частями.
Теорема 1. Если характеристический многочлен матрицы М0 имеет два корня с положительными вещественными частями и четыре корня с отрицательными вещественными частями, то динамическая система (2) имеет по крайней мере один цикл С, содержащийся в области Ш' = Ш \ (Ш П и) и проходящий по ней в соответствии с диаграммой (5).
Следует подчеркнуть, что гиперболичность стационарной точки 50 используется здесь существенным образом для линеаризации системы (2) в окрестности точки 50 с целью построения окрестности и. Такая линеаризация существует на основании теоремы Гробмана-Хартмана [12].
Если же у матрицы линеаризации динамической системы имеются и мнимые собственные числа, то для построения такой окрестности и потребуются старшие члены разложения правых частей уравнений системы в окрестности её стационарной точки. В случаях вырожденности этих старших членов описание фазового портрета динамической системы становится необозримо сложным. В ряде ранних публикаций, посвящённых математическому моделированию подобных генных сетей, см., например, [11], это обстоятельство не было учтено.
2. В дальнейших рассуждениях об условиях существования устойчивого цикла системы (2) мы ограничимся описанием «частично симметричного» слу-
чая к1 = к2 = кз = к. Все остальные положительные параметры вз, а также монотонно убывающие гладкие функции /, как и в [1,3], будут полагаться произвольными.
Следуя Р. Смиту [17] (см. также [5,10]), представим систему (2) в бескоординатной форме:
Х = А ■ X + Ф(Х), (6)
где все нелинейные слагаемые правых частей системы (2) (но не только они) содержатся во втором слагаемом Ф(Х), вектор-функция X(Ъ) имеет координаты (т1,р1,т2,р2,тз,рз), и матрица А имеет постоянные коэффициенты. Именно, пусть
✓
/1(Рз) + прз;
ф
mi(1 - n) + Pi(k - ^1); /2(^1) + npi; m2(1 - n) + P2(k - в2);
/3(P2) + ПР2;
Шэ(1 - n) + Рз(к - вз); здесь n — некоторый положительный параметр, и пусть A = -kE + nD, где
D
0 0 0 0 1 0 0 0 0 0 -10 0 0 00100 0 0 0 -1 0 0 0 0 0 1
0 -1 0 0 0 0 0
Обозначим через x(z) := (zE - A)-1 переходную матрицу для матрицы A, см. [15,17]; здесь z е C. Пусть 0(x) := sup ||x(iш - x)|| для -то < ш < то, x > 0, и пусть Ф' — матрица Якоби отображения ф : R+ ^ R6. В работе [17] было установлено (Теорема 3 и Лемма 6), что: Если у системы (6) матрица A имеет в точности два собственных числа Л1,2, у которых вещественные части больше, чем -р, где р > 0, и в инвариантной области W' выполняется неравенство
1|Ф'11 <%)
-1
(7)
то система (5) имеет по крайней мере один устойчивый цикл в области
3. Доказательство существования устойчивого цикла у динамической системы (2) в области Ш' будет сведено к проверке условия (7). Можно показать, что симметричная матрица Ф' ■ (Ф')х, определяющая квадрат нормы матрицы
Ф', разбивается на три блока:
В
(Л + П)2 (/ + п)(к - вз)
(/' + п)(к - вз) (к - вз)2 + (1 - п)2
В
Вз
(/2 + п)2
(/2 + п)(к - вх)
(/2 + п)(к - вх) (к - вх)2 + (1 - п)2
(/3 + п)2
(/3 + п)(к - в2)
(/3 + п)(к - в2) (к - в2)2 + (1 - п)2
соответствующие парам координатных осей (Отх,Орз), (От2,Орх) и (Отз, Ор2).
Итак, мы получаем оценку
||Ф'||2 = шах{||Вз ||}. з
Корни характеристического многочлена (к + А(А))6 + п6 = 0 матрицы А, у которых вещественные части максимальны, имеют вид Ах,2 = -к+пл/3/2± пг/2, у остальных его корней вещественные части меньше или равны -к, и поэтому в качестве параметра р можно выбрать любое число 7к, где 0 < 7 < 1. Матрица Б в ортогональном базисе ех = (1, 0,1, 0,1, 0), е2 = (0,1, 0,1, 0,1),
ез = (2, /3,-1, 0,-1, - л/3), в4 = (0,-1, л/3, 2, -v/3,-1), е5 = (2, -/3, -1, 0, -1, /3), еб = (0, -1, -/3, 2, /3, -1) распадается на три двумерных блока:
Бех = 2Бе5
е2,
Бе2 = -ех;
2Без
л/3ез -
е4,
2Бе4 = ез + л/3е4;
-л/3е5 - еб, 2Беб = е5 - /3еб.
Следовательно, переходная матрица (гш - р)Е - А = (гш + к - р)Е - пБ распадается на блоки:
я.
з,4 =
я5
5,6 =
Ях ■
гш
гш
гш + к - р п -п гш + к - р
+ к - р - п/3/2 -п/2
+ к - р + п "/3/2 -п/2
гш
гш
п/2
+ к - р - п/3/2 п/2
+ к - р + п/3/2
Чтобы найти 0(р)-х, надо вычислить т£ш минимального квадрата нормы каждого из этих трёх блоков и выбрать из них наименьший. Несложные вычисления показывают, что
||Ях,2||2 = (к - р)2 + шах(ш ± п)2,
я
3,41
[(к - р) - п/3/2]2 + шах(ш ± п/2)2,
2
||#5;6||2 = [(k - Р) + ^v/3/2]2 + max(w ± n/2)2.
Таким образом, нами установлен основной результат настоящей работы.
Теорема 2. Если ki = k2 = k3 ив области W' при некотором п > 0 выполняется неравенство
min ||H2j-i,2j||2 > max ||£в||, j, s = 1, 2, 3, (8)
j s
то система (2) имеет по крайней мере один устойчивый цикл в области W'.
Отметим, что, как было сказано в [17], оценка (7) не является точной, и, следовательно, неравенство (8) является лишь достаточным условием существования устойчивого цикла у системы (2).
В случае ki = k2 = k3 = ki, для которого аналог теоремы 1 был установлен в [1,3], матрица (гш — р)Е — A на двумерные блоки не распадается, и тогда проверка условия (7) становится громоздкой.
Вопросы единственности циклов шестимерных динамических систем (1) и (2) остаются открытыми. У аналогичных динамических систем больших размерностей количество циклов может оказаться довольно большим, см. [2].
Результаты некоторых вычислительных экспериментов с траекториями динамических систем вида (1), (2) приведены в [19-21].
Благодарности
Автор выражает искреннюю благодарность Н.Б. Аюповой за полезные обсуждения, а также А.А. Акиньшину и М.В. Казанцеву за проведение численных экспериментов. Работа поддержана РФФИ, грант 15-01-00745.
Литература
1. Акиньшин А.А., Аюпова Н.Б., Голубятников В.П., Казанцев М.В. Геометрия фазового портрета одной генной сети // Труды 12-й Международной Азиатской школы-семинара «Проблемы оптимизации сложных систем». Новосибирск, ИВМ и МГ СО РАН, 12 - 16 декабря 2016, С. 18-24.
2. Акиньшин А.А., Голубятников В.П. Циклы в симметричных динамических системах // Вестник НГУ. 2012. Т. 12, № 2, С. 3-12.
3. Аюпова Н.Б., Голубятников В.П., Казанцев М.В. О существовании цикла в одной несимметричной модели молекулярного репрессилятора // Сибирских журнал вычислительной математики. 2017. Т. 20, № 2. C. 121-130.
4. Аюпова Н.Б., Голубятников В.П. О единственности цикла в трёхмерной модели молекулярного репрессилятора // Сибирский журнал индустриальной математики. 2014. Т. 17, № 1. С. 3-7.
5. Гайдов Ю.А. Об устойчивости периодических траекторий в некоторых моделях генных сетей // Сибирский журнал индустриальной математики. 2008. Т. 11, № 1. С. 57-62.
6. Glass L., Pasternack J.S. Stable oscillations in mathematical models of biological control systems // Journal of Mathematical Biology. 1978. V. 6. P. 207-223.
7. Глызин С.Д., Колесов А.Ю., Розов Н.Х. Явление буферности в кольцевых генных сетях // Теоретическая и математическая физика. 2016. Т. 187, № 3. С. 560-579.
8. Глызин С.Д., Колесов А.Ю., Розов Н.Х. Существование и устойчивость релаксационного цикла в математической модели репрессилятора // Математические заметки. 2017. Т. 101, № 1. С. 58-67.
9. Глызин С.Д., Колесов А.Ю., Розов Н.Х. Явление буферности в кольцевых цепочках однонаправленно связанных генераторов // Известия РАН, серия математическая. 2014. Т. 78, № 4. С. 73-108.
10. Голубятников В.П., Голубятников И.В., Лихошвай В.А. О существовании и устойчивости циклов в пятимерных моделях генных сетей // Сибирский журнал вычислительной математики. 2010. Т. 13, № 4. С. 403-411.
11. Golubyatnikov V., Likhoshvai V, Ratushnyi A. Existence of closed trajectories in 3D gene networks // The Journal of Three-Dimensional Images. 2004. V. 18, No. 4. P. 96-101.
12. Гробман Д.М. Гомеоморфизм систем дифференциальных уравнений // ДАН СССР. 1959. Т. 128, № 5. C. 880-881.
13. Elowitz M.B., Leibler S. A synthetic oscillatory network of transcriptional regulators // Nature. 2000. V. 403. P. 335-338.
14. Колесов А.Ю., Розов Н.Х., Садовничий В.А. Периодические решения типа бегущих волн в кольцевых генных сетях // Известия РАН, серия математическая. 2016. Т. 80, № 3. С. 67-94.
15. Леонов Г.А. Об устойчивости фазовых систем // Сибирский математический журнал. 1974. Т. 15, № 1. С. 49-60.
16. El Samad H., Del Vecchio D., Khammash M. Repressilators and Promotilators: Loop dynamics in gene regulatory networks // Proc. American Control Conference, 2005. P. 4405-410.
17. Smith R.A. Orbital stability for ordinary differential equations // Journal of differential equations. 1987. V. 69. P. 265-287.
18. Hastings S., Tyson J.J., Webster D. Existence of periodic solutions for negative feedbacks cellular control systems // Journ. Diff. Equations. 1977. V. 25. P. 39-64.
19. Elowitz-Leibler model [Электронный ресурс]. URL: https://maxim-kazantsev. shinyapps.io/ElowitzLeibler/ (дата обращения: 09.03.2017).
20. AndreyAkinshin/model.R [Электронный ресурс] // GitHub. URL: https://gist. github.com/AndreyAkinshin/9cdc1c2dbe71d154ab57 48e8c57e6425 (дата обращения: 09.03.2017).
21. AndreyAkinshin/Elowitz.R // GitHub. URL: https://gist.github.com/ AndreyAkinshin/37f3e68a157 6f9ea1e5c01f2fd64fe5e (дата обращения: 09.03.2017).
THE EXISTENCE OF A STABLE CYCLE IN A MODEL OF A MOLECULAR
REPRESSOR
V.P. Golubyatnikov
Dr.Sc. (Phys.-Math.), Professor, Chief Researcher, e-mail: glbtn@math.nsc.ru
Sobolev Institute of Mathematics, Siberian Branch of the Russian Academy of Sciences
Novosibirsk State University
Abstract. We consider nonlinear 6-dimensional dynamical system which describes a model of functioning of one simple molecular repressilator. We find sufficient conditions of existence of a stable cycle in the phase portrait of this system.
Keywords: nonlinear dynamical system, gene network models, hyperbolic stationary points, cycles, Brouwer's fixed-point theorem, stability..
Дата поступления в редакцию: 09.03.2017
