Задачи теории функционирования генных сетей

ЗАДАЧИ ТЕОРИИ ФУНКЦИОНИРОВАНИЯ ГЕННЫХ СЕТЕЙ
В.А.Лихошвай* , Ю.Г.Матушкин, С.И.Фадеев¹
Институт цитологии и генетики Сибирского отделения Российской академии наук, Новосибирск, 630090.

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

контактный автор likho@bionet.nsc.ru

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

В работе приводятся первые практические результаты по исследованию ГГС. Приводятся примеры ГГС, имеющих множество стационарных решений, предельные циклы и аттракторы сложной природы. В соответствии с типами регуляторных связей в рассмотрение вводятся 4 класса ГГС, а также определяются в них канонические представители. Для канонических симметричных ГГС, формулируется (n,k)-критерий, позволяющий предсказывать в определенной области значений параметров предельные свойства решений соответствующих моделей, т.е. предсказывать число и характер имеющихся у них устойчивых предельных состояний. В случае ГГС первого и второго классов формулируется гипотеза об однозначном соответствии устойчивых состояний и 1-баз структурных орграфов этих сетей, а в случае ГГС третьего и четвертого классов - гипотеза об однозначном соответствии устойчивых состояний и 1m-баз структурных орграфов.

Вместе с этим авторы хотели бы привлечь внимание к проблемам качественного и численного анализа рассматриваемого круга проблем.

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

Введение

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

а) поддерживают в организме устойчивые состояния, характеризующиеся постоянством концентраций веществ (гомеостаз, стационарные состояния);

б) обеспечивают периодические незатухающие колебания концентраций определенных групп веществ (осцилляции, циклы);

в) контролируют необратимые процессы: развитие, рост, дифференцировка, апоптоз.

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

1. МРНК и белки, кодируемые этими генами;

2. пути передачи сигналов от клеточной мембраны к ядру клетки, обеспечивающие активацию генов в ответ на внешние регуляторные стимулы;

3. внешние сигналы, гормоны, метаболиты, которые переносят физиологические регуляторные воздействия, обеспечивая изменение состояний генной сети.

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

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

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

Здесь n – количество динамических переменных модели генной сети, X=(x₁,…,x_n)– вектор динамических переменных, K=(k₁,…,k_m) – вектор параметров, F_i(X_i,K) и G_i(X,K) – рациональные функции, принимающие неотрицательные значения при всех неотрицательных значениях переменных и параметров, такие что уравнение имеет неотрицательное решение относительно x_i при любых заданных значениях параметров K=(k₁,…,k_m) и X_i_′=(x₁,…, x_i_′-1, x_i_′+1,…,x_n). Таким образом, изучение динамики функционирования генных сетей можно рассматривать как задачу изучения свойств систем вида (1).

Как уже отмечалось, ключевую роль в функционировании генных сетей играют отрицательные и положительные регуляторные связи. Для их изучения введем гипотетические генные сети (ГГС), которые состоят из элементов двух типов: элементов матричного синтеза (коротко, m-элементов) и регуляторных механизмов. Под m-элементом мы будем понимать упорядоченную пару, из которых первый элемент будем называть порождающим, а второй – порождаемым или продуктом. В данных объектах перенос информации и вещества осуществляется в направлении от порождающего элемента к продукту. Физическими прообразами m-элементов являются пары (ген, белок), (ген, РНК), (мРНК, белок) и т.д. В них синтез продукта осуществляется под контролем порождающего элемента, при этом сам порождающий элемент не расходуется в процессе синтеза продукта. Для простоты мы не будем вводить в понятие m-элементов пары, в которых перенос информации осуществляется в обратном направлении: РНКДНК, Белок РНК, хотя формальных трудностей это не составляет. Более того, в дальнейшем мы ограничимся рассмотрением только элементов, прообразами которых являются пары (ген, белок). Такие m-элементы назовем генетическими (ГmЭ). Делается это единственно с целью облегчения дальнейшего изложения. Все излагаемое ниже может быть легко обобщено и на m-элементы другой природы.

В генетических элементах продуктами являются белки, а порождающими - гены. Обозначать ГmЭ будем в виде упорядоченной пары (g, p), где g обозначает ген, а p - белок. Два генетических элемента будем считать идентичными, если они имеют одинаковые наименования генов и белков.

Под функционированием ГmЭ будем понимать синтез белка, кодируемого этим элементом. Под активностью функционирования ГmЭ будем понимать скорость синтеза его белка.

