Автоколебания в гипотетических генных сетях



Скачать 168.54 Kb.
Дата01.08.2016
Размер168.54 Kb.





АВТОКОЛЕБАНИЯ В ГИПОТЕТИЧЕСКИХ ГЕННЫХ СЕТЯХ

В.А.Лихошвай*2,3, В.В. Когай1, С.И.Фадеев1
1Институт математики Сибирского отделения Российской академии наук, Новосибирск, 630090.

2Институт цитологии и генетики Сибирского отделения Российской академии наук, Новосибирск, 630090.

3Югорский НИИ информационных технологий, Ханты-Мансийск, 628011
*контактный автор likho@bionet.nsc.ru
Keywords: genetic systems, modeling, delay equation, differential autonomous systems

Resume


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

Results: В работе излагается эффективный метод численного исследования автоколебаний, описываемых автономными системами дифференциальных уравнений специального вида, моделирующих симметричные канонические гипотетические генные сети первого класса (Likhoshvai et al, 2003).
Введение.

Генные сети (ГС) являются структурно сложными пространственными объектами, состоящими из десятков и сотен элементов разной природы и сложности: гены и их регуляторные участки; РНК и белки, кодируемые этими генами; низкомолекулярные соединения, различные комплексы между ферментами и их мишенями и т.д. Ядром ГС являются регуляторные контуры - гены и белки, экспрессия которых подвержена взаимному регулированию. Их наличие обеспечивает ГС уникальную способность адекватно реагировать на изменение внешних условий. Таким образом, выявление у регуляторных контуров возможных режимов функционирования является важной задачей теории генных сетей. Чтобы подойти к решению данной проблемы необходимо проведение систематического исследования закономерностей функционирования регуляторных контуров различного рода структур. Конструктивным шагом в этом направлении представляется выделение из природных ГС некоторого конечного набора стандартных элементов, формализации правил сборки из них теоретических объектов (математических моделей), описывающих регуляторные контуры, и проведение систематического анализа их свойств для выявления общих биологически значимых закономерностей [1].

В данной работе приводится эффективный метод поиска циклических режимов функционирования симметричных гипотетических генных сетей (ГГС) – наиболее простых по формулировке математических объектов, представляемых автономными системами из n дифференциальных уравнений:
(1 )
где положительные параметры, Правые части системы (1) представляют регуляторные связи первого класса. Легко показывается, что фазовые траектории системы (1), выходящие из точки, принадлежащей гиперкубу  с ребром , остаются в . В дальнейшем автономную систему (1) мы будем называть моделью M(n,k).

Подробное изложение свойств симметричных ГГС имеется в работе [2]. В частности, согласно -критерию, сформулированному для симметричных ГГС, существует и такие, что при и автономная система (1), имеет только асимптотически устойчивых стационарных решений, если наибольший общий делитель чисел и равен . Если , то существуют устойчивых предельных циклов, а устойчивые стационарные решения отсутствуют. Общее число стационарных решений модели М(n,k) (как устойчивых, так и неустойчивых) для достаточно больших  и  равно 2d-1.

Один из подходов численного исследования автоколебаний автономных систем общего вида в зависимости от параметров модели состоит в следующем [3-5]. Предположим, что при фиксированном значении параметра модели (например, фиксированном  в системе (1)) начальные данные задачи Коши для автономной системы уравнений таковы, что решение системы, начиная с некоторого момента времени, достаточно хорошо приближает установившийся автоколебательный режим с некоторым периодом. В этом случае возникает возможность уточнить описание автоколебаний, исходя из краевой задачи для рассматриваемой автономной системы с условиями периодичности и трансверсальности. При этом используется метод Ньютона, в котором за начальное приближение берется результат интегрирования задачи Коши. Полученное решение краевой задачи играет роль стартового в методе продолжения по параметру, позволяющему эффективно проводить численное исследование автоколебаний в зависимости от параметра.

Краевую задачу, описывающую автоколебания автономной системы (1) (если они существуют) можно представить в виде:


(2)

, (3)

