Гомоклинические циклы в одной модели генной сети Текст научной статьи по специальности «Физика»

Математические заметки СВФУ Октябрь—декабрь, 2014. Том 21, №4

УДК 517.928.4:517.929.5

ГОМОКЛИНИЧЕСКИЕ ЦИКЛЫ В ОДНОЙ МОДЕЛИ ГЕННОЙ СЕТИ Г. А. Чумаков, Н. А. Чумакова

Аннотация. Гомоклинические циклы, образованные гетероклиническими траекториями и стационарными точками, существуют в основном в динамических системах с симметрией и могут подвергаться различным бифуркациям. В статье изучаются глобальные аспекты модели генных сетей с тремя переменными и численно исследуется бифуркация гомоклинического цикла, состоящего из трех седел и трех гетероклинических орбит. Представляет большой интерес тот факт, что бифуркация гомоклинического цикла характеризуется высокой параметрической чувствительностью к начальным данным. В этом случае при небольших возмущениях концентраций мы сталкиваемся со случайным переходом от релаксационных колебаний, соответствующих наблюдаемому устойчивому предельному циклу, близкому к треугольному гомоклиническому циклу, к различным устойчивым стационарным состояниям. Сценарий разрушения треугольного гомоклинического цикла имеет простой биологический смысл: в окрестности бифуркационных значений параметров в то время, когда происходит генерация релаксационных импульсов в генной сети, малые флуктуационные выбросы концентраций могут «сбивать» систему с колебательного режима и случайным образом стабилизировать ее к существенно различным стокам.

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

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

Пожалуй, еще более важную роль, чем гармонические колебания, играют релаксационные колебания. Предельный цикл в этом случае может иметь форму треугольника. Именно такие автоколебания возникают в простейшей модели первичного жизненного цикла, в которой происходит генерация релаксационных импульсов и чередование быстрых и медленных процессов во времени при постоянных внешних условиях [6].

Работа выполнена при финансовой поддержке Сибирского Отделения Российской академии наук (междисциплинарный проект №80).

© 2014 Чумаков Г. А., Чумакова Н. А.

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

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

Рассмотрим математическую модель, описывающую динамику генной сети с регуляцией стадии деградации [3]:

х = а - (1 + г7 + ум)ж,

у = а - (1+ ж7 + гм)у, (1)

г = а - (1 + у1 + жм)г,

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

Qa = {(ж, у, г) е М3 : ж, у, г е [0, а]},

параметры а, 7 и у — положительные константы. В [8] показано, что куб Qa является инвариантным множеством относительно системы (1).

Будем изучать изменение фазового портрета системы (1) в зависимости от параметра у. Численные результаты приведены для а = 5.908 и 7 = 2.981, если иное не оговорено специально.

При у =1.9 наблюдаем затухающие колебания (рис. 1) в окрестности устойчивого фокуса О, расположенного на диагонали А = {(ж, у, г) е Qa : ж = у = г} (координаты фокуса жд = уд = гд = 1.27352).

ЛЛЛ/\

0.8-1-1-1-1-1-1-

0 10 20 30 40 50 60 70

t

Рис. 1. Затухающие колебания в окрестности устойчивого фокуса О при ^ = 1.9.

Кроме того, в фазовом пространстве имеются шесть неподвижных точек: устойчивые точки А1, А2, А3 типа узла и седловые особые точки Р^ Р2 и Р3, каждая из которых имеет двумерное устойчивое многообразие и одномерное неустойчивое многообразие (рис. 2).

6- ' 5432-1 -