Элементами второго типа являются регуляторные механизмы (РМ). Под РМ мы будем понимать процесс, посредством которого некоторая совокупность генетических элементов осуществляет регуляцию активности некоторого фиксированного генетического элемента. Вхождение конкретного ГmЭ (g₁,p₁) в РМ другого конкретного ГmЭ (g₂,p₂) будем называть регуляторной связью (РС). В этом случае ГmЭ (g₁,p₁) будем называть регулятором (g₂,p₂). Если (g₁,p₁) идентичен (g₂,p₂), то регуляторную связь будем называть авторегуляторной.

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

3. Системы, описывающие функционирование ГГС

3.1. Описание конститутивного синтеза белка в отсутствие регуляции.
Рассмотрим сначала ГmЭ (g,p), который имеет некоторый конститутивный уровень функционирования и не подвержен регуляции другими генетическими элементами. Его функционирование выражается в синтезе белка p. А эффективность функционирования – в скорости синтеза белка p. Хорошо известно, что синтез белка в живых клетках слагается из значительного количества последовательно протекающих стадий. Общее количество промежуточных стадий может быть достаточно большим. Для разных генов оно является различным, так как они кодируют РНК и белки различной длины. Для гена, с которого синтезируется РНК длиной 1000 нуклеотидов и белок длиной 300 аминокислот минимальное количество промежуточных стадий, только этапов элонгации транскрипции и трансляции, будет составлять 1900, так что реальное количество их значительно больше.

Будем считать, что количество промежуточных стадий равно n_g. Обозначим каждую промежуточную стадию через x_i_,,. Конститутивный синтез белка p (переменная x) будем описывать системой дифференциальных уравнений:

(2)

где  - константа скорости конститутивного синтеза, – константы скорости перехода из i-го промежуточного состояния в (i+1)-ое промежуточное состояние,  - константа деградации белка p. Как видим, промежуточные стадии (переменные x_i) синтеза белка в отсутствие регуляции мы описываем линейными уравнениями. В реальности (в живых системах) промежуточные стадии протекают по более сложным нелинейным механизмам. Однако наш опыт моделирования показывает, что при некоторых естественных ограничениях, которые обычно выполняются для биохимических систем, они могут быть хорошо приближены линейными системами (2). (Более подробно об этом говорится в разделе 4.2.)

3.2. Описание регуляции генетического элемента
Будем полагать, что в ГГС регуляция активности генетического элемента направлена на изменение (в сторону уменьшения или увеличения) активности функционирования его порождающего элемента. Будем полагать, что генетические элементы-регуляторы осуществляют функцию регулирования посредством кодируемых ими белков. На формальном уровне наличие регуляции выражается в том, что в системе (2) параметр  становится функцией, зависящей от концентраций белков, кодируемых генетическими элементами-регуляторами данного генетического элемента.

Описывать регуляцию активности ГmЭ будем следующим образом. Считаем, что порождающий элемент g может пребывать в нескольких альтернативных состояниях. Обозначим эти состояния через u_i, i=1,…,n_u. Считаем, что нам известна схема перехода состояний друг в друга. Эту схему будем задавать ориентированным графом, в котором вершины представляют альтернативные состояния, а дуги – переходы между ними. Каждый переход осуществляется с некоторой скоростью v_i_,_j, которая зависит от концентраций ряда белков-регуляторов, что выражается в конкретном виде зависимости v_i_,_j от концентраций белков-регуляторов. По смыслу v_i_,_j  0 при любом наборе неотрицательных переменных. Припишем каждой существующей дуге (u_i,u_j) в качестве веса выражение v_i_,_j, а каждой i-й вершине припишем переменную u_i, которая будет иметь смысл вероятности пребывания регуляторного участка в i-м состоянии. Составим следующую алгебраическую систему из n_u уравнений на неизвестные u_i:

. (3)
Здесь

- множество вершин, из которых дуги входят в i-ю вершину,

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

. (4)
Исходя из биохимической природы процессов, протекающих в живых клетках, для задания v_i_,_j мы ограничиться использованием полиномиальных выражений от концентраций продуктов генетических элементов-регуляторов:

, (5)

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

Система (3) с уравнением баланса (4) является линейной относительно неизвестных. Она всегда имеет единственное решение , u_i0, которое при необходимости может быть выписано в аналитическом виде.

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

(6)

Выражение A будем называть моделью регуляции генетического элемента. Заменим в системе (2)  на A. В результате мы приходим к системе, которая описывает функционирование регулируемого генетического элемента.

4. Построение моделей функционирования ГГС