при s = 0. (4)
Здесь T - период, подлежащий определению, равенство (3) – условие периодичности решения, равенство (4) – один из возможных вариантов условия трансверсальности, когда производная первого компонента при обязана быть равна нулю. Так как определение зависимости решения краевой задачи от параметра, описывающей автоколебания методом продолжения по параметру, не связано с проблемой устойчивости то метод позволяет находить как устойчивые, так и неустойчивые периодические решения, если последние имеют место в рассматриваемой математической модели. Устойчивость можно в последствии проверить, например, используя алгоритм вычисления наибольшего собственного числа матрицы монодромии.

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


Численное исследование автоколебаний автономных систем.

Краевая задача (2)-(4) является частным случаем более общей проблемы, поиска периодической вектор-функции x(s), x(s)=x(s+1), с компонентами , дающую решение краевой задачи в зависимости от скалярного параметра q автономной системы уравнений с достаточно гладкими правыми частями в некоторой области их определения:



(5)

, . (6)

Здесь f(x,q) – вектор-функция векторного аргумента x и скалярного аргумента q с компонентами . Если обозначить через y и F составные векторы



,
то краевой задаче (5) можно придать “стандартный” вид:
, , (7)
где вектор-функция представляет краевые условия (6).

Пусть y = y(s,q) – решение краевой задачи (7) при некотором значении параметра q. Для численного определения y(s,q) воспользуемся методом “множественной стрельбы” [3,5]. Согласно методу “множественной стрельбы” разобьем отрезок [0,1] по s на m частей:


0 = s1 < s2 < … < sm < sm+1 = 1. (8)

Обозначим через pk сеточное значение вектор-функции y(s,q) в k – ом узле сетки (8): pk = y(sk,q), k = 1,2,…. m+1. На каждом из отрезков [sk, sk+1], k = 1,2,…. m+1, рассмотрим следующие серии задач Коши в виде векторного и матричного уравнений:



, (9)

(10)

(11)
где - единичная матрица. Решения серии (9) на каждом из отрезков [sk, sk+1] , обозначим через y(s,pk,q), решения серии (10) – через Y(s,pk,q), решения серии (11) – через u(s, pk,q). Очевидно, серия (9) будет давать решение краевой задачи (7) при выполнении краевых условий и условий непрерывности вектор-функции y(s,q) в узлах сетки:
Ф1 = g(p1,pm+1,q) = 0,

Ф2 = y(s2,p1,q) - p2 = 0

………

Фm+1 = y(sm+1,pm,q) - pm+1 = 0.



или

Ф(p,q) = 0, (12)


где p и Ф – составные векторы размеров N, N = (n + 1)(m+1). Серии задач Коши (10) и (11) позволяют выписать (N x N)-матрицу производных Фp вектор-функции Ф(p,q) и столбец производных Фq по параметру q:
Фp=,
Фq =.
Здесь для удобства используются обозначения векторных аргументов вектор-функции g(y(0),y(1),q) : a = y(0), b = y(1).

Таким образом, формально проблема численного исследования решения краевой задачи (7) в методе “множественной стрельбы” сводится к численному исследованию решения системы нелинейных уравнений (12) с параметром q, определяющей вектор-функцию p = p(q). Организация процесса продолжения по параметру изложена в работах [6-8] . Число разбиений отрезка [0,1] по s определяется “приемлемой жесткостью” задач Коши (9)-(11), которая достигается за счет выбора величины отрезка [sk, sk+1], например, из условия: s[si, si+1], max  D, D( si+1- si ) 1. При этом столбцы матрицанта Y, удовлетворяющего условиям (10) на итерациях по методу Ньютона, остаются почти ортогональными при s = si+1, в полном соответствии с идеей метода ортогональной прогонки С.К.Годунова [9]. Для эффективности метода в случае больших N существенен учет при вычислениях структуры матрицы Фz.


Уравнения с запаздывающими аргументами, описывающие автоколебания модели M(n,k).

Как показывают численные эксперименты, периодические решения краевой задачи (5),(6) проведенные для разных значений n и k, обладают свойством частичной симметрии, которое выражается в том, что все переменные разбиваются на d равномощных групп, в которых траектории идентичны с точностью до сдвига по фазе. При этом d может принимать значения 1,2,3... Если d=1, то цикл обладает полной симметрией, и называется симметричным. В общем случае цикл будем называть d-симметричным. Также, численно установлено, что для устойчивых циклов d=нод(n,k).

Свойство частичной симметрии дает возможность описывать автоколебания автономной системы (1) краевой задачей для системы из d уравнений с запаздывающими аргументами, векторное представление которой имеет вид:
(13)