Д1....Г"""

ур, ...........

р":■ ■

-Т. 2 -.......Ч; 7

Рис. 2. Глобальный фазовый портрет системы (1) с устойчивыми узлами Ах, А2, А3, седловыми стационарными состояниями Рх, Р2, Р3 и устойчивым фокусом О при ^ = 1.9.

Р

Р

А

3

0

2

4

2

4

6

6

У

X

Чтобы не усложнять рисунок, в окрестности точек Р^, г = 1, 2, 3, мы изобразили только расположение сепаратрис на двумерном ведущем многообразии. Расчеты показывают, что в окрестности седловой точки, например, Р2 = (1.1617, 2.0532, 0.5436), матрица линеаризованной системы уравнений имеет собственные значения {-14.140, -7.015, +2.324}. В этом случае устойчивое одномерное подпространство, соответствующее наименьшему собственному значению, является быстрым подпространством. Отметим, что соответствующие значения для точек Р1 и Р3 получаются в результате циклической подстановки координат в силу симметричности системы (1). Поэтому изображающая точка в фазовом пространстве быстро приходит в окрестность показанного на рис. 2 локального двумерного многообразия в окрестности каждой из точек Р^.

Далее, точка А3 = (5.8529, 0.0303, 0.1988) — устойчивый узел, для которого матрица линеаризованной в его окрестности системы уравнений имеет собственные значения {-194.93, -29.75, -0.956}. Здесь ведущим является многообразие, касающееся собственного подпространства, отвечающего наименьшему по модулю собственному значению. Поэтому изображающая точка в фазовом пространстве с увеличением времени быстро приходит в окрестность двумерного многообразия и следует вдоль показанной на рис. 2 гетероклинической траектории (сеператрисы соответствующей седловой точки).

2. Бифуркация Андронова — Хопфа

В [3] обнаружена бифуркация Андронова — Хопфа в системе (1) и численно показано существование периодического решения при у = у*.

Используя теорию бифуркации Андронова — Хопфа периодического решения [9] и метод продолжения семейства периодических решений по параметру [10], мы уточнили бифуркационное значение параметра у* и построили максимальное семейство периодических решений (рис. 3).

Семейство периодических решений Г (у) при у* < у < у* зарождается из стационарной точки при у = у* и «влипает» при у = у* в гомоклинический

Р2

1 4 1 6 1 8 у

Рис. 3. Максимальное семейство Гпериодических решений, зависящих от параметра для ^ф < ^ < и гомоклинический цикл при ^ = , где = 1.9998 и = 2.27659792 — бифуркационные значения.

цикл

Г * = Р1 и 712 и Р2 и 723 и Рз и 731.

(2)

Эти периодические траектории получены с помощью метода расчета структурно устойчивых (грубых) периодических траекторий, описанного в [11,12].

Рис. 4. Почти гармонические колебания при ^ = 2.01 > с периодом То = 2.586176.

Рис. 5. Релаксационные колебания при ^ = 2.2765979 < с периодом То = 21.417515.

При увеличении параметра у автоколебания изменяются от гармонических в окрестности у = у* до релаксационных при приближении параметра к критическому значению у = у* (рис. 4, 5). Видно, что на циклах большой амплитуды система достаточно быстро переходит из окрестности одной стационарной точки р в окрестность следующей седловой точки, где задерживается на некоторое время. Поэтому на рис. 5 выделяются интервалы быстрого изменения у-координаты решения и плато, на которых скорость изменения решения мала и, соответственно, у(4) изменяется медленно в окрестности седловой стационарной точки. Чем ближе значение параметра у к бифуркационному, тем шире становятся плато на соответствующей интегральной кривой и тем больше период автоколебаний.

3. Глобальная бифуркация при / = /*.

Треугольный гомоклинический цикл Г*

Рассмотрим три седловые точки Д, г = 1, 2, 3, в М3, для которых обозначим двумерное устойчивое многообразие через Шв(Р^) и одномерное неустойчивое многообразие - через ШИ(Р^) (см. рис. 2).

При критическом значении у = у* точки Рх, Р2 и Р3 связаны между собой так, что образуется замкнутый контур, составленный из седел и их сепаратрис. Соединительные кривые ^ = ШИ(Р^) П Шв(Рц) служат примерами седловых соединений или гетероклинических орбит, т. е. фазовая траектория при 4 ^ —то асимптотически приближается к седлу Р^, а при 4 ^ +то — к точке Рц. Эти орбиты возникают в фазовом пространстве системы (1) в результате наличия симметрии относительно циклической подстановки х ^ у ^ г ^ х и являются следствием глобальной бифуркации трехмерного устойчивого предельного цикла Г (у), период которого неограниченно возрастает при у ^ у* (рис. 6). Топологический предел Г (у) при у ^ у* обозначен в (2) через Г * и является гомоклиническим циклом, который состоит из трех гетероклинических орбит и трех седловых состояний равновесия.

Рис. 6. Изменение периода автоколебаний в зависимости от параметра

На рис. 7 показан предельный цикл с мультипликаторами Лх = 1, Л2 = 0.0254, Аз = 7.02 х 10-8 при значении параметра у « у*. Расчеты показывают, что устойчивое многообразие, соответствующее мультипликатору Л2 этого цик-

Рис. 7. Линейные приближения ориентируемого устойчивого многообразия треугольного периодического решения, соответствующего мультипликатору А2.

ла почти треугольной формы, является ориентируемым. Это обстоятельство упрощает численное построение глобального фазового портрета (рис. 8).

На рис. 8, 9 показана динамическая неопределенность в рассматриваемой модели, которая возникает при исчезновении области устойчивости предельного цикла на верхней границе максимального семейства периодических решений.

Рис. 8. Глобальный фазовый портрет с треугольным гомоклиническим циклом Г * при ^ = .

С увеличением параметра у при переходе через бифуркационное значение у* все три соединения ^уц разрушаются (см. рис. 9). Тогда неустойчивые многообразия ШИ(Р^) будут иметь компоненты, приближающиеся при 4 ^ к одному из устойчивых узлов вблизи угловых точек (а, 0, 0), (0, а, 0) или (0, 0, а) куба Сепаратрисы седел Р^, приходящие из неустойчивого фокуса, лежат на поверхностях, которые разделяют области притяжения устойчивых угловых стационарных точек А.

4. Параметрическая чувствительность решений к начальным данным

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

Рис. 9. Разрушение гетероклинических орбит, идущих из седла в седло при ^ = 2.28 >

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

Поскольку такие возмущения неизбежны, возникает вопрос о том, какое влияние могут оказать эти случайные отклонения на характер поведения изучаемых нами динамических систем и параметрическую чувствительность решений к начальным данным и параметрам модели. Возвращаясь к системе (1), обратим внимание на особенности поведения траекторий с близкими начальными данными, когда однопараметрическое семейство динамических систем, зависящее от параметра у, пересекает границу раздела разных областей невырожденных систем при у = у* или у = у*.

При возрастании параметра у в окрестности у* состояние равновесия из устойчивого становится неустойчивым, однако изображающая точка остается в малой е-окрестности состояния равновесия. При обратном изменении параметра, когда состояние равновесия опять становится устойчивым, изображающая точка снова возвращается в состояние равновесия. Система ведет себя обратимо, т. е. граница области устойчивости в пространстве параметров является «безопасной».

При переходе через бифуркационное значение у* система (1) имеет негрубый особый элемент — треугольный гомоклинический цикл Г*, состоящий из седел и гетероклинических орбит, идущих из одного седла в другое. Для у = у* система (1) не является структурно устойчивой. Малые возмущения будут переводить ее в области невырожденных систем с различными фазовыми портретами. При у = у* гомоклинический цикл Г* устойчив «изнутри». В частности, все решения из малой окрестности неустойчивого фокуса, расположенного на диагонали А куба притягиваются к Г*. При уменьшении параметра у < у* из Г * рождается устойчивый предельный цикл. Соответствующие перестройки фазовых портретов показаны на рис. 2, 3 и 8, 9.

Отметим, что описанная граница раздела областей невырожденных систем

является «опасной», т. е. это такая граница, сколь угодно малое нарушение которой влечет за собой переход системы в новое состояние, из которого нельзя вернуть систему в прежнее состояние при обратном изменении того же управляющего параметра у. В связи с этим в [13] выдвинута следующая аксиома неопределенности: при переходе через опасную границу

1) изображающая точка покидает малую окрестность и (Г*) гомоклиниче-ского цикла Г* ,

2) выход ее может происходить по любой траектории, покидающей и (Г *).

Далее возникает вопрос: куда при переходе через опасную границу «перескакивает» изображающая точка? Особенность системы (1) состоит в том, что после прохождения опасной границы у = у* новый установившийся режим системы указывается неоднозначно.

Опасную границу будем называть динамически неопределенной, если по крайней мере две траектории, выходящие из и (Г *), имеют разные стоки. В случае распада гомоклинического цикла Г* новым установившимся режимом будет один из трех стоков, расположенных вблизи угловых точек (а, 0, 0), (0, а, 0) или (0, 0, а) куба поэтому опасная граница является динамически неопределенной.

Следовательно, если при у = у* изображающая точка находится в и (Г *), малой окрестности Г*, то из аксиомы неопределенности следует, что малые возмущения начальных данных могут приводить к качественному различию в поведении решений системы (1). Кроме того, при у = у* — <г, <г > 0, области притяжения предельных циклов становятся очень малыми при стремлении <г к нулю. Это означает, что случайные возмущения, присутствующие всегда, например, при численном интегрировании, будут выбрасывать изображающую точку из этой области еще до того, как область притяжения полностью исчезнет. Выбор стока, в который «перескакивает» изображающая точка в фазовом пространстве, в этой ситуации носит случайный характер.

Таким образом, мы описали ситуацию в фазовом пространстве, которая характеризуется высокой параметрической чувствительностью к начальным данным. При численном анализе системы (1) в окрестности значения у = у* мы будем наблюдать, что решения с близкими и даже почти равными начальными данными будут случайным образом притягиваться к существенно различным стокам.

Описанная ситуация имеет простой биологический смысл: малые флуктуа-ционные выбросы концентраций в то время, когда происходит генерация релаксационных импульсов в генной сети, могут сбивать систему с колебательного режима и случайным образом стабилизировать ее.

Отметим, что качественное описание максимальных семейств периодических решений системы (1) в зависимости от других параметров, 7 или а, аналогично приведенному выше. При увеличении значения параметра периодическое решение появляется в результате локальной бифуркации Андронова — Хопфа, затем амплитуда колебаний увеличивается, форма периодической орбиты в фазовом пространстве изменяется, период автоколебаний растет, и на другой границе существования автоколебаний максимальное семейство заканчива-

6 5 4

z

3

2

1 2

Рис. 10. Семейство периодических решений в зависимости от параметра а при y = 2.981, ß = 2 и 5.908 < а < 24.5.

Рис. 11. Семейство периодических решений в зависимости от параметра y при а = 5.908, ^ = 2 и 2.98 < Y < 4.087.

ется глобальной бифуркацией влипания периодической орбиты в треугольный гомоклинический цикл. Интервал существования периодических решений по параметру y сравним с интервалом существования автоколебаний по параметру у, однако интервал существования периодических решений по параметру а много больше, чем по параметрам y или у. На рис. 10, 11 показана эволюция предельного цикла при варьировании а и y.

Авторы выражают благодарность Р. Л. Гомбоеву за участие в проведении численных экспериментов.

ЛИТЕРАТУРА

1. Колчанов Н. А. Регуляция транскрипции генов эукариот: базы данных и компьютерный анализ // Молекулярная биология. 1997. № 31. С. 581—583.

2. Колчанов Н. А., Игнатьева Е. В., Подколодная О. А., Лихошвай В. А., Матушкин Ю. Г. Генные сети // Вавилов. журн. генетики и селекции. 2013. Т. 17, № 4/2. С. 833—850.

3. Golubyatnikov V., Likhoshvai V., Ratushny A. Oscillatory gene networks and Hopf bifurcation // 7th Int. Conf. Human and Computer-2004. 2004. P. 72-77.

4. Волокитин Е. П. О предельных циклах в простейшей модели гипотетической генной сети // Сиб. журн. индустр. математики. 2004. Т. 7, № 3. С. 58-65.

5. Xiao M., Cao J. Genetic oscillations deduced from Hopf bifurcation in a genetic regulatory network with delays // Math. Biosci. 2008. V. 215, N 1. P. 55-63.

6. Кастлер Г. Возникновение биологической организации. М.: Мир, 1967.

7. Гукенхаймер Дж., Холмс Ф. Нелинейные колебания, динамические системы и бифуркации векторных полей. М.; Ижевск: Ин-т компьютер. исслед., 2002.

8. Golubyatnikov V. P., Likhoshvai V. A., Fadeev S. I., Matushkin Yu. G., Ratushny A. V., Kolchanov N. A. Mathematical and computer modeling of genetic networks // 6th Int. Conf. Human and Computer—2003. University of Aizu, Japan. 2003. P. 200—205.

9. Kuznetsov Yu. A. Elements of applied bifurcation theory. New York: Springer-Verl., 2004.

10. Чумаков Г. А., Лашина Е. А., Чумакова Н. А. Максимальные семейства периодических решений кинетической модели гетерогенной каталитической реакции // Вестн. Новосиб. гос. ун-та. Сер. Математика, механика, информатика. 2005. Т. V, № 4. С. 3—20.

11. Chumakova N. A., Chumakova L. G., Kiseleva A. V., Chumakov G. A. Computation of periodic orbits in a three-dimensional kinetic model of catalytic hydrogen oxidation // Selcuk J. Appl. Math. 2002. V. 3, N 1. P. 3-20. (http://www5.in.tum.de/selcuk/sjam).

12. Chumakov G. A., Lashina E. A., Chumakova N. A. On estimation of the global error of numerical solution on canard-cycles // Math. Comput. Simul. 2014. http://dx.doi.org/10.1016/ j.matcom.2014. 10.003 (in press).

13. Баутин Н. Н., Шильников Л. П. Поведение динамических систем вблизи границ области устойчивости состояний равновесия и периодических движений («опасные» и «безопасные» границы) // Бифуркация рождения цикла и ее приложения. Дополнение I. М.: Мир, 1980. С. 294-316.

Статья поступила 30 октября 2014 г.

Чумаков Геннадий Александрович Институт математики им. С. Л. Соболева СО РАН, пр. Академика Коптюга, 4, Новосибирск 630090; Новосибирский гос. университет, ул. Пирогова, 2, Новосибирск 630090 chumakov@math.nsc.ru

Чумакова Наталия Алексеевна Институт катализа им. Г. К. Борескова СО РАН, пр. Акад. Лаврентьева 5, Новосибирск 630090; Новосибирский гос. университет, ул. Пирогова, 2, Новосибирск 630090 chum@catalysis.ru