4.1. Описание функционирования ГГС автономными системами.
Возьмем теперь некоторое количество ГmЭ (g_j,p_i), j=1,…,m, i=1,…,n. Очевидно, что по смыслу nm. Так как генетический элемент полностью определяется порождающим элементом g_j, то договоримся, что ГmЭ имеет тот же номер j. Считаем, что для каждого g_j задан механизм A_j и совокупность белков, участвующих в регуляции активности функционирования генетического элемента (g_j,p_i) является подмножеством всех белков, синтезируемых рассматриваемыми генетическими элементами. Чтобы записать систему дифференциальных уравнений, описывающую функционирование ГГС, состоящей из выбранной совокупности генетических элементов и заданных регуляторных механизмов, для каждого генетического элемента построим системы (2) с учетом соответствующего регуляторного механизма A_j.

Теперь договоримся, что если два генетических элемента кодируют один белок, то последние переменные соответствующих систем (переменные xв системе (2)) являются идентичными. Для всех таких переменных составим балансные уравнения вида:

“скорость наработки белка x” = ”сумма скоростей наработки белка x по всем генетическим элементам, кодирующих белок x” – “скорость распада белка x”.
К этим уравнениям добавим все уравнения из систем регуляции всех генетических элементов, исключая из них все последние уравнения.

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

, (7)
которую будем называть моделью ГГС. Здесь n – количество различных белков, m – количество различных порождающих элементов, n_j_,_i – количество промежуточных стадий в синтезе i-го белка j-м порождающим элементом,