, , .
где - вектор запаздываний с компонентами , задаваемый в долях периода T, а определяется правыми частями автономной системы (1).

Дискретная модель краевой задачи (13) в виде системы нелинейных уравнений с параметром q может быть построена способом, аналогичным [10].


Примеры численного исследования автоколебаний

Приведем некоторые результаты численного исследования автоколебаний модели M(n,k). На рис.1 представлено решение краевой задачи (5),(6) модели M(6,4) при указанных значениях параметров и периодом колебаний T = 5.75. Каждая из групп характеризуется своей амплитудой. На рисунке указана последовательность сдвига кривых в каждой группе. При этом сдвиг соседних графиков в каждой группе переменных равен T/3. Эквивалентная рассматриваемой краевой задаче система из двух уравнений с запаздывающим аргументом имеет вид:


,

(14)


Здесь u1(s) = x1(s), u2(s) = x2(s). Тогда, согласно сдвигу графиков компонент на рис.1, x3(s) = u1(s-2/3), x4(s) = u2(s-2/3), x5(s) = u1(s-1/3), x6(s) = u2(s-1/3).

Другой тип автоколебаний той же модели представлен на рис.2. Здесь все компоненты имеют одну и ту же амплитуду. Сдвиг соседних графиков равен 1/6 периода. Соответствующая краевая задача для уравнения с запаздывающим аргументом имеет вид:


u(s) = x1(s),

(15)







Рис.1. 2-симметричные автоколебания модели М(6,4).

Рис.2. Симметричные автоколебания модели М(6,4).

На рис.3 приведены в виде графиков результаты численного исследования автоколебаний модели M(6,4) в зависимости от параметра , полученные методом продолжения по параметру. Графики с отметкой (1) представляют ветви зависимости периода T и амплитуды A колебаний компонентов двух групп. Как следует из рисунка, за областью изменения  от 0 до точки поворота при , в которой автоколебания отсутствуют, следует область , в которой при одном и том же значении  существуют два предельных цикла (отметки (1), (2)). При возникает предельный цикл с одинаковыми по амплитуде компонентами (отметка (3)). Значение, при котором амплитуды равны нулю, соответствует особой точке на диаграмме стационарных решений модели M(6,4). В особой точке пересекаются симметричное решение, т.е. решение, имеющее одинаковые компоненты, и два частично симметричных решений, у которых свои одинаковые значения имеют четные и нечетные компоненты. При этом все стационарные решения неустойчивы, если , что объясняет возникновение автоколебаний модели M(6,4) при достаточно больших значениях .

В качестве еще одного примера рассмотрим автоколебания модели M(5,3). Согласно -критерия, здесь следует ожидать при достаточно больших значениях параметров  и  один устойчивый симметричный предельный цикл. Отметим, что модель M(5,3) имеет только один стационар. Начиная с некоторого значения   0, стационарное решение теряет устойчивость, переходя в автоколебания.







Рис. 3. Зависимость периода T и амплитуды автоколебаний A модели M(6,4) от параметра



На рис.4-5 представлены два численно выявленных предельных

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






Рис.4 Симметричные автоколебания модели М(5,3), характеризующиеся одинаковыми порядками следования переменных по фазе и периодами, но разными амплитудами.

Рис.5. Другая форма симметричного автоколебания модели М(5,3) при тех же значениях параметров, что и на рис.4.

Автоколебания в системе возникают при  > 1=2.44 и не связаны с бифуркацией Андронова-Хопфа(см. рис 6).

Кривым, указанным на рис.4, соответствует уравнение:
(16)
Случаю (рис.5), соответствует уравнение
(17)







Рис.6. Зависимость периода T и амплитуды колебаний A модели M(5,3) от параметра .

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

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

Очевидно, проведение исследования автоколебаний модели M(n,k) методом продолжения по параметру решения краевой задачи для системы уравнений с запаздывающими аргументами представляется более эффективным подходом, поскольку в этом случае проблема сводится к d уравнениям, d
Определение устойчивости автоколебаний.

Остановимся на некоторых известных фактах в связи с проблемой устойчивости автоколебаний [11,12]. Рассмотрим автономную систему уравнений



