О численном исследовании автоколебаний в гипотетических генных сетях Текст научной статьи по специальности «Математика»

Вычислительные технологии Том 10, № 3, 2005

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

В. В. КогАй, С. И. Фадеев Институт математики им. С.Л. Соболева СО РАН,

Новосибирск, Россия e-mail [email protected], [email protected]

В. А. ЛихошВАй Институт цитологии и генетики СО РАН, Новосибирск, Россия Югорский НИИ информационных технологий, Ханты-Мансийск, Россия

e-mail [email protected]

Identification of cycle operations in gene networks is an important problem of bioinformatics that analyzes structure and functioning regularities of gene networks and their regulator circuits. In the article we formulate an effective method of calculation of the auto-oscillations governed by differential autonomous systems of special form which model symmetric canonical hypothetical gene networks of the first class.

Введение

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

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

дение систематического анализа их свойств для выявления общих биологически значимых закономерностей [1].

В настоящей работе описывается эффективный метод поиска циклических режимов функционирования симметричных гипотетических генных сетей (ГГС) — наиболее простых по формулировке математических объектов, представляемых автономными системами из п дифференциальных уравнений:

<іхі а , ,

и = т+ж- х‘- г=1'2...............п’ (1)

где

¿1 = хП + ХП—1 + ••• + Х1-к+3 + Х1-к+2>

¿2 = Х1 + ХП + ••• + Х1-к+4 + ХП-к+3>

г -- х^ + х^ + + х^ + х^ •

¿п п— 1 + хп— 2 + ••• + хп—к+2 + хп— к+1;

а, в, 7 — положительные параметры; 1 < к < п. Правые части системы (1) представляют регуляторные связи первого класса. Легко показывается, что фазовые траектории системы (1), выходящие из точки, принадлежащей гиперкубу 2 с ребром а, остаются в 2. В дальнейшем автономную систему (1) мы будем называть моделью М1 (п,к).

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

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

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

dxi ^ / а \ dT

0<8<1 & — т(г+й; — х;) ■ и — ° г — 12’-’п; (2)

Хі(0) = Хі(1);

(3)

а

——--------Х1 - 0 при 8 - 0. (4)

1 + в^1

Здесь Т — период, подлежащий определению; равенство (3) — условие периодичности решения; равенство (4) — один из возможных вариантов условия трансверсальности, когда производная первой компоненты при 8 — 0 обязана быть равна нулю.

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

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

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

dx ^ . dT

0 < 8 < l, — — Т^(х^ — 0;

¿Х1 , .

х(0) = х(1), —— (0) = 0^

ав

(5)

(6)

Здесь / (х, д) — вектор-функция векторного аргумента х и скалярного аргумента д с компонентами /1 (х, д), f2(x, д), • • • , /„(х, д). Если обозначить через у и ^ составные векторы

Х , ^ = -Т/-

Т 0

то краевой задаче (5) можно придать “стандартный” вид:

0 < в < 1, ^ ^(у,д), ^(у(0),у(1),д) = 0,

(7)

где вектор-функция д(у(0), у(1), д) представляет краевые условия (6).

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

0 — 81 <82 < . . . < 8т < 8т+1 — 1.

Обозначим через рк сеточное значение вектор-функции у(8,д) в к-м узле сетки (8)

Р

у(вк ,д), к = 1, + 1

у

На каждом из отрезков [вк, «к+1], к = 1, 2, • • • , т, рассмотрим следующие серии задач Коши в виде векторного и матричного уравнений:

« є [вк,вк+1], оОв = р(у,д)’ у(^)

Р

От = (у.д)« + (у.д). «(«к-) = о,

(9)

(10)

(11)

где I — единичная матрица. Решения серии (9) на каждом из отрезков [«к, «к+1 ] обозначим через у(т,рк, д), решения серии (10) — через У(т,рк,д), решения серии (11) — через «(т,рк,д). Очевидно, серия (9) будет давать решение краевой задачи (7) при выполнении краевых условий и условий непрерывности вектор-функции у(т, д) в узлах сетки:

Ф1 = ^(р1,рт+1,д) = 0,

Ф2 = у(«2,р1, д) - р2 = 0,

Фз = у(«з,р2,д) - р3 = 0,

или

Фт+1 = у(«т+1,Рт,д) - Рт+1

ф(р, д) = 0,

(12)

где р и Ф — составные векторы размера N, N — (п + 1)(т + 1). Серии задач Коши (10) и (11) позволяют выписать ^ х N)-матрицу производных Фр вектор-функции Ф(р,д) и столбец производных Фд (р, д) по параметру д:

