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

Вычислительные технологии

Том 11, № 5, 2006

ИССЛЕДОВАНИЕ ЦИКЛИЧЕСКИХ КОНТУРОВ ГЕННЫХ СЕТЕЙ С ОТРИЦАТЕЛЬНЫМ ТИПОМ

РЕГУЛИРОВАНИЯ

М.А. Клишевич, В. В. КогАй Институт математики СО РАН, Новосибирск, Россия e-mail: mansi@ngs.ru, kogai@math.nsc.ru

The core of gene networks (GN) is regulatory circuits — genes and proteins with mutually regulated expression. Thus, detection of possible operation modes of regulatory circuits is an important problem of systems biology. The work presents an efficient method for computing symmetric periodic solutions for cyclic regulatory circuits for M(n, 2) model.

Введение

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

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

dx1 а

x1,

¿Ь 1 + Хп ¿х2 а

= Т+Х1 -Х2' (1)

¿Хп аа

¿Ь 1 +

которую будем называть моделью М(п, 2), описывающей наиболее простые по строению циклические контуры молекулярно-генетических систем [4]. Здесь п — количество генетических элементов в системе. Каждый генетический элемент х^, г = 1,... ,п, кодирует

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

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

В работе [3] представлены следующие выводы.

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

2. В системе (1) всегда имеется одна симметричная точка покоя, которая при достаточно малых а, 7 является единственной и устойчивой стационарной точкой.

3. При достаточно больших а, 7 симметричная точка теряет устойчивость.

4. При достаточно больших а, 7, если п — четное, в системе возникают две частично симметричные устойчивые точки покоя.

5. При достаточно больших а, 7, если п — нечетное, в системе (1) существует один симметричный устойчивый предельный цикл.

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

При изучении периодических решений в модели гипотетических генных сетей М(п, 2) найдены симметричные решения, все компоненты которых идентичны по форме и отличаются друг от друга только сдвигом по фазе:

Ргт т = — Т , п

где р — некоторое натуральное число с учетом ограничения р < п, так как решение для р = кп+г соответствует решению для р = г, г < п, причем значение тп должно быть кратно периоду Т. Таким образом, любое симметричное периодическое решение со сдвигом т удовлетворяет уравнениям

х^^) = х( — т), г = 2, ...,п, хп(£) = х\(1 — т), которые при подстановке в (1) дают систему из п идентичных уравнений

dxl(t) а . .

— Х2 (¿), Хп^)-

dt dx2(t) 1 + х1 ^ — а т)

dt 1 + х2 ^ — т)

dxn(t) а

dt 1 + хП(г — т)

Таким образом, достаточно исследовать уравнение

dx(t) а

dt =1 + х^ (£ — т)

х(т) (2)

с параметром

- Р™ . п - 1

Т = - т, Р = 1,...,-

п п

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

На рис. 1 изображено симметричное периодическое решение для модели М(16, 2). Вторая компонента периодического решения сдвинута влево относительно первой на 7/16

7

периода (всего кривых 16), поэтому запаздывание равно Г = — Т. На рис. 2 изображено

7

16

решение для уравнения (2) с запаздыванием Г = — Т. Первая компонента на рис. 1 совпадает с решением на рис. 2, т. е. для построения полного решения достаточно размножить кривую на рис. 2 со сдвигом Г влево.

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