(18)
которая при некотором значении параметра q имеет T-периодическое решение, описываемое вектор-функцией x = x(t,q), x(t,q) = x(t+T,q), Обозначим через A(t) матрицу производных правых частей системы (1), , а через U(t) – матрицант:

Матрица M = U(T) называется матрицей монодромии, а ее собственные числа – мультипликаторами.

Как известно, матрица монодромии имеет мультипликатор, равный 1, которому соответствует собственный вектор v(T). Согласно теореме Андронова-Витта [11], если мультипликатор, равный 1, имеет кратность 1, а все остальные мультипликаторы лежат внутри единичной окружности на комплексной плоскости, то периодическое решение системы (21) является устойчивым по Ляпунову. Если хотя бы один из мультипликаторов лежит вне единичной окружности, то соответствующее периодическое решение неустойчиво. Следовательно, определение устойчивости периодического решения сводится к нахождению наибольшего по абсолютному значению собственного числа матрицы монодромии [11,12], либо к дихотомии спектра матрицы монодромии окружностью единичного радиуса [13]. Так как у матрицы монодромии всегда есть одно собственное число, равное 1, то вначале необходимо применить к ней метод исчерпывания [14,15], выделив единичное собственное число, с целью получить матрицу размерности на единицу меньше, спектр которой уже не содержит выделенного собственного числа.

Метод исчерпывания состоит в следующем. Пусть A  (nxn)-матрица с элементами aij, i,j = 1,2,…,n, а  - одно из её собственных чисел , которому соответствует собственный вектор z с компонентами z1, z2,…,zn. При этом будем считать, что компонент zk вектора z отлична от 0. Тогда ((n-1)x(n-1)) – матрица B с элементами bij,


имеет те же собственные числа, что и матрица A за исключением выделенного собственного числа .

Как уже отмечалось, в матрице монодромии M выделенным является собственное число, равное 1, которому соответствует, согласно (22), собственный вектор v(T). Обозначим через L матрицу, полученную в результате применения метода исчерпывания к матрице M. Для вычисления наибольшего собственного числа max матрицы L воспользуемся степенным методом, предполагая, что max – вещественное [16-18]. В степенном методе по заданному начальному вектору u строится последовательность векторов с компонентами При этом для достаточно больших значений формула


,

дает приближенное значение max, которое мало зависит от номера .





Рис.7. Графики наибольшего по модулю собственного числа  матриц монодромии модели M(6,4) в зависимости от параметра .
На рис.7 приведены результаты исследования устойчивости периодических решений модели M(6,4) в виде графиков зависимости нибольшего по модулю собственного числа (обозначение ) матрицы монодромии от параметра . Неустойчивым автоколебаниям соответствуют знчения Из сопоставления рис.3 и рис.7 следует, что ветвь с отметкой (1) представляет устойчивые, а ветвь с отметкой (2) – неустойчивые автоколебания. Автоклебания с отметкой (3) также являются неустойчивыми.

На рис.8 приведены графики зависимости наибольшего по модулю собственного числа  матрицы монодромии от параметра ,. соответствующие модели M(5,3). Из сопоставления рис.6 и рис.8. следует, что автоколебания, устойчивые для достаточно больших значений , сохраняют устойчивость с уменьшением  вплоть до точки поворота при . Далее, после точки поворота при автоколебания становятся неустойчивыми. В области автоколебания вновь устойчивы, а при автоколебания теряют устойчивость.

Отметим другой подход к исследованию устойчивости периодических решений автономной системы в зависимости от параметра, который связан с определением нормы решения дискретного матричного уравнения Ляпунова [13].



Рис.8. Графики наибольшего по модулю собственного числа  матриц монодромии модели M(5,3) в зависимости от параметра .

Заключение.

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



Для примера рассмотрим уравнение
. (19)
Из результатов, представленных выше, следует, что данное уравнение при определенных значениях параметров , , а,  имеет траектории, идентичные циклическим траекториям системы (1) при k=3. В то же время при =50, =50, а=25, =0.031 численный расчет демонстрирует следующее поведение решения, представленное на рис.9.