д«(р1,рт+1,д) 0 • •0 де (р1,рт+1,д)

У(«2 ,р1,д) -/ • •0 0

фр(р,д) = 0 У(«з,р2,д) • •0 0

0 0 • • У («т+1,рт,д) -I

ФР(Р,д)

д, (р1,рт+1,д) «(«2,р1, д)

«(тз,р2,д)

м(«т+1,Рт ,д)

Здесь для удобства используются обозначения векторных аргументов вектор-функции

д(у(0),у(1),д): а = у(0) в = у(1).

Таким образом, формально проблема численного исследования решения краевой задачи (7) в методе “множественной стрельбы” сводится к численному исследованию решения системы нелинейных уравнений (12) с параметром д, определяющим вектор-функцию Р = р(д). Здесь для фиксированного значения параметра “пристреливаются” (т +1) сеточных значений вектор-функции у = у(т,д), причем процесс пристрелки организуется в виде итераций по методу Ньютона, что требует на каждой из итераций определения вектора Ф и матрицы Б = [Фр, Фд] из серии задач Коши (9)-(11). Кроме того, привлечение

0

метода продолжения по параметру для построения зависимости р — р(д) и, следовательно, зависимости решения краевой задачи (7) от параметра д предусматривает вычисление производных решения системы (12) по параметру, т. е. обращения к матрице Б после завершения итераций по Ньютону.

Согласно теореме о неявных функциях, система (12) определяет в ^ + 1)-мерном пространстве (р, д) гладкую пространственную кривую Е — график вектор-функции р — р(д), если в окрестности Е ранг матрицы Б равен N. В этом случае численное построение Е методом продолжения по параметру основано на равноправии аргументов вектор-функции Ф(р,д): любая из компонент вектора р и д может играть роль параметра ^ системы (12), а оставшиеся N аргументов составляют вектор неизвестных г, если в некоторой окрестности изменения ^ соответствующая матрица производных Фг будет невырожденной.

Остановимся на организации процесса продолжения по параметру [6, 7, 8]. При регулярном выборе параметра ^ (регуляризации), который мы будем называть текущим параметром, всякий раз решается система нелинейных уравнений

Ф(г,^) —0, det Фг — 0. (13)

Методом Ньютона находятся решение г — г(^) системы (13) (в частности, д — д(^), если д не является текущим параметром), а также производная решения гм(^) по параметру ^ из системы линейных алгебраических уравнений

Фг (г,^)^ — -Фм, (14)

после того как найдено решение системы (13). Далее, среди компонент вектора г и д определяется очередной текущий параметр V, который соответствует максимальной по модулю компоненте вектора гм, нормировкой находится производная решения (13) уже по параметру V. Таким образом, при найденном в результате параметризации значении параметра V известны как решение ¿(V), так и производная (V). Это позволяет, сде-

лав достаточно малый шаг Дv по параметру V, задать начальное приближение решения системы нелинейных уравнений со значением текущего параметра V + Д V. Величина шага Дv автоматически регулируется, например, ограничением на число итераций в методе Ньютона, за которое норма невязок приближения должна быть меньше заданного числа. В стартовой позиции требуется как задать начальное приближение, так и указать текущий параметр среди компонент вектора р и д.

Отметим, что значениям параметра д, при которых det Фр — 0, отвечает “точка поворота” пространственной кривой Е. В точке поворота выполняется условие: производная параметра д по текущему параметру ^ равна нулю. (Очевидно, что в этом случае д не является текущим параметром.) Таким образом, в точке поворота пространственная кривая имеет вертикальную касательную, что означает либо бифуркацию числа решений системы (12) и, следовательно, решений краевой задачи (7), либо точку перегиба.