- константа скорости перехода s-й промежуточной стадии в (s+1)-ю промежуточную стадию при синтезе i-го белка j-м порождающим элементом, _i – константа распада i-го белка, P_j – множество номеров белков, которые кодируются j-м порождающим элементом, R_j – множество номеров белков, которые регулируют активность функционирования j-го порождающего элемента, G_i - множество номеров порождающих элементов, которые кодируют i-й белок. Класс этих систем будем обозначать через МГГС.
4.2. Описание функционирования ГГС системами дифференциальных уравнений с запаздывающими аргументами
Как упоминалось выше, в реальных генных сетях синтез белков слагается из большого количества промежуточных стадий, суммарное время протекания которых ограничено, а каждая промежуточная стадия протекает достаточно быстро. Для описания промежуточных стадий мы используем линейные уравнения вида (2), а для описания механизмов регулирования – нелинейные выражения (6). Последние, с учетом (5), представляют рациональные выражения от концентраций соответствующих белков-регуляторов. Также легко видеть, что значительное количество уравнений (7) являются дифференциальными уравнениями с простейшей по виду линейной правой частью. Благодаря их наличию системы из МГГС, как правило, имеют очень большую размерность. Понятно, что большая размерность существенно затрудняет проведение анализа таких систем. Возникает вопрос, а нельзя ли как-нибудь понизить b[ размерность? Оказывается, стандартный вид линейных уравнений, входящих в системы (7), действительно позволяет сделать важный шаг по пути упрощения их и перейти к новому классу моделей, в которых число переменных в точности равно числу различных белков.

Наряду с системой (7) из МГГС будем рассматривать в качестве модели ее функционирования систему уравнений с запаздывающими аргументами следующего вида:

(8)
Здесь все параметры и обозначения, кроме _j_,_i , имеют тот же смысл, что и для (7), а

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

Отметим, что системы (8) являются в некотором роде минимальными системами, которые удовлетворяют требованию адекватности описания функционирования ГГС. Действительно, нижняя граница сложности вида такой системы вытекает из необходимости учета в ней всех регуляторных механизмов. Поэтому в системах, описывающих функционирование ГГС, должны присутствовать динамические переменные, имеющие смысл концентраций белков-регуляторов и механизмы регулирования их скоростей синтеза. С другой стороны, в таких системах должны быть минимизированы описания второстепенных деталей функционирования ГГС. Системы (8) полностью удовлетворяют перечисленным условиям, так как содержат в качестве динамических переменных только концентрации белков-регуляторов и не содержат никаких иных деталей описания кроме механизмов и запаздываний. Но запаздывания являются важным элементом описания функционирования ГГС, поскольку изменение величины запаздывания может приводить к кардинальному изменению поведения ГГС в целом.

В качестве примера рассмотрим уравнение

(9)
которое описывает автоингибирование генетического элемента. Очевидно, при любых неотрицательных значениях  и , как при =0, так и при  > 0, уравнение (9) имеет только одно стационарное решение, которое находится из уравнения
. (10)
Однако, при других значениях уравнение (9) может иметь в качестве устойчивого режима функционирования предельный цикл. Такая ситуация, например, реализуется, если при =5 и =3 мы начнем постепенно увеличивать значение . Тогда, начиная с критического значения * =0.9643, стационарное решение теряет устойчивость, что приводит к возникновению автоколебания. Можно показать, что критические значения ^*для уравнения (9) вычисляются по формуле:
^* = , , .

Обозначим класс систем (8) через МГГСЗ. Между системами из МГГС и МГГСЗ существует тесная взаимосвязь. Она выражается в том, что системы из МГГСЗ являются предельными случаями систем из МГГС. Исходя из опыта численного анализа систем из МГГС и МГГСЗ мы высказываем следующие утверждение, требующее доказательства.

Рассмотрим систему (7) из МГГС, для которой в нулевой момент заданы начальные условия:
(11)
Обозначим через - концентрации i-х белков в момент времени t, найденные как решение задачи Коши (7), (11) при фиксированных n_j_,_i и . Здесь вектор (n_1,_i,…,n_m_,_i) вынесен в верхний индекс i-го белка, чтобы подчеркнуть, что решение найдено при фиксированных по всем генетическим элементам количеств промежуточных стадий синтеза белков.

Возьмем теперь систему (8) из МГГСЗ, которая соответствует системе (7) из МГГС. Пусть в ней запаздывания , рассмотренные в совокупности по всем генетическим элементам, имеют k различных значений. Обозначим их через t_r, r=1,…,k. Не ограничивая общности, мы можем считать, что t₁<…k. Обозначим через N_r множество генетических элементов (g_j,p_i), имеющих запаздывание t_r. По определению для любой пары (j,i) из (8) полагаем

, x_i(t-_j,i)=0 при t-_j,i<0, x_i(0)=0.
Через x₁(t),…, x_n(t) обозначим функции, которые получаются следующим образом. Если t₁>0, то полагаем x_i(t)=0 для всех 0t1, i=1,…,n. Затем переходим к последовательному решению задач (7), (12) на интервалах [t₁, t₂),…, [t_k_-1, t_k), [t_k, ). Если t₁=0, то сразу переходим к последовательному решению задачи (7), (12) на интервалах [t₁, t₂),…, [t_k_-1, t_k), [t_k, ). А именно, задав в момент времени t₁ начальное условие

,
решаем задачу (7),(12) на интервале [t₁,t₂), откуда находим

.
Затем переходим к решению задачи (7),(12) на интервале [t₂,t₃) с начальным условием

что позволяет определить

.
Данную процедуру повторяем для каждого следующего интервала, вплоть до последнего. В результате получаем искомые функции x_i(t), определенные на всей временной оси.

Сформулируем теперь утверждение. Предположим, что n_j_,_i , для всех j и i, таких что пара (g_j,p_i) является генетическим элементом, устремляются к бесконечности так, что выполняются следующие условия:

Тогда для любого t  0 имеем:

.

5. Отображение ГГС в виде диаграмм и ориентированных орграфов

5.1. Построение структурно-функциональных диаграмм ГГС.
Пусть имеется некоторая гипотетическая генная сеть. Важными ее характеристиками являются

1. количество генетических элементов, из которых она состоит;

2. количество различных белков, синтезируемых генетическими элементами (напомним, что некоторые генетические элементы могут кодировать идентичные белки);

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

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

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

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

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

Первое выражение имеет вид

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

-ый активный регулятор, h_i_,_j_,_k_,_u, _i_,_j_,_k_,_u – степени мультимеризации, _i_,_j_,_v_,_k ,_i_,_j_,_v_,_k – неотрицательные константы,

P_j_,_v_,_k – множество белков, входящих в активный регулятор.

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

Второе формальное выражение имеет вид

.
Оно применяется для описания независимого воздействия нескольких активных регуляторов на генетический элемент. Простейшим примером в данном случае могут служить промоторы, которые пересекаются с оператором репрессии и с последним могут взаимодействовать несколько различных репрессоров. Например, в системах регуляции активности промоторов P_R, P_L, Prm фага  -репрессор и антиреперессор могут независимо связываться с сайтом O_R1 (а также с O_R2, O_R3, O_L1, O_L2, O_L3) [4,5].

Третье формальное выражение

.
Оно учитывает совместное воздействие нескольких неконкурирующих независимых регуляторных событий. Пример такой регуляции мы снова находим у промоторов Prm, P_R, P_L, фага . Так, участки ДНК, обозначаемые O_R1, O_R2, O_R3 (O_L1, O_L2, O_L3, соответственно) составляют опереторные участки, входящие в систему регуляции промоторов Prm, P_R (P_L, соответсвенно) фага  и участки не пересекаются в том смысле, что с каждым из них -репрессор или антирепрессор могут связываться независимо от того, заняты другие сайты репрессорами или нет [4,5].

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

.
Дополнительно будем считать, что для j-го генетического элемента, кодирующего i-й белок, имеется запаздывание в его синтезе, вызванное многостадийностью процесса. Запаздывание будем обозначать через _j. Наличие запаздывания выражается в том, что в механизме регулирования i-го белка в момент времени t выражение A_i_,_jберется в момент времени t-_j

Введем теперь следующий перечень ограничений на конструируемые ГГС.

1. Все запаздывания приняты равными 0.

2. Рассматриваются только отрицательные регуляторные связи.

3. Все генетические элементы кодируют различные белки.

4. Генетический элемент имеет только один регуляторный сайт.

5. Активные регуляторы являются гомомультимерами белков одного вида.

6. Для каждого генетического элемента имеется только по одному активному регулятору. Активные регуляторы являются гетеромультимерами белков разного вида.

7. Все активные регуляторы действуют только на один сайт генетического элемента.

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

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

Класс 1. ГГС строятся с ограничениями 1, 2, 3, 4, 5, 7. Их функционирование описывается системами следующего вида:

. (13)
Класс 2. ГГС строятся с ограничениями 1, 2, 3, 5, 8. Соответствующие системы имеют следующий вид:

. (14)
Класс 3. ГГС строятся с ограничениями 1, 2, 4, 6 и описываются системами следующего вида:

. (15)
Класс 4. ГГС строятся с ограничениями 1, 2, 3, 4, 6. Они описываются системами следующего вида

. (16)
Численные расчеты показывают, что ГГС всех классов могут иметь самые различные режимы устойчивого функционирования – точки покоя, циклы, а также более сложные режимы функционирования. Для иллюстрации приведем несколько примеров.

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

.
Здесь

– концентрации продуктов синтеза,

 0 и  > 1 – параметры модели. Численный анализ выявил наличие двух параметрических областей. В первой области имеется единственное устойчивое стационарное решение, к которому стремятся все решения, независимо от начальных данных. Во второй области модель имеет два устойчивых стационарных решения. Таким образом, уже самые простые модели ГГС могут иметь нетривиальные режимы поведения.

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

.
Соответствующий численный анализ выявил существование параметрической области, в которой все решения задачи Коши при задании любых асимметричных начальных условий ( т.е условий, для которых не выполняются соотношения x_i=x_j для любых i,j=1,…,n) стремятся к единственному устойчивому предельному циклу( рис. 3г).

Подчеркнем, что описанные режимы функционирования рассмотренных ГГС найдены в генных сетях аналогичных конструкций, полученных искусственно [6,7]. Это является важным свидетельством адекватности предложенного математического аппарата.

На рис. 4а приводится пример структурного графа, в котором семь вершин и восемнадцать дуг. Построим на нем ГГС первого класса. Соответствующая система уравнений имеет вид :
,,,
,

,,,
В результате численного анализа найдена параметрическая область, >3, >3, в которой модель имеет одно устойчивое стационарное решение, три устойчивых предельных циклов и пять глобально неустойчивых предельных циклов. Каждый неустойчивый предельный цикл является устойчивым в некотором собственном гиперпространстве размерности меньше 7 (рис 4в.). При этом выполняются соотношения: для первого цикла  x₁(t) = x₅(t+) = x₃(t+ 2) = x₇(t+3) = x₂(t+4) = x₆(t+5), для второго цикла  x₁(t) = x₆(t+)=

x₃(t+2) = x₅(t+3) = x₂(t+4) = x₇(t+5), для третьего цикла  x₁(t)=x₇(t +)=x₃(t +2)=x₆(t +3)=x₂(t+4)= x₅(t+5).

Наконец, рассмотрим еще один пример. На рис.5а. приведен структурный граф, содержащий девять вершин и двадцать четыре дуги. Соответствующая для ГГС первого класса система уравнений имеет вид:

Численный анализ позволил выявить параметрическую область >3, >3, в которой имеется одно устойчивое стационарное решение, один устойчивый предельный цикл и два неустойчивых предельных цикла. При этом есть основания полагать, что устойчивые колебания при выбранных параметрах являются квазипериодическими.
7. Канонические ГГС и гипотезы об их поведении
Системы уравнений (13)-(16) описывают функционирование ГГС четырех классов. Как видно из структуры правых частей, данные системы содержат большое количество параметров. Это сильно затрудняет непосредственное проведение качественного анализа их поведения во всем параметрическом многообразии, в частности затрудняет выявление связи структуры гипотетической генной сети с ее динамическими свойствами. Поэтому естественно первоначальное исследование ГГС проводить, ограничившись минимально возможным набором параметров. С этой целью введем в рассмотрение класс канонических гипотетических генных сетей.

Канонической ГГС называется генная сеть, описываемая одной из систем (13)-(16) при следующих ограничениях на параметры: _i= _i,k= , _i = _i,k= 1, _i,k=, _i= 1. Таким образом, в канонических генных сетях оставлено всего два независимых параметра:  и . Тем не менее, этого вполне достаточно с биологической точки зрения, так как  имеет смысл конститутивной скорости синтеза белков, а,  имеет смысл степени нелинейности регуляции. Можно ожидать, что именно в канонических ГГС должно максимально проявляться влияние структурной организации на формирование возможных устойчивых режимов их функционирования. Это ожидание вполне подтверждается численными исследованиями, которые мы обобщаем ниже в виде гипотез-критериев.

При формулировании критериев будем пользоваться терминологией из теории графов, излагаемой в [8]. Пусть имеется ориентированный граф G(V,W), где V – множество вершин, а W – множество дуг. В силу важности приведем определение понятия 1-базы, изменив только обозначения (см.[8], стр.237):

1-базой называется минимальный набор U таких попарно несмежных вершин, что любая вершина орграфа G(V,W) или принадлежит U, или смежна из некоторой вершины U.

Введем понятие 1m-базы. Множество U вершин орграфа назовем 1m-базой, если и только если все вершины из множества V\U попарно несмежны, все вершины, имеющие нулевые полустепени заходов, лежат в U и любая вершина из U, имеющая ненулевую полустепень захода, смежна хотя бы из одной вершины, содержащейся в V\U.

7.1 Критерий о числе устойчивых состояний канонических ГГС классов 1-4.
Пусть дан структурный орграф G. Пусть M_i(G) обозначает модель канонической ГГС i-го класса, построенную на орграфе G. Тогда существуют такие ₀ и ₀, что при любых >₀ и >₀

1. для канонических ГГС первого и второго классов количество имеющихся у них устойчивых состояний равняется числу 1-баз орграфа G;

2. для канонических ГГС третьего и четвертого классов количество имеющихся у них устойчивых состояний равняется числу 1m - баз орграфа G.
Далее рассматриваются только первый класс.
3. если каждая вершина G принадлежит хотя бы одной 1-базе, то M₁(G) не имеет других устойчивых состояний;

4. если каноническая гипотетическая генная сеть 1-го класса, построенная на ограничении структурного орграфа G на совокупность вершин, которые не содержатся ни в одной его 1-базе, не имеет устойчивых колебательных режимов функционирования ни при каких значениях параметров  и , то и исходная каноническая ГГС S₁(G) не имеет устойчивых колебательных режимов функционирования.

5. если каноническая гипотетическая генная сеть 1-го класса, при достаточно больших значениях параметров  и , содержит устойчивые колебательные режимы функционирования, то каноническая гипотетическая генная сеть, построенная на ограничении структурного орграфа G на совокупность вершин, которые не содержатся ни в одной его 1-базе, также содержит при достаточно больших значениях параметров  и  устойчивые колебательные режимы функционирования.
Можно показать, что если орграф G(V,W) не имеет вершин с нулевой полустепенью захода, то некоторое подмножество его вершин U тогда и только тогда является 1-базой, когда V\U является 1m-базой. Для таких орграфов пункты 1 и 2 являются эквивалентными.

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

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

Пункт 4 критерия дает только достаточное условие отсутствия устойчивых колебательных режимов. Наличие устойчивых колебаний у канонической ГГС, построенной на ограниченном орграфе еще не гарантирует наличие таких режимов у исходной ГГС.

7.2 Описание глобальных режимов функционирования канонических ГГС, конструируемых на симметричных структурных графах.
Орграф G(V,W) назовем симметричным, если и только если все вершины имеют одинаковые полустепени исходов и существует такое расположение вершин орграфа на окружности, что все достижимые из u вершины, располагаются непосредственно после нее, считая по ходу часовой стрелки. В этом случае множество D_i номеров вершин смежных в вершину u_i в системах (13)-(16) записывается в следующем виде:

D_i={mod_n(j) j=i-1,…, i-k-1},

Так как симметричные орграфы полностью характеризуются числом вершин n и полустепенью исхода (k-1), то будем данные орграфы обозначать через G(n,k). На рис. 2б и 3б приводятся примеры простейших симметричных орграфов. Смысл использования в качестве обозначения не значения полустепени, а значения, увеличенного на единицу, становится ясным из формулируемого ниже (n,k)-критерия.

Обозначим через М_i(n,k) каноническую ГГС i-го класса, построенную на орграфе G(n,k). Имеет место следующий (n,k)-критерий существования устойчивых предельных режимов функционирования моделей M_i(n,k).

Пусть d - наибольший общий делитель n и k. Если d=k, то существуют такие ₀ и ₀, что при любых >₀ и >₀каждая из моделей М_i(n,k) имеет ровно k устойчивых стационарных решений, а все другие предельные режимы функционирования М_i(n,k) в указанной параметрической области неустойчивы.

Если dk, то существуют такие ₀ и ₀, что при любых >₀ и >₀модель М₁(n,k) имеет d устойчивых предельных циклов, При этом в указанной параметрической области другие устойчивые режимы функционирования гипотетической генной сети 1 класса невозможны.
8. О некоторых свойствах канонических ГГС первого класса.
Для понимания свойств гипотетических генных сетей важно не только знать о наличии или отсутствии у них устойчивых режимов функционирования, но также важно знать о наличии неустойчивых точек покоя и циклов, а также важно знание отличительных характеристик предельных режимов и умение находить их численно и аналитически. Приведем описание предельных свойств для наиболее изученного нами случая канонической ГГС первого класса.
8.1. Структура точек покоя M₁(n,k)
1. При любых  и  точки покоя удовлетворяют соотношениям x_i = x_i₊_jd, i=1,…,d, j=1,…,n/d-1.

2. Для любой точки покоя существуют, по крайней мере, два номера 1ijd, такие что x_i=x_j Таким образом, модель М₁(n,k) не имеет асимметричных точек покоя, т.е. таких, что для любых 1ijd выполняется неравенство x_ix_j.

3. В параметрической области >₀ и >₀ для каждого r1<…< i_r d существует ровно две различные точки покоя, которые удовлетворяют соотношениям:

.

4.Существует ровно одна симметричная точка покоя, которая удовлетворяет соотношению . Других точек покоя не имеется.

5.При dk в параметрической области >₀ и >₀ все точки покоя являются неустойчивыми.

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

8.2. Структура предельных циклов M₁(n,k)
Eсли dk, то

точка

, лежит в области притяжения собственного устойчивого предельного цикла, в котором для любого l=1,…,d выполняются соотношения

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

Следующая процедура дает алгоритм выхода на целую группу предельных режимов функционирования канонических ГГС произвольного класса, в том числе и для M₁(n,k).

Пусть имеется ориентированный орграф G(U,D), где U – множество вершин, D – множество дуг. Будем считать, что все вершины из U пронумерованы в некотором порядке от 1 до n. Обозначим через D_i множество номеров вершин, из которых дуги входят в вершину u_i: D_i={j(u_j,u_i)D}. Функцию соответствия вершине ее номера в выбранной нумерации, обозначим через I(u). Далее, разобьем множество вершин на непересекающиеся непустые подмножества: U=U₁…U_k, k1, и обозначим его через . Пусть _i обозначает множество пар _i={(s,U_sD_i)U_sD_i}, где U_sD_i -мощность множества U_sD_i, i=1,…,n., и пусть для разбиения  выполняется следующее условие: для всех j=1,…,k, и всевозможных двух вершин u и v из U_j имеем _i₍_u₎ _i₍_v₎. Возьмем некоторое -разбиение вершин графа G(U,D), которое удовлетворяет введенным выше условиям. По данному разбиению построим орграф G(V,W), в котором вершинами являются множества U_j, j=1,…,k, а дугами входящими в Uj являются дуги, исходящие из вершин, номера которых стоят первом месте пар, входящих в множество _j. Назовем этот граф -порожденным, а операцию порождения обозначим через P(G(U,D), ).

Обозначим каноническую ГГС i-го класса, построенную на орграфе G(U,D), через S_i(G(U,D)), а систему уравнений, соответствующую S_i(G(U,D)) – через M_i(G(U,D)) с переменными x_j, j=1,…,n. Заменим все x_j в каждом уравнении системы M_i(G(U,D)) на соответствующее y_l ,где l является номером множества U_l, содержащего x_j. В результате, как легко убедиться, система M_i(G(U,D)) преобразуется в систему, в которой имеется по U_lидентичных дифференциальных уравнений относительно y_l, l=1,…,k. Возьмем по одному такому уравнению и составим новую систему, которую обозначим через M_i(P(G(U,D), )). ГГС i-го класса, которая описывается данной системой будем обозначать через S_i(P(G(U,D), )) и будем называть -порожденной.

Имеется следующий простой алгоритм построения предельного режима функционирования S_i(G(U,D)) по заданному предельному режиму функционирования S_i(P(G(U,D), )). Пусть векторная функция (y₁,…,y_k) является предельным режимом функционирования M_i(P(G(U,D), ))). Это может быть точка покоя, цикл или аттрактор более сложной природы. Он может быть устойчивым или неустойчивым по своей природе. Зададим векторную функцию (x₁,…,x_n) следующим образом: . Тогда вектор-функция (x₁,…,x_n) будет представлять предельный режим функционирования того же типа канонической ГГС S_i(G(U,D)).
Для ГГС первого класса S₁(G(U,D)) мы предполагаем, что вектор-функция (x₁,…,x_n), представляя глобально неустойчивый предельный режим функционирования, будет устойчивым на многообразии .