Рис. 9. График решения уравнения (26) при =50, =50, а=25, =0.031 и . По оси ординат отложено значение переменной t, по оси ординат – значение функции x в точке t.
Предварительный анализ показывает, что траектория не является периодической или квазипериодической и, возможно, мы имеем странный аттрактор. Поскольку уравнение (19) может интерпретироваться также как модель, описывающая биосинтез белка с двумя отрицательными авторегуляциями, то можно заключить о возможности хаотического синтеза белка уже в простейших биологических системах. Практическое применение полученные результаты могут иметь в области синтеза генных сетей с заданными свойствами, а также могут быть использованы при решении практических задач в областях биотехнологии, биотерапии, генной инженерии и фармакогенетики.

Благодарности

    Работа была частично поддержана РФФИ (гранты No 02-04-48802, 02-07-90359, 03-01-00328, 03-04-48506, 03-07-96833, 03-04-48829, 04-01-00458) и СО РАН (Интеграционный проект No. 119), Программой фундаментальных исследований Президиума РАН и отделений РАН (№10.4), государственным контрактом 10002-251/П-25/155-270/200404-082.


Литература.
1. Лихошвай В.А., Матушкин Ю.Г., Фадеев С.И. (2003) Задачи теории функционирования генных сетей. Сибирский журнал индустриальной математики. 4, 64-80.
2. Лихошвай В.А., Фадеев С.И. (2003) О гипотетических генных сетях. Сибирский журнал индустриальной математики. 4, выпуск 3(15), стр 134-153. Новосибирск, ИМ СО РАН.
3. Когай В. В., Фадеев С. И. (2001) Применение продолжения по параметру на основе метода множественной стрельбы для численного исследования нелинейных краевых задач. Сибирский журнал индустриальной математики. 4, С.83-101.
4. Когай В. В. (2002) Применение продолжения по параметру для численного исследования периодических решений автономных систем обыкновенных дифференциальных уравнений // Вестник НГУ, 2, С.40-48.
5. Fadeev S. I., Kogai V. V.Using parameter continuation based on multiple shooting method for numerical research of nonlinear boundary value problems // International J. of Pure and Appl. Mathematics. Vol.14. \No.4. 2004. 467-498p.
6. Фадеев С.И. О решении системы трансцендентных уравнений с параметром методом Ньютона. // Вычислительные системы. – Новосибирск, 1985.- Вып.108. – С.78-93

7.Холодниок М.,Клич А., Кубичек М., Марек М.: Методы анализа нелинейных динамических моделей. М.: Мир, 1991. 320 с.


8. Фадеев С. И., Покровская С. А., Березин А. Ю., Гайнова И. А.. Пакет программ STEP для численного исследования систем нелинейных уравнений и автономных систем общего вида. Новосибирск: Изд-во НГУ. 1998. 188с.
9. Годунов С.К. Обыкновенные дифференциальные уравнения с постоянными коэффициентами. Т.~1. Краевые задачи. Новосибирск: Изд-во НГУ - 1994. - 264с.
10 Фадеев С.И. Программа численного решения нелинейных краевых задач для систем обыкновенных дифференциальных уравнений с параметром. // Вычислительные методы линейной алгебры.- Новосибирск. – Наука. – Сибирское отделение.- 1990, С.104-198
11. Бибиков Ю. Н. Курс обыкновенных дифференциальных уравнений. М.: Высшая школа. 1991. 304с.
12. Понтрягин Л. С. Обыкновенные дифференциальные уравнения. М.: ГИФМЛ. 1961. 312с.
13. Годунов С. К. Современные аспекты линейной алгебры. Новосибирск: Научная книга. 1997. XXVI+390с.
14. Коллатц Л. Задачи на собственные значения. М.: Наука. 1968. 504с.
15. Годунов С. К., Антонов А. Г., Кирилюк О. П., Костин В. И. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. Новосибирск: Наука. Сиб. отд-ние. - 1988. 456с.
16. Семендяев К. А. О нахождении собственных значений и инвариантных многообразий матриц посредством итераций // ПММ. Т.7. N.3 1943. С.193-222.
17. Фаддеев Д. К., Фаддеева В. Н. Вычислительные методы линейной алгебры. М.: ГИФМЛ. 1960. 656с.
18. Крылов В.И., Бобков В.В., Монастырный П.И. . Вычислительные методы высшей математики. Минск. Издательство Высшая школа. 1972. 584 с.




Каталог: hgnet -> stat


Поделитесь с Вашими друзьями:


База данных защищена авторским правом ©uverenniy.ru 2019
обратиться к администрации

    Главная страница