Число разбиений отрезка [0,1] по в определяется “приемлемой жесткостью” задач Коши (9)—(11), которая достигается за счет выбора величины отрезка [вд, вд+1], например, из условия:

в е [в*,вт], тах ||(у,д)|| < Д Д(вт - в*) ^ 1.

При этом столбцы матрицанта У, удовлетворяющего условиям (10) на итерациях по методу Ньютона, остаются почти ортогональными при в — в*+1 в полном соответствии с идеей метода ортогональной прогонки С. К. Годунова [9]. Для эффективности метода в случае больших N существен учет при вычислениях структуры матрицы Фг.

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

Как показывают численные эксперименты, периодические решения краевой задачи (5), (6) обладают следующим свойством. Пусть О — наибольший общий делитель чисел п и к, О = к. Тогда компоненты вектор-функции ж(з), являющейся решением краевой задачи (5),

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

О«

0 < « < 1, — = ТС(и(«),«(« - тА «(« - т2), • • •,«(« - тк-1), д),

* О« (15)

«(0) = «(1), —— (0) = 0,

О«

где т — вектор запаздываний с компонентами т1, т2, • • • , Тк-1, задаваемый в долях периода Т; вектор-функция С(и(з),«(« - т^1), • • • ,«(« - Тк-1 ),д) определяется правыми частями автономной системы (1). В силу периодичности решения для каждой компоненты т) вектора т должно выполняться условие

-т) < « < 0, «(«) = «(1 + «), (16)

что позволяет рассматривать решение краевой задачи только на отрезке [0, 1] по «.

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

0 < « < 1, Я(и(з), д) = С(и(з),«(« - т}),«(« - т2), • • •,«(« - тк-1), д)•

Тогда на заданном разбиении (8) отрезка [0,1] по « краевую задачу (15) можно представить в виде

«к + 1

м(«к+1) - «(«к )= J ^«(з^д)^, к = 1, 2,^^^,т,

«к

«(0) = «(1), —^ (0) = 0^

О«

Заменим интеграл на приближенную квадратурную формулу, в которой для вычисления промежуточных значений подынтегральной функции воспользуемся приближенным представлением «(«) локальным сплайном. Используя далее условие периодичности решения, мы приходим к дискретной модели краевой задачи в виде системы нелинейных уравнений с параметром д относительно сеточных значений вектора неизвестных на отрезке [0,1] по «.

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

Приведем некоторые результаты численного исследования автоколебаний модели М1(п, к). На рис. 1 представлено решение краевой задачи (5), (6) модели М1(6, 4) при указанных

Рис. 1. Автоколебания модели М\(6, 4), в которых компоненты решения х(в)

образуют две группы.

значениях параметров с периодом колебаний Т = 6.21. Каждая из групп характеризуется своей амплитудой. На рисунке указана последовательность компонент в каждой из групп. Сдвиг соседних графиков компонент в каждой из групп равен Т/3. Эквивалентная рассматриваемой краевой задаче система из двух уравнений с запаздывающим аргументом имеет вид

0 < « < 1,

-и1 = Т г________________________а_____________________________( )

-з 1+ в (и2 (з — 1/3)+ и( (в — 1/3)+ и1 (з — 2/3)) 1 ’

-и2 = Тг________________________________________________________________________а_( )] (17)

-з 1+ в (и! (з) + и2 (з — 1/3)+ и? (з — 1/3)) 2 ’

и1(0) = и1(1), и2(0) = и2 (1) , (0) = 0.

-з

Здесь и1(з) — обозначение одной из компонент первой группы, например и1(з) = х1(з); и2(з) — обозначение одной из компонент второй группы, например и2(з) = х2(з). Тогда, согласно сдвигу графиков компонент на рис. 1, запишем