Как легко заметить, симметричные структурные графы G(n,k) имеют следующие -разбиения. Пусть n=n₁n₂, n₁>k. Тогда следующее разбиение множества вершин орграфа G_n_,_k на n₁ подмножеств U_j={j,j+n₂,…,j+(n₁-1)n₂}, j=1,…,n₁, является -разбиением. Это позволяет легко строить неустойчивые режимы функционирования канонических симметричных ГГС сетей. Тем не менее, хотя последнее сильно облегчает проблему описания всех режимов функционирования канонических симметричных ГГС сетей, полностью ее не решает, так как нельзя исключить, что в системах могут быть и другие неустойчивые режимы функционирования, которые нельзя построить с помощью описанной процедуры.

Заключение

В работе описывается класс теоретических объектов - ГГС, исследование свойств которых представляется важным для понимания закономерностей функционирования природных генных сетей, которые обеспечивают выполнение практически всех жизненно важных функций организмов. Предложен математический аппарат для исследования ГГС. Показано, что ГГС могут иметь устойчивые стационарные и циклические режимы функционирования, а также аттракторы более сложной природы. Численные расчеты показывают, что для канонических ГГС при достаточно больших  и  их предельные свойства определяются только структурно-функциональными связями, которые устанавливаются посредством регулирования активности одних генетических элементов другими. Последнее выразилось в формулировании ряда эмпирических критериев, которые связывают предельные свойства канонических гипотетических генных сетей со свойствами соответствующих структурных графов. Критерии выполняются при некоторых достаточно больших значениях параметров моделей. Применительно к реальным генным сетям это означает, что необходимым условием существования в них необходимого количества стационарных и/или циклических режимов является наличие определенной минимальной сложности (нелинейности) процессов регуляции активности генетических элементов сети. Необходимая сложность может достигаться мультимеризацией белков-репрессоров и/или наличием достаточно большого числа промежуточных стадий, ведущих от гена к белку-регулятору.

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