(х1 (г) (х2 (г)

а

1 + (г - г)

а

гда — Ж2(г)-

— х1(г), Г = РТ, п

(3)

Рассмотрим уравнение (2), которое само является моделью авторепрессилятора одиночного генетического элемента, способного подавлять собственную экспрессию. Для нахождения циклов уравнения (2) запишем краевую задачу, которая будет давать как устойчивые, так и неустойчивые периодические решения [5, 6, 7]. Приведем задачу к единичному отрезку, сделав замену

Г = гт, г = тТ, х(гт) = и(г),

и получим уравнение с нормированным по периоду Т временем

(и(г)

(г

=Т

а

1 + и^ (г — т)

— и(г)

(4)

Рис. 1. Симметричное периодическое решение модели М(16, 2) при а = 25, 7 = 10, т = 7/16. Время £ нормировано по периоду Т.

Рис. 2. Первая компонента симметричного периодического решения М(16, 2), а = 25, 7 = 10.

Добавим условие периодичности для (4)

и(0) - и(1) = 0 (5)

и условие трансверсальности

а

и(0) - 1+^4—7) =(6)

В данном случае условие (6) определяет, что при £ = 0 функция и(£) имеет нулевую производную.

Таким образом, приходим к краевой задаче (4)-(6) с параметрами а, 7, т для неизвестной функции и(£) и периода Т.

1. Дискретная модель краевой задачи для уравнения с запаздывающим аргументом

Будем искать решение в виде кусочного эрмитова сплайна [8]. Для этого введем сетку

0 = ¿1 < ¿2 < ... < ^ = 1

и сведем задачу к системе нелинейных уравнений, проинтегрировав обе части равенства на интервале [¿¿,¿¿+1]:

¿¿+1 ¿¿+1

Иг+1 — и» = аТ J f (и(£ — т— Т I г = 1,..., N — 1,

t¿ t¿

введем обозначение

¿¿+1 ¿¿+1

^ = Иг+1 — и» + Т / — аТ / f(и(£ — т= 0, г =1,...,Ж — 1, (7)

где T — период, т — запаздывание, f (x) =

1 + xY

Заменим интегралы в равенстве (7) квадратурными формулами. Функцию u(t) будем аппроксимировать кубической параболой Эрмита S(t) класса C1:

t = ti + his, 0 ^ s ^ 1, hi = ti+i - ti, i =1,...,N - 1,

U = u(ti) = Si(ti), ui+i = u(ti+i) = Si(ti+i), ui = u'(ti) = Si(ti), ui+1 = u'(ti+i) = Si(ti+i). Кубическая парабола и ее производная имеют вид

Si(t) = (1 - s)ui + sUi+i + s(1 - s)[(1 - s)a + sbi], (8)

hiSi(t) = Ui+i - Ui + (1 - 2s)[(1 - s)a + sbi] + s(1 - s)(b - a), где параметры a и bi определяются из условий

hiui = Ui+i - Ui + ai, hiui+i = Ui+i - Ui - bi,

t

t

т. е.

а = —г - —г+1 + Лг—г, 6 = —г+1 - Иг - а + 6 = Лг(—г -

Отсюда следуют формулы

4+1

5= Л(—г + мг+1) + Л(аг + 6) = У(иг + —+ 12 (—г - —г+1^ (9)

—г = 5(¿г + Л) = ^(мг + «¿+1) + ^г + 6) = 1 (—г + —г+1) + ^Т (—г - И^). (10)

2 2 о 2 о

Согласно формуле (9)

■¿+1

Лг Л2

— СО-: = у(—г + Мг+1) + 12(иг - иг+1)-

(11)

Второй интеграл в (7) приближенно заменим квадратурной формулой Симпсона:

¿+1

Лг

/ (—^ - т))-& = -6 / («(¿г - т)) + 4/ ^ 2

¿г + ¿г+1

- Т + / (—(¿г+1 - Т)) • (12)

Значения функции —(¿) в неузловых точках можно вычислять либо по формуле кубической параболы, построенной по точкам —г-1, —г, —г+1, —г+2, либо по формуле параболы Эрмита (8) с использованием значений —г, —г+1, и г, Мы используем второй вариант,

потому что он более устойчивый. Так как значения гг = -— (¿г) не представляется возможным определить явно по формуле (4), мы их находим из системы:

^N+1 = г1 - Т +2 = г2 - Т

а

1 + —^ (¿1 - т)

а

1 + —•у (¿2 - т)

-1 = гМ-1 - Т

- —1 —2

а

1 + —^ (¿м-1 - т)

: 0, : 0,

—N-1

(13)

к = г1 - гм = 0,

где значение —(¿^ - т) определяется по формуле (8), значение индекса г такое, что:

г Ч ¿Й - т - ¿г ¿й - т е [¿г, ¿г+1), 5 = —-— для случая ¿й - т ^ 0;

¿г+1 - ¿г г ч 1 + ¿Й - т - ¿г

1 + ¿Й - т е [¿г,¿г+1), 5 =----- для случая ¿Й - т < 0.

¿г+1 - ¿г

Выпишем уравнения, которые соответствуют краевым условиям (5), (6):

= —1 - — М = 0, +1 = г1 = 0.

г

г

г

0

Представим уравнения (7) с учетом интегралов (11), (12) и уравнений (13), (14) как систему из (2Ж+1)-го нелинейного уравнения, содержащую скалярные параметры а, 7, т:

= «2 - «1 + Т

—аТ — 6

у («1 + «2) + ^!(гл - ^

+

1+ и № - т) ■ ! + т

^2 = «3 - «2 + Т

^2 -аТ— 6

2

у («2 + ^ (г2 - г3)

+

+

1 + ^ (¿2 - т)

1+ в - т > ' 1 + т

+

1+ п1 (¿3 - т)

^м-1 = «м - «м-1 + Т

км-1

к2Л

(«м-1 + «м) + (гм-1 - гм)

12

аТ

км-1

6

1

+

4

1+ ^ (¿м-1 - т) ' / ¿м-1 + ¿м

1 + «ч -1--т

= - «м = 0 +1 = г1 - Т

2

+

1 + ^ (¿м - т)

а

+2 = г2 - Т

1 + ^ (¿1 - т)

а

1 + ^ (¿2 - т)

«1

- «2

-1 = -1 - Т

= г1 - гм = 0, , ^2м+1 = г1 = 0,

а

1 + ^ (¿м-1 - т)

= 0, = 0,

- «м-1

где значения в неузловых точках «(¿) определяются по формуле кубической параболы (8), построенной на узлах ¿¿,¿¿+1, £ € [¿¿,¿¿+1).

Перепишем последнюю систему в векторном виде:

^(и, г, Т, т, а, 7) = 0.

(15)

4

1

1

0

4

1

1

0

2

1

0

0

Таким образом, формально проблема численного исследования решения краевой задачи (4)-(6) сводится к численному исследованию системы (15) с параметрами а, 7 и т, определяющими вид правых частей системы относительно векторов г и неизвестного периода Т. Решение находится стандартными методами с использованием процедуры продолжения по текущему параметру, учитывающей множественность решений [5, 6, 9].

2. Адаптация сетки

При продолжении по параметру решения системы нелинейных уравнений, аппроксимирующей дифференциальную краевую задачу, возникает, вообще говоря, проблема адаптации первоначально заданной сетки по Ь к изменению параметра т. Для решения этой задачи мы используем аналитическое решение, которое задача (4)-(6) имеет при 7 = те. Рассмотрим краевую задачу (4)-(6) при 7 = те:

и' (Ь) =

а — и(Ь), если и(Ь — т) < 1,

а

— — и(Ь), если и(Ь — т) = 1, 2

(16)

—и(Ь),

если

и(Ь — т) > 1,

и(0) — и(Т) = 0,

и(0) = шт и(Ь).

Щ0,Т ]

(17)

(18)

Для задачи (16)-(18) можно выписать точное решение, для этого проинтегрируем правые части (16):

и(Ь) =

_ \ и = а — С1е если и(Ьу — т) = 1, Ь £ (Ьу, Ьу+1) ^ и(Ь — т) < 1,

и2 = С2е *, если и(Ьу — т) = 1, Ь £ (Ьу,Ьу+1) ^ и(Ь — т) > 1.

(19)

Любое решение (16) представляет собой непрерывную кривую, составленную из участков, показанных на рис. 3. Для вычисления характеристик Ь1, Ь2, Ь3, Т и констант С1, С2 в (19) будем использовать условие периодичности, непрерывность и(Ь) в точке Ь2, условие периодичности, равенство единице в точках Ь1, Ь3 и равенство запаздываний величине т:

и1 (0) = и2 (Т), и1 (Ь2)= и2(Ь2), и1(Ь1) = 1, и2(Ь3) = 1, Ь2 — Ь1 = т, Т — Ьз = т.

(20)

Рис. 3. Аналитическое решение п(Ь) задачи Рис. 4. Решение краевой задачи (4)-(6) при (16)-(18). а = 25, 7 = 10, т = 7/16 и аналитическое

решение (22) при том же значении т.

Пусть т и а — некоторые произвольные параметры (с учетом ограничения а > 1), они же есть фиксированные неотрицательные числа. Решая (20), находим

С = (а - 1)е41, С2 = (а - (а - 1)е-тУ2,

(21)

а - е-т ает - 1 (ает - а + 1)(ает - 1)

¿1 = 1п-—, ¿2 = 1п-—, ¿з = Т - т, Т = 1п---.

а - 1 а - 1 а - 1

Прямой проверкой убеждаемся, что для любых а > 1, 7 > 0 выполняются естественные соотношения 0 < ¿1 < ¿2 < ¿3 <Т.

Таким образом, мы можем в явном виде выписать решение задачи (16)—(18):

«(*)= а- (а- х-л 1 ;е (22)

а - (а - 1)е*1-4, г е [0, ¿2], [а - (а - 1)е-т]в*2-4, Ь е [¿2,Т].

На рис. 4 приведены решения краевых задач (4)-(6) и (16)-(18). Так как при достаточно больших значениях а и 7 они близки, аналитическое решение задачи (16)—(18) может быть использовано для построения неравномерной сетки.

Рассмотрим способ вариации узлов сетки, позволяющий повысить точность аппроксимации численного решения краевой задачи [10], основанный на методе выбора узлов интерполяции функций с заданной погрешностью интерполяции [11]. Согласно этому методу функция и(£) приближается эрмитовым кубическим сплайном [11] класса С1 на каждом из интервалов [0, ¿2] и [¿2,1]. При этом на каждом из элементарных отрезков , ^+1] сплайн представлен кубической параболой, построенной по значениям функции и ее производной в узловых точках. Узлы сетки

0 = & < ... < 0 = ¿2 < 0+1 < ••• = 1, Ы = £+1 - & (23)

определяются из равенств

Ы

где е — задаваемая погрешность интерполяции.

£ <*>

1/4

= (384е)1/4, (24)

Рис. 5. Решение для а = 25, 7 = 10, т = 2.

При построении сетки (23) используется эквивалентное условие

1

т— 1

с1/4 £ нг

г=1

£ <*>

1/4

(Ж - 1)е1/4 « С1/4

¿4и . .

1/4

<И

или

1

Ж- 1

(25)

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

Необходимо отметить, что при малых т решения задач (4)-(6) и (16)—(18) начинают различаться. Поэтому при т < 3 используется равномерная сетка, которая дает хорошие результаты ввиду небольших значений градиентов решения (рис. 5).

3. Построение диаграммы решений методом продолжения по параметру

На основе продолжения решения по параметру т строится диаграмма решений в виде зависимости Г(т) = Т(т)/т, которая описывает отношение периода Т к значению запаздывания т. Диаграмма решений модели М(п, 2) с фиксированными параметрами а и 7 позволяет определять количество симметричных периодических решений для произвольного п и численно находить их.

Рассмотрим систему (1) с фиксированными п, а и 7. Пусть в ней имеется симметричный цикл. Это означает, что существует такое натуральное р, 1 ^ р ^ п, что п-я

о

1

ь

о

Рис. 6. Диаграмма решений задачи (4)—(6)при а = 25, 7 = 10.

переменная сдвинута вправо относительно первой переменной по фазе на величину, равную т = pT/n. Следовательно, в уравнении (2) при том же значении параметра а и т = pT/n существует циклическая траектория, идентичная по форме первой переменной рассматриваемого цикла системы (1). Обратное также верно. Поэтому ответ на вопрос о существовании или отсутствии для некоторого фиксированного 1 ^ p < n цикла в системе (1) зависит от того, пересекает прямая n/p график функции F(т) или нет (см. рис. 4).

Функция F(т) позволяет подсчитать все симметричные периодические решения системы (1) для данной зависимости при любом n, не решая системы. Если учесть, что F(т) > 2, то для того чтобы подсчитать количество решений, достаточно построить [n — 1/2] прямых y = n/p, p = 1,..., [n — 1/2], на графике F(т) и для каждой подсчитать количество точек пересечения. Например, для модели M(5, 2) имеем три решения (рис. 6).

Продолжение влево по параметру возможно только до некоторого значения т*, зависящего от а и y , которое может быть вычислено по формулам [3]

Для т < т* задача (4)-(6) не имеет решений.

4. Численные результаты

На рис. 6 для модели М(п, 2) приведена диаграмма решений ^(¿) при а = 25, 7 = 10. Используя данную зависимость, можно определять количество всех симметричных периодических решений для любого п. Однако вопрос о единственности такого рода зависимости требует дальнейшего изучения.

На примере модели М(5, 2) подсчитаем количество симметричных периодических решений и найдем их качественно. Для этого построим прямые у = 5 и у = 5/2 на графике (рис. 6). Из рисунка видно, что количество точек пересечения равно трем, т.е. модель М(5, 2) имеет три симметричных периодических решения. Значения т, которые соответствуют точкам пересечения, определяют параметры т для задачи (4)-(6), отвечающие периодическим симметричным решениям модели М(5, 2), что позволяет строить множественность решений для модели М(5, 2). Определим значения т, соответствующие точкам пересечения:

т1 = 6.3313, т2 = 0.7503, т3 = 0.4181.

На рис. 7 приведены в виде графиков результаты численного исследования периодических решений для соответствующих параметров т. Графики представляют зависимость «(¿) , где время ¿ нормировано по периоду Т.

Таким образом, для параметров а = 25, 7 =10 мы можем построить таблицу решений модели М(п, 2) для различных значений п.

Количество симметричных периодических решений при а = 25

Значение

Количество генетических элементов в системе

Y 2 3 4 5 6 7 8 9 10 11 12

10 0 1 2 3 1 2 3 4 4 3 4

15 0 1 2 3 3 2 3 4 4 5 6

а

Рис. 7. Решение задачи (4)-(6) при а = 25, 7 = 10: а — т = 6.3313; б — т = 0.7503; в — т = 0, 4181.

Цф

-

5- \

4- \

2-

-1— 10

Рис. 8. Диаграмма решений задачи (4)-(6) при а = 25, 7 = 15.

Рис. 9. Несимметричное решение модели М(12, 2) а = 25, 7 = 10.

Рис. 10. Несимметричное неустойчивое периодическое решение модели М(12, 2) при а = 25, 7 =10 в фазовой плоскости (^2,^1).

Если построить зависимость Г(т) при а = 25, 7 =15 (рис. 8), то вид кривой сохранится, но увеличится значение Г(т) в точке поворота с ^ 5.5 до ^ 6.5. В таблице отображено количество решений в зависимости от п для а = 25, 7 =15.

Можно заметить, что, например, для модели (6, 2) при 7 =10 количество решений на два меньше, чем при 7 = 15. Отметим, что для определения числа решений модели М(п, 2) достаточно знать параметры в2 (рис. 8).

При исследовании модели М(12, 2) численно найдено периодическое решение, в котором все компоненты разбиваются на две группы, в каждой из которой компоненты отличаются только на постоянный сдвиг по фазе, равный 5/6 (рис. 9). Такие решения найдены также при исследовании моделей М(14, 2) и М(16, 2). Краевая задача для уравнения (3) решается аналогичным методом. Это позволяет сделать вывод о том, что в моделях М(п, 2) могут быть и несимметричные периодические решения. Однако найденное решение является неустойчивым, что проверяется интегрированием (рис. 10).

Заключение

Изучение свойств моделей М(п, 2) — актуальная задача теории генных сетей. В частности, данные модели описывают при п = 2 молекулярный триггер, а при п = 3 — репрессиля-тор — конструкции, реализованные генно-инженерными методами [12, 13]. Исследованная в работе модель представляет значительный интерес при решении задач конструирования более сложных молекулярно-генетических систем [14]. Поэтому важно обнаружить в простейших регуляторных контурах все типы динамического поведения.

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

Численные исследования показывают, что при четном n в моделях M(n, 2) могут существовать частично-симметричные замкнутые траектории, все компоненты которых делятся на две группы. В каждой группе траектории идентично сдвинуты одна относительно другой на некоторую фазу. Все такие периодические решения удовлетворяют системе уравнений с запаздыванием (3), которая представляет модель молекулярного триггера, в котором учитываются времена синтеза белков-регуляторов. Следовательно, в молекулярном триггере можно ожидать наличие асимметричного цикла. Исследование данного цикла показывает, что он полуустойчив. Поэтому переходный режим из одной точки покоя в другую (переключение молекулярного триггера) может протекать по колебательной траектории, наматывающейся на данный цикл сколь угодно долго. Визуально данная траектория выглядит как периодическая траектория. Но так как цикл не является устойчивым, за счет флуктуаций система (3) рано или поздно перейдет в устойчивый режим (рис. 10). В дальнейшем мы намерены изучить параметры осцилляции более сложных моделей молекулярного триггера с целью оценки возможности его существования в природных системах.

Выражаем благодарность С.И. Фадееву и В.А. Лихошваю за обсуждение работы и ценные замечания.

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

[1] КолчАнов Н.А., АнАнько Е.А., КолпАков Ф.А. и др. Генные сети // Молекулярная биология. 2000. T. 34, № 4. C. 449-460.

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

[3] ЛихошвАй В. А., Матушкин Ю. Г., Фадеев С. И. Задачи теории функционирования генных сетей // Там же. С. 64-80.

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

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

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

[7] КогАй В.В., Фадеев С.И., ЛихошвАй В.А. О численном исследовании автоколебаний в гипотетических генных сетях // Вычисл. технологии. 2005. Т. 10, № 3. С. 56-71.

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

[9] 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. Math. 2004. Vol. 14, N 4. P. 467-498.

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

[11] Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980. 352 с.

[12] Elowitz M. B., Leibler S. A synthetic oscillatory network of transcriptional regulators // Nature. 2000. Vol. 403(6767). P. 335-338.

[13] Gardner T.S., Cantor C.R. et al. Construction of a genetic toggle switch in Escherichia coli // Nature. 2000. Vol. 403(6767). P. 339-342.

[14] Sprinzak D., Elowitz M.B. Reconstruction of genetic circuits // Nature. 2005. Vol. 438. P. 443-448.

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