х3(з) = и1(з — 2/3), х4(з) = и2(з — 2/3), х5(з) = и1(з — 1/3), х6(з) = и2(з — 1/3).

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

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

0 < з < 1,

-и = Тг________________________________________________________а_( )]

-з г1+ в (и7 (з — 1/6)+ и^ (з — 1/3)+ и^ (з — 1/2)) (18)

и(0) = и(1), —1 (0) = 0.

-з

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

составляют одну группу.

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

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

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

ная с некоторого значения а > 0, стационарное решение теряет устойчивость, переходя в автоколебания. На рис. 4, 5 представлены предельные циклы, каждый из которых характеризуется своей амплитудой и периодом, а также порядком следования компонент. Автоколебания возникают при а > 5.74.

Последовательности компонент, указанных на рис. 4, соответствует краевая задача, которая имеет вид п(в) = ж1(з),

0 < в < 1,

^П = Тг____________________а____________________( )]

¿в г1 + в(п7(в - 1/5)+ п^(в - 2/5)) п(в)1, (19)

п(0) = п(1), --1 (0) = 0.

¿в

Другому случаю (рис. 5) соответствует краевая задача

0 < в < 1,

¿и = Т г___________________а____________________(

¿в [1 + в(п7(в - 1/5)+ п^(в - 3/5)) п(в)1, (20)

п(0) = п(1), --1 (0) = 0.

¿в

Рис. 4. Автоколебания модели М1(5, 3). Компоненты решения х(з) принадлежат одной группе.

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

ю-

9-

8

7-

6-

5-

4-

3-

Э=1

У = 5

2-

1-

а, = 2.45 а2= 5.74 а3= 6.52

0 2 а, 4 агбаз 8 10 12 14

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

На рис. 6 представлены результаты численного исследования автоколебаний модели М1(5, 3) в виде графиков зависимостей периода Т и амплитуды А от параметра а. Отметим явление удвоения периода при а = а2 = 6.52.

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

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

4. Определение устойчивости автоколебаний

Остановимся на некоторых известных фактах в связи с проблемой устойчивости автоколебаний [11, 12]. Пусть А(£) — непрерывная Т-периодическая (пхп)-матрица, А(£) = А(£+Т), а и(¿) — матрицант, т. е. и(¿) является решением задачи Коши однородного матричного уравнения

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

Рассмотрим автономную систему уравнений

— = А(()и, и (0) = /.

которая при некотором значении параметра д имеет Т-периодическое решение, описываемое вектор-функцией х = х(£,д), х(£,д) = х(£ + Т, д). Обозначим через А(£) матрицу производных правых частей системы (1), определенную на решении х(£,д):

А(£) = /х(х(^,д),д).

Продифференцировав по £ тождество

¿х

«(*) = ^ М) = / (х(*,д),д),

получаем

¿г (£) = /х(х(^,д),д)^(^) = А(Ф(£).

Таким образом, производную решения (21) по £ можно рассматривать как решение векторного однородного уравнения с матрицей А(£), следовательно,

•и(£) = и (¿)г>(0).

Поскольку г>(£) — Т-периодическая вектор-функция,

■и(Т) = и (Т )г>(0) = Ми(Т). (22)

Тем самым показано, что матрица монодромии М имеет мультипликатор, равный 1, которому соответствует собственый вектор г>(Т).

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

Метод исчерпывания состоит в следующем. Пусть А — (п х п)-матрица с элементами а-,, г,3 = 1, 2,... , п, а А — одно из ее собственных чисел, которому соответствует собственный вектор г с компонентами г1,г2,... , гп. При этом будем считать, что компонента гк вектора г отлична от нуля. Тогда ((п — 1) х (п — 1))-матрица В с элементами

г-

Ь-, = а,----- аку, г,3 = 1, 2,...,п, г = к, 3 = к,

гк

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

Как уже отмечалось, в матрице монодромии М выделенным является собственное число, равное 1, которому соответствует согласно (22) собственный вектор г>(Т). Обозначим

через Ь матрицу, полученную в результате применения метода исчерпывания к матрице М. Для вычисления наибольшего собственного числа Атах матрицы Ь воспользуемся степенным методом, предполагая, что Атах — вещественное [16, 17, 18]. В степенном методе по заданному начальному вектору и строится последовательность векторов Ьки, к = 1, 2,... , с компонентами (Ьки),, I = 1, 2,... , п — 1. При этом для достаточно больших значений к формула

Ак = (Ьк+1и)* к = 1 2

Атах (Ьки), , к 12,"',

тах

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

На рис. 7 приведены результаты исследования устойчивости периодических решений модели Мі (6, 4) в виде графиков зависимости наибольшего по модулю собственного числа (обозначение А) матрицы монодромии от параметра а. Неустойчивым автоколебаниям соответствуют значения |А| > 1. Из сопоставления рис. 3 и 7 следует, что ветвь с отмет-

Рис. 7. Графики наибольшего по модулю собственного числа Л матриц монодромии модели М1(6, 4), характеризующих устойчивость автоколебаний в зависимости от параметра а.

1 2 3 4 5 6 7 а

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

кой (1) представляет устойчивые, а ветвь с отметкой (2) — неустойчивые автоколебания. Автоколебания с отметкой (3) также являются неустойчивыми.

На рис. 8 приведены графики зависимости наибольшего по модулю собственного числа А матрицы монодромии от параметра а, соответствующие модели М1(5, 3). Из сопоставления рис. 6 и 8 следует, что автоколебания, устойчивые для достаточно больших значений а, сохраняют устойчивость с уменьшением а вплоть до точки поворота при а = а1. Далее, после точки поворота при а1 < а < а2 автоколебания становятся неустойчивыми. В области а2 < а < а3 они вновь устойчивы, а при а > а3 — теряют устойчивость.

Отметим другой подход к исследованию устойчивости периодических решений автономной системы в зависимости от параметра, который связан с определением нормы решения дискретного матричного уравнения Ляпунова [13]. Согласно теореме Ляпунова, если все собственные числа (п х п)-матрицы А расположены внутри круга на комплексной плоскости, то матричное уравнение

Н — А*ЯА = С, (23)

известное как дискретное матричное уравнение Ляпунова, однозначно разрешимо при любой матрице С. При этом, если С = С * > 0 — эрмитова положительно определенная матрица, то и решение Н будет также эрмитовой положительно определенной матрицей. Верно и обратное утверждение. Если Н = Н* > 0 — решение матричного уравнения Ляпунова, где С = С * > 0, то все собственные числа матрицы А расположены внутри круга единичного радиуса. При этом решение (23) можно представить в виде суммы бесконечного ряда:

Н = С + А*СА +(А*)2СА2 + ... + (А*)к САк + ... (24)

Обычно полагают С = I.

С приближением собственных чисел матрицы А к окружности единичного радиуса норма решения матричного уравнения (23) будет неограниченно расти. Считается, что если в процессе вычисления суммы ряда (24) норма частичной суммы окажется больше наперед заданного достаточно большого числа, то матричное уравнение неразрешимо с вычислительной точки зрения.

Таким образом, для определения устойчивости периодического решения требуется найти матрицу Ь, применив к матрице монодромии метод исчерпывания. Затем вычисляется норма решения дискретного матричного уравнения Ляпунова

X — Ь*ХЬ = I (25)

с использованием для этой цели представления решения (25) в виде бесконечного ряда:

X = I + Ь*Ь + (Ь*)2Ь2 + ... + (Ь*)к Ьк + ...

Периодическое решение будет устойчивым, если ||X|| < ф, где Q — достаточно большое число.

Заключение

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

1,6 —|

I------1------1-----1------1------1------1-----1------1------1-----1------1------г

100 150 200 250 300 350 400 450

500

Рис. 9. График решения уравнения (26) при а = 50, 7 = 50, а = 25, т = 0.031 и х(£) = 0.9, ^ < 0. По оси абсцисс отложено значение переменной по оси ординат — значение функции х в точке ¿.

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

Для примера рассмотрим уравнение

^ж

а

1 + ж7 (£ — т) + ж7 (£ — ат)

— ж.

(26)

Как уже отмечалось, данное уравнение при определенных значениях параметров а, 7, а, т имеет траектории, идентичные циклическим траекториям системы (1) при к = 3. В то же время при а = 50, 7 = 50, а = 25, т = 0.031 численный расчет демонстрирует колебания сложной формы, представленные на рис. 9. Предварительный анализ показывает, что, видимо, имеет место странный аттрактор. Поскольку одной из интерпретаций уравнения (26) является описание биосинтеза белка с двумя отрицательными авторегуляциями, не исключается возможность хаотического синтеза белка уже в простейших биологических системах.

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

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

Список литературы

[1] ЛихошвАй В.А., Матушкин Ю.Г., Фадеев С.И. Задачи теории функционирования генных сетей // Сиб. журн. индустриальной математики. Новосибирск: ИМ СО РАН. 2003. T. 6, № 2(14). C. 64-80.

[2] ЛихошвАй В.А., Фадеев С.И. О гипотетических генных сетях // Сиб. журн. индустриальной математики. Новосибирск: ИМ СО РАН. 2003. T. 6, № 3(15). C. 134-153.

[3] КогАй В.В., Фадеев С.И. Применение продолжения по параметру на основе метода множественной стрельбы для численного исследования нелинейных краевых задач // Сиб. журн. индустриальной математики. Новосибирск: ИМ СО РАН. 2001. T. 4, № 1(7). С. 83-101.

[4] КогАй В.В. Применение продолжения по параметру для численного исследования периодических решений автономных систем обыкновенных дифференциальных уравнений // Вест. НГУ. Новосибирск: Изд-во НГУ. 2002. T. 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 // Intern. J. of Pure and Appl. Mathematics. Bulgaria. Sophia. 2004. Vol. 14, N 4. P. 467-498.

[6] Фадеев С.И. О решении системы трансцендентных уравнений с параметром методом Ньютона // Вычисл. системы. Новосибирск: ИМ СО РАН. 1985. Вып. 108. С. 78-93.

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

[8] Фадеев С.И., Покровская С.А., Березин А.Ю., Гайнова И.А. Пакет программ STEP для численного исследования систем нелинейных уравнений и автономных систем общего вида. Новосибирск: Изд-во НГУ, 1998. 188 с.

[9] Годунов С.К. Обыкновенные дифференциальные уравнения с постоянными коэффициентами. Т. 1: Краевые задачи. Новосибирск: Изд-во НГУ, 1994. 264 с.

[10] Фадеев С.И. Программа численного решения нелинейных краевых задач для систем обыкновенных дифференциальных уравнений с параметром // Вычисл. методы линейной алгебры. Новосибирск: Наука. Сиб. отд-ние, 1990. С. 104-198.

[11] Бивиков Ю.Н. Курс обыкновенных дифференциальных уравнений. М.: Высш. шк., 1991. 304 с.

[12] Понтрягин Л.С. Обыкновенные дифференциальные уравнения. М.: ГИФМЛ, 1961. 312 с.

[13] Годунов С.К. Современные аспекты линейной алгебры. Новосибирск: Науч. книга, 1997. 390 с.

[14] Коллатц Л. Задачи на собственные значения. М.: Наука, 1968. 504 с.

[15] Годунов С.К., Антонов А.Г., Кирилюк О.П., Костин В.И. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. Новосибирск: Наука. Сиб. отд-ние, 1988. 456 с.

[16] Семендяев К.А. О нахождении собственных значений и инвариантных многообразий матриц посредством итераций // ПММ. 1943. Т. 7, № 3. С. 193-222.

[17] Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. М.: ГИФМЛ, 1960. 656 с.

[18] Крылов В.И., Бовков В.В., МОНАСТЫРНЫЙ П.И. Вычислительные методы высшей математики. Минск: Высш. шк, 1972. 584 с.

Поступила в редакцию 19 ноября 2004 г., в переработанном виде — 24 декабря 2004 г.