Благодарности: Авторы благодарят Н.А. Колчанова, А.В. Ратушного, Г.В. Демиденко, Кагая В.В. за плодотворное обсуждение результатов работы. Работа получила частичную финансовую поддержку грантов РФФИ 02-04-48802-а, 02-07-90359-в, 01-07-90376-в и Междисциплинарного интеграционного проекта СО РАН "Моделирование фундаментальных генетических процессов и систем".

СПИСОК ЛИТЕРАТУРЫ.

Kolchanov N.A. (1997). Регуляция транскрипции генов эукариот: базы данных и компьютерный анализ. Mol. Biol. (Msk.), 31, 581–583.
Kolchanov N.A., Anan'ko E.A., Kolpakov F.A., Podkolodnaya O.A., Ignat'eva E.V., Goryachkovskaya T.N., Stepanenko I.L. (2000) Gene Networks Mol. Biol. (Msk.), 34, 449-460.
Thomas R., Thieffry D., Kaufman M. (1995) Dynamycal behavior of biological regulatory networks-I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state. Bulletin of Mathematical Biology, 57, 247-276.
Johnson AD, Poteete AR, Lauer G, Sauer RT, Ackers GK, Ptashne M.lambda Repressor and cro--components of an efficient molecular switch. Nature 1981 294(5838):217-23. Sugano M., Okamatsu H., Ide T. // Agricultural and Biological Chemistry. 1978. V. 42(11). P. 2009-2014.
Ptashne. M. A Genetic Switch: Phage Lambda and Higher Organisms. Oxford: Blackwell Scientic, 1992. 217 p.
Gardner TS, Cantor CR, Collins JJ. (2000) Construction of a genetic toggle switch in Escherichia coli. Nature, 403, 339-342.
Elowitz MB, Leibler S. (2000) A synthetic oscillatory network of transcriptional regulators. Nature. 403, 335-338.
Харари Ф. Теория графов. М.: Мир, 1973, 300с.

Рисунки к статье.

Рис.1 Диаграммы и орграфы гипотетических генных сетей

Рис.2 Простейшая гипотетическая генная сеть.

Рис.3 Гипотетическая генная сеть, имеющая циклический режим функционирования

Рис.4. Пример генной сети, составленной из семи генетических элементов.

Рис.5. Пример генной сети, составленной из девяти генетических элементов. При >3, >3 численно выявляется одна устойчивая точка покоя и один незатухающий режим функционирования, предположительно квазицикл.