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

На правах рукописи

ГУРСКИЙ Виталий Валериевич

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ГЕННЫХ СЕТЕЙ В РАННЕМ РАЗВИТИИ ДРОЗОФИЛЫ

Специальность 05.13.18—-математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Санкт-Петербург - 2009

003464229

Работа выполнена в Секторе теории твердого тела Учреждения Российской академии наук «Физико-технический институт им. А. Ф. Иоффе РАН»

Научный руководитель:

Доктор физико-математических наук, профессор А. М. Самсонов

Официальные оппоненты:

Доктор физико-математических наук, профессор А. П. Киселев

Кандидат физико-математических наук С. А. Руколайне

Ведущая организация: Федеральное государственное образовательное

учреждение высшего профессионального образования «Московский государственный университет им. М. В. Ломоносова»

Защита диссертации состоится « I (Л 2009 г. в 1 ^ часов на заседании

диссертационного совета Д 212.229.13 Санкт-Петербургского государственного

политехнического университета по адресу: 195251 Санкт-Петербург, Политехническая ул., д. 29, корпус £ аудитория

Ж

Автореферат разослан « » ч^^ ^РчЛЛ 2009 г.

Ученый секретарь диссертационного совета Д 212.229.13

. 'Г/

доктор технических наук, профессор /■'у^-— Б- С. Григорьев

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

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

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

Научная новизна результатов исследования складывается из следующих фактов:

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

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

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

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

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

На защиту выносятся следующие положения:

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

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

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

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

Апробация работы. Основные результаты диссертации докладывались на Международных конференциях "Days on Diffraction" (С.-Петербург, 2000, 2001; 2003; 2004), на Европейских междисциплинарных школах по нелинейной динамике для систем и анализу сигналов "Euroattractor" (Варшава, Польша, 2000; 2002), на Международной междисциплинарной школе по моделированию транскрипционных регуляторных сетей (Амблетез, Франция, 2002), на Московской международной конференции по

вычислительной молекулярной биологии "МССМВ-2003" (Москва, 2003), на Международной ежегодной конференции BMES (США, 2003), на Политехнических симпозиумах «Молодые ученые - промышленности Северо-Западного региона» (С.Петербург, 2004; 2007), на Международной конференции по системной биологии «From Molecules & Modeling to Cells (SysBio 2005)» (Гозау, Австрия, 2005), на Международной конференции по вычислительной системной биологии "WCSB 2005" (Тампере, Финляндия, 2005), на Международных симпозиумах по сетям в биоинформатике "ISNB" (Амстердам, Голландия, 2006; 2007), на Международной конференции по биоинформатике регуляции и структуры генома "BGRS 2006" (Новосибирск, 2006), на Международной конференции по нанобиотехнологиям (С.-Петербург, 2006), на Международных симпозиумах "Atomic Cluster Collisions: Structure and Dynamics from the Nuclear to the MesoBioNano scales" (ISACC'07, Дармштадг, Германия, 2007; ISACC'08, С.-Петербург, 2008), на Российско-Немецкой конференции по системной биологии (Москва, 2008), на Международной конференции по системной биологии «ICSB-2008» (Гетеборг, Швеция, 2008), на семинаре кафедры зоологии Кембриджского университета (Кембридж, Великобритания, 2008), на семинаре кафедры прикладной математики и статистики Университета Стони Брук (Стони Брук, США, 2003), на семинарах ФТИ им. А.Ф. Иоффе РАН (С.-Петербург, 2005-2006 гг.), на семинарах отдела компьютерной биологии СПбГПУ (С.-Петербург, 2002-2008 гг.).

Публикации. По теме диссертации опубликовано 6 статей в реферируемых научных журналах, одобренных ВАК, и 4 статьи в научных сборниках (без учета аннотаций докладов на конференциях).

Структура диссертации. Диссертация состоит из введения, четырех глав, заключения, списка цитируемой литературы и приложения. Основной текст работы содержит 148 страниц, включая 34 рисунка, 12 таблиц, заключение и список литературы. Приложение содержит 20 страниц, включая 8 рисунков.

Краткое содержание работы

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

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

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

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

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

Континуальные уравнения, описывающие динамику концентраций иа(х,1) белков, кодируемых генами а в момент I в точке х на центральной оси эмбриона, выводятся из закона сохранения массы и закона Фика для диффузии и имеют вид:

где и = u(x,t) — вектор концентраций и" = ua(x,t) для N белков, 5я — скорость синтеза а-го белка, Р* — скорость его химического распада, и If — коэффициент диффузии а-го белка между ядрами. При постановке начально-краевой задачи система уравнений (1) снабжается начальными условиями и краевыми условиями Неймана.

Также исследуется дискретная модель для концентраций и'(I) белка от а-го гена в /-м ядре эмбриона:

=5; (и,, I) - г (и;) ■+ Ь° [«, - и;)+- и°)], а = 1,..., N,;=1,..., М(0, (2)

где через и, = u,(t) обозначен вектор концентраций белков от всех генов в ядре /. Коэффициенты D" и DP в (1) и (2) имеют разную размерность. Уравнения (2) не являются дискретной версией уравнений (1) по следующей причине. Клеточная структура эмбриона представлена в (2) явным образом через зависимость от номеров i ядер будущих клеток. Поэтому перестроение этой структуры в ходе деления ядер выражается в том, что число ядер М меняется со временем и, соответственно, число уравнений в (2). Таким образом, уравнения (2) представляют динамическую систему, размерность фазового пространства которой есть функция времени и удваивается в результате каждого деления ядер.

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

В рамках принятых предположений о функциях и уравнения континуальной модели (1) предложены в следующем виде [3, 4, 5, 7, 10]:

^ = xWadtT"bub+m"u™(x) + h°)-rU° + D"^-, a = l,...,N, (3) St J 8x

где N = 5 — число генов (белков), R" = const, монотонно возрастающая нелинейная функция g(y) = 0.5(1 +уу2) описывает плавное переключение уровня синтеза от нулевого значения до максимального и называется «функцией регуляции», аргумент функции g в скобках в (3) описывает сумму всех факторов, регулирующих синтез а-го белка. В число этих факторов входят концентрации всех белков uh(x,t), которые влияют на синтез с весом Т'1' (коэффициент взаимодействия между генами а и b: если ген b активирует ген а, то Тъ > 0; если репрессирует, то Т'ь < 0; если не оказывает никакого влияния, то Гь = 0), а также заданная концентрация материнского белка Bicoid «Bcd(x),

влияющая на синтез с весом т". Постоянная h" описывает среднее влияние остальных, не учтенных явно, факторов, слагаемое Хаи в (3) представляет линейное приближение для скорости распада белка Р", функция /(/) описывает влияние деления ядер в эмбрионе, которое приводит к изменению общего числа копий генов и, следовательно, изменяет среднюю скорость синтеза белков от этих генов. Уравнения (3) снабжаются начальными условиями и краевыми условиями Неймана на ограниченном интервале центральной оси эмбриона. В качестве объекта исследования выбраны гены hb, Кг, gt, kni и eve поскольку экспрессия этих генов может считаться независимой от других генов в рассматриваемом интервале центральной оси эмбриона.

Значения пятидесяти параметров р = {/?', Th, т", h", X", If), 1 < а, Ь < 5, в модели (3) априори не известны и были найдены путем решения обратной задачи моделирования в ходе численной минимизации следующего функционала:

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

Ранее активность описанных генов в эмбрионе дрозофилы исследовалась только в рамках дискретной модели (2), в которой явно учтена клеточная структура эмбриона. Однако оставалось неясным, насколько важна дискретная природа клеточной структуры эмбриона и деления ядер с точки зрения математического моделирования. С целью ответить на этот вопрос исследовалась континуальная модель (3), поскольку в такой модели клеточная структура эмбриона аппроксимируется континуальной средой. Были рассмотрены три варианта системы (3) с разными функциями /(/), моделирующими деление ядер в эмбрионе с разной степенью точности. В каждой из этих моделей значения параметров были получены минимизацией функционала (4). В результате решения обратной задачи моделирования во всех моделях были получены решения, описывающие экспериментальные функции концентраций с удовлетворительной точностью (рис. 1).

(4)

а* П

ч 10 /11

ю о 112 з

Кг ИЬ \ 9' кш ече

V / А

А \ / А 1 и м

* О 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 10 0.2 0.4 0 6 08 10 0.2 0.4 0.6 0.8 1

X

Рис. 1. Решения и° (пунктирные кривые) в континуальной модели (3) для трех разных функций х(0 в сравнении с экспериментальными данными (сплошные кривые). Названия белков указаны на графиках. Верхний ряд графиков соответствует раннему циклу развития эмбриона, нижний ряд — более позднему. Безразмерная пространственная координата изменяется вдоль центральной оси эмбриона.

Таким образом, показано, что явный учет клеточной структуры эмбриона и деления ядер не является необходимым при моделировании пространственно-временного распределения концентраций («паттернов») белков.

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

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

*й>4>4>4>- (5)

где М — заданное конечное число ядер в системе, 5Д? — число линейно устойчивых стационарных точек уравнений (2) при данном М, Л^ — число всех стационарных аттракторов, реализующихся в системе, в которой начальное число ядер выбиралось равным М12к. Таким образом, описывает число аттракторов в системе с к последовательными делениями и вычисляется путем решения уравнений (2) с большим количеством случайных начальных условий.

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

аттракторов. Этот вывод является следствием переменной размерности системы (2) во времени; он проиллюстрирован на рис. 2 в частном случае.

Рис. 2. Иллюстрация эффекта сокращения числа аттракторов в модели (2) для случая М = 4. Показаны два случая: в первом

1 деление 2 деления

система приходит к М ядрам в результате одного деления, во втором — двух делений. В обоих случаях число линейно устойчивых стационарных точек составляет 5д-у = 16

(звездочки на рисунке) В первом случае из 16 стационарных точек только 4 реализуются как аттракторы (закрашенные ^^^^ ^ ^ ^ ^ звездочки) при варьировании начальных условий в модели; во ^^^^ ^^^^ втором — только 2. Остальные стационарные состояния (не

закрашенные звездочки) недостижимы для динамики в модели при любых начальных условиях в обоих случаях.

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

Рассмотрена дискретная модель экспрессии шести генов, основанная на следующих уравнениях [9]:

А

du°(t) dt

= R"g{jtT°buî -«,") +«+1-о], (6)

где белки нумеруются через а = 1,..., 6, ядра в эмбрионе нумеруются через i = 1,..., M(t). Решения уравнений (6) биологически обоснованны во временном интервале 0 < t < т, где г определяется экспериментально измеренной длительностью биологически важного промежутка времени. Уравнения (6) снабжаются фиксированными начальными условиями.

Значения параметров в (6) получались путем решения обратной задачи моделирования, то есть в ходе минимизации функционала, определяющего квадрат разности между решением уравнений (6) и экспериментальными данными. С целью улучшения параметрической устойчивости решений в модели (6) минимизация функционала осуществлялась с дополнительным требованием квазистационарности решений на временах t> т. Устойчивость решений оценивалась вычислением отклонения возмущенных решений от экспериментальных данных в норме ¿2-

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

правильную динамику в биологически важный период времени 0 < / < г (рис. 3). Для большинства параметров их значения в наборах р:„ з = 1,2,3, близки друг к другу, однако значения некоторых параметров отличаются в /э, довольно значительно.

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

минимизации функционала качества, в сравнении с экспериментальными данными (пунктирные кривые) для трех

О 0.5 1 0 0.5 1 0 0.5 ^ 1 0 0.5 1 0 0.5 1 0 0.5 1

моментов времени. По оси абсцисс отложены масштабированные

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

Была исследована устойчивость решений при найденных значениях параметров к возмущениям трех типов: возмущениям начальных условий, возмущениям решений и° в нескольких фиксированных моментах времени / > 0, и возмущениям значений параметров. Характеристики возмущений первых двух типов были оценены по разбросу экспериментальных данных в ходе анализа концентраций всех белков во многих эмбрионах.

Возмущенные решения в модели (6) с р = х = 1,2,3, сравнивались с возмущенными решениями в модели (6) со значениями параметров р, полученными ранее в результате решения обратной задачи без дополнительного требования квазистационарности решения на больших временах. Для возмущений всех трех типов решения в модели с р =р, демонстрировали более устойчивое поведение.

В качестве примера исследования параметрической устойчивости опишем схему и результаты возмущений решений на временах / > 0. В интервале 0 < I < т фиксировались 24 момента времени 1-„ в которые решения возмущались следующим образом:

где и° — невозмущенное решение, (м,")рс„ — возмущенное решение, £ — случайная

переменная с нормальным распределением е ~ Щ0,а) с нулевым средним и фиксированным стандартным отклонением а. Способ возмущения решений в (7) был

выбран в качестве простой модели логнормального распределения концентраций белков, наблюдаемого в эксперименте. Анализ экспериментальных кривых дал приблизительный разброс 0.5 < а < 1 для разных белков. В каждый момент и и для каждого фиксированного а производилось 500 случайных возмущений согласно формуле (7). Результаты вычислений представлены на рис. 4.

г = 0.1

£ 1

Г = 0.5

= 1.0

¡6

3 4

X г га

& 3

5

6

1.

<"0.5

АЗ .7 5.* 4. 6 : Л'2

«г -м

40 60 £

номер ядра /

Рис. 4. Сравнение решений, возмущенных по формуле (7), в моделях двух типов, (а): Возмущенное решение для белка Кг в модели (6) со значениями параметров из наборов р\-ръ для трех значений а. На графиках показаны средние значения концентрации (и')^ для Кг и

стандартные отклонения в каждом ядре в момент ( = т. (б)-. То же, что в (а), но для значений параметров, полученных решением обратной задачи без дополнительного требования квазистационарности решений на временах / > г. (в): Зависимость среднеквадратичного отклонения возмущенных решений от экспериментальных данных (в долях от такого же отклонения невозмущенных решений) от степени Е квазистационарности невозмущенных решений на временах ( > т при а = 0.5 в моделях двух типов. Е вычисляется как среднеквадратичное отклонение невозмущенных решений при < = г от аттрактора в модели. Через А1-АЗ обозначены результаты для значений параметров Р\~Ръ, цифрами 1-7 на графике обозначены результаты для семи известных наборов значений параметров, полученных решением обратной задачи без дополнительного требования квазистационарности решений на временах (> т.

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

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

X, % от длины эмбриона х, % от длины эмбриона

Рис. 5. Эффект малой позиционной вариабельности решений на примере белка НЬ: йс1" « ¿х8"1.

(а): Экспериментальные функции и®"^х), полученные из 45 эмбрионов дрозофилы. Позиционная вариабельность функции обозначенная через йх®"1, вычисляется как разброс точек пересечения фиксированного значения концентрации с данными функциями, где фиксированное значение равно значению «усредненной» функции 1^с\х) в точке, отмеченной пунктирной линией.

(б): Решения для белка НЬ при (= т в модели (6), в которой в качестве ы0"1^) поочередно взяты 45 функций из (а). Величина определяется как разброс пространственных положений точек перегиба решений из (б). Везде вместо дискретного индекса / (номер ядра) используется континуальная переменная х, равная положению ядра на центральной оси эмбриона (в % от длины эмбриона).

Позиционная вариабельность решений была исследована в континуальном приближении и в пределе ступенчатой функции регуляции: g = Н, где Н(у) — функция Хевисайда. Получены аналитические стационарные решения в континуальной модели экспрессии четырех генов, а именно решения следующей системы уравнений [8]:

К-Н [ У Г V (х) + т'иш (х) ■+ V" (х) + И"| - Л V* (х) + В" = о, (8)

) Ох

где индекс а = 1,...,4 нумерует концентрации белков от рассматриваемых генов, 1"{х) — внешний фактор, отличный от ыБс<1(х), остальные обозначения описаны выше. Уравнения (8) снабжаются краевыми условиями Неймана на ограниченном интервале на центральной оси эмбриона (интервал 0 <х < 1 для безразмерной пространственной переменной).

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

«"(*,?) = Со (*,?) + с? (х, д) ехр(у°х) + с°2 (х, д)ехр(-/°х), а = 1,..., 4, (9)

где через ^ = обозначен вектор положений д, всех 5 точек перегиба решений для

всех а, у" = лМ°/£>", и коэффициенты с°к(х,д) суть кусочно-постоянные функции от х. Значения д, были фиксированы равными пространственным положениям точек

перегиба функций концентраций белков из эксперимента. Из графиков для белков НЬ, Кг, Gt и Kni из нижнего ряда на рис. 3 можно определить S = 10 таких положений, нумеруемых в порядке их следования на единичном интервале переменной х. Можно показать, что по построению решений (9) положения q¡ должны удовлетворять равенствам:

Yj*'xku,'(q,,q) + m,-wu*"i(ql) + V(q¡) + h"w j = l,..,10, (10)

к

где a(j) есть однозначная функция s —► a(s), связывающая номер s положения q¡ с номером (или названием) белка а, чье стационарное решение имеет точку перегиба в q¡ (например, а(3) = "НЬ"; положение пунктирной линии на рис. 56 совпадает с qi). Функции ub(x,q) в (10) суть решения (9).

При фиксированных q¡ равенства (10) можно интерпретировать как уравнения для значений параметров р = {i?1, Ть, mН\ Xa, Lf}. Были найдены несколько наборов рк значений параметров, удовлетворяющих равенствам (10), с помощью метода разностной эволюции. Таким образом, при значениях параметров рк стационарные решения (9) существуют с точками перегиба в q¡, то есть, эти решения качественно совпадают с экспериментальными функциями концентраций рассматриваемых белков.

Варьирование (10) приводит к линейной алгебраической системе уравнений для вариаций Sqs вида:

f!(q,p,A,l)Sq, + f?(q,p,A,l,SA,Sl) = 0, S = l,...,10; (И)

где А и I — параметры аппроксимации функции wBcd(x): «Bcd(*) = Аехр(-/х); и /J — известные функции.

Система (11) была решена при найденных значениях параметров р = p¡, и при SA = а4, 81 = с/, где а4 и а1 суть стандартные отклонения во множестве значений А и /, соответственно, полученных из 45 экспериментальных функций nBcd(x). Полученные решения Sqs были использованы для вычисления так называемых коэффициентов позиционной вариабельности K¡ для каждой точки перегиба: К, = Sq,/SxBcá(qs), s = 1,...,10, где йхВс<| оценивалось из эксперимента.

Из всех наборов значений параметров p¡t три набора (p¡- р3) соответствовали очень малым значениям коэффициентов K¡. В частности, для положения q¡ точки перегиба решения для белка НЬ (пунктирная линия на рис. 56) соответствующий коэффициент фильтрации при трех наборах значений параметров pi-p3 равен: А'з = 0.009, 0.112 и 0.217. Эти результаты согласуются с экспериментальными оценками.

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

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

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

- исходя из экспериментальных данных, пространственный интервал разбивается на отрезки (О.!, 05+0 преимущественной локализации одного белка (рис. 6);

л:

- исходя из анализа литературных данных, предполагаются следующие ограничения на параметры 1аЪ, характеризующие тип взаимодействия между генами:

— в уравнениях (3) полагается 1У = 0 и £ = Н, где Н— функция Хевисайда.

В принятых предположениях система (3) расцепляется, и х становится параметром. В такой системе для каждого фиксированного * может быть 2Л' различных аттракторов с компонентами либо 0, либо К'ИС\ где N — число генов (белков). Для построения решений, качественно сходных с рис. 6, необходимо найти значения параметров, при которых в каждом интервале (<2а, &+]) реализуется аттрактор с компонентой /У/!0 (максимальное значение концентрации) для соответствующего белка и 0 для остальных белков. Доказано, что для того, чтобы в модели реализовывались аттракторы только с одной компонентой, равной КЧХ", и нулевыми остальными компонентами, достаточно выполнения следующих ограничений на параметры:

Рис. 6. Экспериментальные функции концентрации белков Hb, Кг, Gt, Kn¡ и Til в момент t = т. Пространственный интервал условно разбит на интервалы (Q¡, Q,+1) преимущественной локализации только одного белка, название которого приведено для

ü Q1 Q2 Q3 Q4 Ó5 l каждого интервала.

T°°> 0, Г'<0, Va,А;

(12)

Г" —+rb¿7-+m°uBal(x) + ha <0 Va*b,Vx.

go jí \ '

R-

Обозначим через ¿'"(л,/) = ^ Т°>'ик(х,1) + т°иьл(х) + И' аргумент функции

регуляции в модели (3), и = . С учетом (12) доказывается, что для того,

чтобы в каждом (£>,, £?)-и) выполнялось и" —► /?'/л" при I —> °о для соответствующего белка а и и4 —> О для всех остальных Ь, достаточно потребовать, чтобы 5,„(л)> 0 внутри ((¿¡, бл-О и $а(х)<0 вне этого интервала. Например, для белка НЬ (а = 2) соответствующие ограничения имеют следующий вид (см. также рис. 6):

So(x) > о, (0,0)u(a,a); 502(X) < о, xe(Q,Q4)u(Qs,l).

(14)

Аналогичные ограничения выписываются для остальных белков.

Приближенные значения параметров Л" и Lf вычислялись из требования равенства максимальных значений экспериментальных концентраций значениям Fi'/l" и путем подгонки стационарных решений вида (9) к экспериментальным функциям концентраций (рис. 6) в локальных пространственных областях. Приближенные значения остальных параметров (Vh, m", h") фиксировались согласно ограничениям (12), (13), (14) и аналогам (14) для других белков.

Рис. 7. (а): Решение в модели (3) при приближенных значениях параметров и 1 = т. (б): Решение в той же модели (сплошные кривые), но при значениях параметров, полученных в результате применения метода наискорейшего спуска с приближенными значениями «о параметров в качестве начальной точки, в сравнении с экспериментальными функциями концентраций (пунктирные кривые), (в) и (г): То же, что в (а) и (б), соответственно, но для дискретной модели (6). Концентрация TU не показана в силу того, что в дискретной модели она (вместе с концентрацией Cad) использовалась в качестве заданного внешнего фактора.

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

в процессе оптимизации. Результирующее решение показано на рис. 76 в сравнении с экспериментальными данными. Метод также применялся к решению уравнений дискретной модели (6); были получены удовлетворительные результаты (рис. 7в,г).

В Заключении сформулированы результаты и основные выводы работы.

Выводы

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

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

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

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

Основные результаты диссертации опубликованы в следующих работах:

1. Samsonov, A.M. Exact solutions to a nonlinear reaction-diffusion equation and hyperelliptic integrals inversion [Текст] / A.M. Samsonov, V.V. Gursky // Journal of Physics A: Mathematical and General. - 1999. - V. 32. - C. 6573-6588.

2. Gursky, V.V. How gap genes make their domains: An analytical study based on data driven approximations [Текст] / V.V. Gursky, J. Reinitz, A.M. Samsonov // Chaos. - 2001. -V. 11. - № 1. - C. 132-141.

3. Gursky, V.V. Mathematical modeling of pattern formation due to the segmentation gene expression in Drosophila [Текст] / V.V. Gursky, J. Reinitz, A.M. Samsonov // Proceedings of the 1st European Interdisciplinary School on Nonlinear Dynamics for Systems and Signal Analysis "Euroattractor 2000". - Germany: Pabst Science Publishers, 2001.-C. 540-544.

4. Gursky, V.V. Pattern formation and nuclear divisions are uncoupled in Drosophila segmentation: Comparison of spatially discrete and continuous models [Текст] / V.V. Gursky, J. Jaeger, K.N. Kozlov, J. Reinitz, A.M. Samsonov II Physica D: Nonlinear Phenomena. - 2004. - V. 197. - C. 286-302.

5. Samsonova, M.G. A survey of gene circuit approach applied to modelling of segment determination in fruit fly [Текст] / M.G. Samsonova, A.M. Samsonov, V.V. Gursky, C.E. Vanario-Alonso // Multiple aspects of DNA and RNA: from Biophysics to Bioinformatics. -Elsevier, 2005.-C. 305-323.

6. Gursky, V.V. Cell divisions as a mechanism for selection in stable steady states of multi-stationary gene circuits [Текст] / V.V. Gursky, K.N. Kozlov, A.M. Samsonov, J. Reinitz // Physica D: Nonlinear Phenomena. -2006. -V. 218. -№ 1. -C. 70-76.

7. Самсонова, М.Г. Системный подход к исследованию развития организмов [Текст] / М.Г. Самсонова, В.В. Гурский, К.Н. Козлов, A.M. Самсонов // Научно-технические ведомости СПбГТУ. - 2006. - № 2. - С. 222-234.

8. Gursky, V.V. Approximate stationary attractors in Drosophila gap gene circuits in the limit of steep-sigmoid interactions [Текст] / V.V. Gursky, A.M. Samsonov, J. Reinitz // Proceedings of the 5th International Conference on Bioinformatics of Genome Regulation and Structure. - Novosibirsk: Institute of Cytology and Genetics, 2006. - V. 2. - C. 125-127.

9. Гурский, В.В. Модель с асимптотически устойчивой динамикой для сети генов gap в дрозофиле [Текст] / В.В. Гурский, К.Н. Козлов, A.M. Самсонов, Дж. Рейниц // Биофизика. - 2008. - Т. 53. - № 2. - С. 235-249.

10. Samsonov, A.M. On modelling of gene expression patterns in Drosophila embryo by gene circuit method [Текст] / A.M. Samsonov, V.V. Gursky, K.N. Kozlov, J. Reinitz // Proceedings of the 2nd International Symposium "Atomic Cluster Collisions: structure and dynamics from the nuclear to the biological scale". - London: Imperial College Press, 2008. -C. 426-440.

Лицензия ЛР № 020593 от 07.08.97

Подписано в печать 16.02.2009. Формат 60x84/16. Усл. печ. л. 1,0. Уч.-изд. л. 1,0. Тираж 100. Заказ 0047.

Отпечатано с готового оригинал-макета, предоставленного автором, в типографии Издательства Политехнического университета. 195251, Санкт-Петербург, Политехническая ул., 29.

Введение '

1 Биологические определения и модели генной экспрессии

1.1 Генная экспрессия в многоклеточных организмах.

1.2 Описание системы сегментационных генов в растущем эмбрионе дрозофилы

1.3 Вывод континуальных уравнений.

1.4 Генная сеть.

1.5 Учет деления ядер в математической модели

1.6 Уравнения дискретной модели.

2 Континуальная модель сети сегментационных генов hb, Кг, gt, kni и eve. Роль клеточной структуры и деления ядер

2.1 Континуальная модель экспрессии сегментационных генов Кг, hb, gt, kni и eve в дрозофиле.

2.1.1 Уравнения модели.

2.1.2 Получение значений параметров путем решения обратной задачи математического моделирования.

2.1.3 Описание экспериментальных данных

2.1.4 Описание метода минимизации функционала качества.

2.1.5 Результаты оптимизации.

2.1.6 Решение в модели с оптимальными значениями параметров

2.2 Динамическая роль деления ядер.

2.2.1 Эффект сокращения числа аттракторов.

2.2.2 Ограничение от правила деления.

2.2.3 Ограничение от динамики в ранних циклах.

3 Дискретная модель сети генов gap. Исследование механизмов параметрической устойчивости

3.1 Модель с асимптотически устойчивой динамикой для сети генов gap

3.1.1 Описание модели.

3.1.2 Оптимизация параметров с дополнительным требованием квазистационарности решений на больших временах.

3.1.3 Обоснование использования предположения о квазистационарности на основе анализа экспериментальных данных.

3.1.4 , Результаты моделирования.

3.1.5 Устойчивость к возмущению начальных условий.

3.1.6 Устойчивость к возмущению решений в отдельные моменты времени

3.1.7 Устойчивость к возмущению значений параметров.

3.2 Устойчивость границ областей экспрессии при вариации функции itBcd(:r)

3.2.1 Позиционная вариабельность в эксперименте и в модели.

3.2.2 Упрощающие предположения и аналитические стационарные решения

3.2.3 Уравнения для вариаций границ кинков при вариации -uBcd(x)

3.2.4 Вычисление коэффициентов позиционной вариабельности.

4 Комбинированный метод решения обратной задачи моделирования

4.1 Анализ упрощенных уравнений.

4.2 Контроль начальной стадии эволюции решения типа слоистого паттерна

4.3 Применение к сети gap-генов в дрозофиле.

4.3.1 Получение приближенных значений параметров

4.3.2 Уточнение значений параметров путем оптимизации методом оптимального наискорейшего спуска.

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

Физико-математическое моделирование биологических процессов переживает период бурного развития, тесно связанный с возникновением новых междисциплинарных научных областей, таких как биоинформатика и системная биология [72, 13]. Суть данных научных областей состоит в применении методов математического моделирования и вычислительных алгоритмов к экспериментальным биологическим данным для описания процессов, происходящих в живых системах (в клетках, тканях или организмах).

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

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

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

Другой важной задачей математического моделирования в системной биологии является математическое описание механизмов, ответственных за устойчивое развитие организмов (робастность) [41]. Большинство процессов в клетках происходят на фоне "шумов" различной природы. Данные шумы могут определяться как внутренними причинами, связанными со стохастической природой фундаментальных молекулярных процессов, так и внешними, связанными с флуктуацнями внешней среды. Адекватное математическое описание механизмов биологической робастности является актуальным для увеличения предсказательной силы моделей.

Одним из наиболее изученных многоклеточных организмов является плодовая мушка дрозофила [42]. Получено большое количество новых экспериментальных данных по концентрациям белков в эмбрионе дрозофилы, которые требуют интерпретации и количественного объяснения в рамках математических моделей [65, 20]. Актуальной в последние годы также стала дискуссия о степени точности позиционирования сегментов эмбриона в ходе раннего развития дрозофилы в условиях изменчивых внешних факторов [33, 26]. За формирование этих сегментов отвечают специальные гены, называемые генами сегментации, или сегментационпыми генами. Понимание робастности, наблюдающейся в эксперименте, требует математической интерпретации.

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

Научная новизна результатов диссертации складывается из следующих фактов:

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

2. Впервые описан математический механизм сокращения числа аттракторов в динамических системах, в которых размерность фазового пространства изменяется во времени.

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

4. Предложен новый способ подбора начальной точки в пространстве параметров для задачи оптимизации значений параметров в модели активности генов сегментации.

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

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

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

Уравнения модели основаны на балансе между тремя процессами: синтезом белка, его диффузией в эмбрионе и химическом распадом молекул белка. Существуют различные методы формализации этих процессов и, следовательно, различные виды уравнений [63, 17]. Соответствующие модели можно условно разделить на три основных класса, связанных с тремя различными математическими подходами к формулировке уравнений. Первый класс ("логические" модели) включает модели генной регуляции, уравнения которых основаны па булевой логике и ее различных модификациях [67, 23, 39, 11, 66]. В логических моделях состояние гена описывается булевой переменной, значение 1 которой соответствует состоянию гена, в котором он активен (кодирует белок), а значение 0 — состоянию, в котором ген неактивен (белок не синтезируется). Существуют также обобщения логических моделей, в которых логическая переменная может иметь более, чем два дискретных значения. Динамика такой переменной оштсывпется с помощью булевых функций. Логические модели позволяют лишь качественно описывать экспрессию генов, то есть решения для функций концентрации белков, кодируемых генами, представляют собой дискретные уровни в пространстве и времени.

Во втором классе моделей концентрации белков (продуктов экспрессии генов) являются непрерывными функциями времени и дискретными функциями пространственной переменной; динамика таких функций описывается с помощью обыкновенных дифференциальных уравнений (см., например, [6, 13, 62]). К этому же классу можно отнести модели, в которых концентрации белков также непрерывны по пространству и являются решением уравнений в частных производных (см., например, [46, 6, 62, 43]). Функции синтеза и распада белка в уравнениях строятся на основе имеющихся знаний об общих механизмах процесса экспрессии генов и химической кинетики. Диффузия моделируется как правило на основе закона Фика или дополнительно с учетом возможных нелокальных эффектов.

В третий класс входят стохастические модели экспрессии генов [52, 38, 40]. Данные модели являются наиболее детальными, поскольку сам процесс экспрессии гена осуществляется на основе стохастических молекулярных механизмов. Стохастический подход наиболее важен в ситуациях, когда биологические функции белков, кодируемых генами, осуществляются при очень маленьких концентрациях этих белков, поскольку именно при малых концентрациях стохастические эффекты значительны. В практике построения стохастических моделей генных сетей можно выделить два подхода. В первом подходе стохастические уравнения получаются просто добавлением в детерминированные уравнения случайной функции, определяющей стохастический характер процессов синтеза, распада и диффузии молекул белка. Экспрессия гена во втором подходе моделируется как набор элементарных реакций (присоединение к ДНК, диссоциация и гак далее) с заданными вероятностями случиться [22].

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

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

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

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

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

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

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

Во второй части главы исследуется другой вид параметрической устойчивости решений. Уравнения дискретной и континуальной моделей экспрессии генов сегментации зависят от заданного внешнего фактора uBcd(x), определяющего известную пространственную функцию, которая аппроксимирует распределение концентрации материнского белка Bed ("бикоид") вдоль центральной оси эмбриона. Изучается вариация решений модельных уравнений при вариации функции -иВсс1(ж). Выводится система алгебраических уравнений для вариаций 8ql пространственных положений (qi) точек перегиба решений. Решения Sqi этой системы определяют позиционную вариабельность решений модельных уравнений. Показывается, что Sqi << &rBcd, где <5xBcd есть вариация пространственного положения некоторого фиксированного значения функции iiBcd(a;). Этот результат совпадает с экспериментальными оценками.

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

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

Основные результаты диссертации докладывались на Международных семинарах "Day on Diffraction" (С.-Петербург, 2000, 2001; 2003; 2004), на Европейских школах по нелинейной динамике для систем и анализу сигналов "Euroattractor" (Варшава, Польша, 2000; 2002), на Международной школе по моделированию транскрипционных регуляторных сетей (Амблетез, Франция, 2002), па Московской международной конференции по вычислительной молекулярной биологии "МССМВ-2003" (Москва, 2003), на Международной ежегодной конференции BMES (США, 2003), на Политехнических симпозиумах "Молодые ученые — промышленности Северо-Западного региона" (С.Петербург, 2004; 2007), на Международной конференции по системной биологии "From Molecules h Modeling to Cells (SysBio 2005)" (Гозау, Австрия, 2005), на Международной конференции по вычислительной системной биологии "WCSB 2005" (Тампере, Финляндия, 2005), на Международных симпозиумах по сетям в биоинформатике "ISNB" (Амстердам, Голландия, 2006; 2007), на Международной конференции по биоинформатике регуляции и структуры генома "BGRS 2006" (Новосибирск, 2006), на Международной конференции по нанобнотехнологиям (С.-Петербург, 2006), на Международных симпозиумах "Atomic Cluster Collisions: Structure and Dynamics from the Nuclear to the MesoBioNano scales" (ISACC 2007, Дармштадт, Германия, 2007; ISACC 2008, С.-Петербург, 2008), на Российско-Немецкой конференции по системной биологии (Москва, 2008), на Международной конференции по системной биологии "ICSB-2008" (Гетсборг, Швеция, 2008), на семинаре кафедры зоологии Кембриджского университета (Кембридж, Великобритания, 2008), на семинаре кафедры прикладной математики и статистики Университета Стони Брук (Стони Брук, США, 2003), на семинарах ФТИ им.

А.Ф. Иоффе РАН (С.-Петербург, 2005-2006), на семинарах отдела компьютерной биологии СПбГПУ (С.-Петербург, 2002-2008).

Основные результаты диссертации опубликованы в работах [60, 29, 30, 27, 62, 28, 9, 31, 2, 61].

Заключение

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

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

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

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

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

Далее сформулированы основные выводы из полученных результатов.

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

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

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

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

Описанный в диссертации механизм выбора аттракторов может иметь важную интерпретацию в контексте биологической эволюции. Основываясь на проведенном математическом анализе, можно предположить, что соотношение расписания делений ядер и времени экспрессии генов сегментации в эмбрионе дрозофилы является неслучайным, но представляет собой результат эволюционного отбора. В ходе раннего развития в эмбрионе происходят 14 последовательных делений ядер, при этом белки, кодируемые генами сегментации, впервые появляются в ядрах только в 13-м цикле. Можно предположить, что в ходе эволюционного отбора получили преимущество особи с таким поздним появлением белков в том числе благодаря фактору, связанному с описанным механизмом выбора числа аттракторов в последнем 14-м цикле. Если бы белки появлялись раньше, то генная сеть была бы ограничена слишком малым количеством стационарных состояний в 14-м цикле, доступных для динамики концентраций белков. Особи с малым количеством таких состояний могут быть эволюционно ущербными, поскольку сложнее найти те значения параметров генной сети, которые гарантируют существование аттрактора с необходимым пространственным распределением белков в эмбрионе в конце 14-го цикла. С другой стороны, белки не появляются в 14-м цикле, поскольку, согласно предлагаемой гипотезе, в этом случае число разрешенных стационарных состояний было бы слишком большим. Существование таких особей может быть также невыгодно с эволюционной точки зрения, поскольку в математических терминах такие системы характеризуются фазовым пространством, "нарезанным" на очень большое количество областей притяжения соответствующих аттракторов. Такие системы, вообще говоря, должны быть параметрически менее устойчивыми, поскольку высока вероятность того, что при возмущении решения оно попадет в область притяжения другого аттрактора, и, сиедовательно, качественно изменится поведение решения на больших временах. Таким образом, появление белков в 13-м цикле при общем числе циклов, равном 14-и, можно рассмотреть как компромисс между скудостью и разнообразием динамических режимов функционирования генной сети в описанном выше смысле. Данный компромисс осуществляется с помощью механизма деления ядер (клеток) н описывается в терминах построенных математических моделей.

Из результатов первой части главы 3 следует, что требование квазистационарности решений на больших временах при решении обратной задачи моделирования позволяет улучшить параметрическую устойчивость решений. При исследовании модели вида (3.1)~(3.2), в которой значения параметров равны одному из десяти наборов (модели А\—Аз, Вг-В7), было показано, что наборы значений параметров, при которых решение в модели демонстрирует квазистационарное поведение после конца 14-го цикла (/ > т), соответствуют решениям, параметрически более устойчивым, чем наборы, при которых решение меняется значительно после конца 14-го цикла. Ограничение на поведение решений на больших временах может использоваться для решения обратной задачи моделирования при построении более робастных моделей в тех случаях, когда моделируемый процесс происходит в ограниченный промежуток времени и заканчивается до момента установления решений в модели.

Параметрическая устойчивость решений в построенных в диссертации моделях А\— Аз частично соответствует экспериментально полученным оценкам вариабельности концентраций белков, кодируемых генами сегментации. Данное соответствие имеет место для возмущений решений по формуле (3.7). Средне-квадратичное отклонение возмущенных по этой формуле решений в моделях А^Аз от усредненных (по эмбрионам) экспериментальных данных из рис. 3.8 для а = 0.1 п а = 0.5 находится в пределах разброса, наблюдаемого в разных эмбрионах в эксперименте [65, 20]. С другой стороны, возмущение начальных условий в А\-Аз приводит к отклонению решений от усредненных экспериментальных данных (рис. 3.6), которое больше, чем в эксперименте. Наконец, разброс решений при возмущении значений параметров (рис. 3.10) соответствует экспериментальным данным при а = 0.01.

Даже в тех случаях, когда параметрическая устойчивость решений в моделях А\-Аз хуже, чем в эксперименте, она всегда лучше, чем в моделях, в которых решение не является квазистационарным после конца 14-го цикла (модели В\~В7). В этих случаях использование требования квазистационарности при решении обратной задачи моделирования повышает параметрическую устойчивость решений в модели, но не в достаточной степени, и требуются дополнительные исследования других возможных механизмов, отвечающих за устойчивость.

Во второй части главы 3 показано, что, в согласии с экспериментом, разброс пространственных положений точек перегиба решении в модели экспрессии генов семейства gap оказывается значительно меньше, чем аналогичная позиционная вариабельность функции концентрации белка Bed (;/1!cd(ж)), являющейся внешним фактором в модели. Данный эффект "фильтрации ошибки позиционирования" математически описывается на уровне стационарных решений в модели со ступенчатой функцией регуляции и может быть охарактеризован как специальный вид параметрической устойчивости решений модели.

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

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

Метод хорошо работает как в дискретной, так и в континуальной моделях (рис. 4.4). Данный факт подтверждает выводы главы 2 о том, что оба подхода (дискретный и континуальный) одинаково применимы к моделированию экспрессии генов сегментации в дрозофиле. В более общем контексте метод может применяться к математическим моделям вида (2.1)—(2.3), в которых решение имеет вид "слоистого паттерна". Данный термин означает, что разные компоненты решения локализованы в пространственных областях, поч си непересекающихся друг с другом. Также очевидно, что метод применим к моделям с многомерной пространственной переменной. Единственным ограничением является наличие достаточного количества параметров, чтобы удовлетворить ограничения на параметры, сформулированные па первом шаге метода.

Автор глубоко благодарен научному руководителю А. М. Самсоттову за постановку задач и постоянное внимание к работе, а также всестороннюю помощь в науке и в жизни. Также автор выражает огромную признательность коллегам по научной работе: К. Н. Козлову, А. Матвеевой, Е. М. Мясниковой, А. Н. Писареву, Е. Г. Пустельниковой, М.Г. Самсоновой, С.Ю. Сурковой, Й. Йегеру, Ману, Д. Репницу. Автор искренне благодарен С. А. Вакуленко, чьи идеи лежат в основе материала, изложенного в главе 4 и во второй части главы 3.

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

2. В. В. Гурский, К.Н. Козлов, A.M. Самсонов, Дж. Рейниц, Модель с асимптотически устойчивой динамикой для сети генов gap в дрозофиле, Биофизика 53(2): 235249, 2008

3. И. Ф. Жимулев, Действие генов в раннем развитии дрозофилы, Соросовский образовательный журнал 7:30-34, 1998

4. К.Н. Козлов. Л. В. Петухов, М.Г. Самсонова, A.M. Самсонов, Применение метода наискорейшего спуска 1-го порядка для нахождения феноменологических параметров модели генных сетей, Тр. СПбГТУ. Прикладная математика, 485: 73-83, 2002

5. К. Н. Козлов, А. М. Самсонов, Новый метод обработки экспериментальных данных на основе теории оптимального управления, Журнал технической физики 73(11):6-14, 2003

6. Дж. Марри. Нелинейные дифференциальные уравнения в биологии. Лекции о моделях, М., Мир, 1983, 399 стр.

7. И. А. Молотков, С. А. Вакуленко. Сосредоточенные nejiunewibie волны, JL, Изд. Ленинградского гос. университета, 1988, 240 стр.

8. Э. Рис, М. Стернберг. Введение в молекулярную биологию. От клеток к атомам, М., Мир, 2002

9. М. Г. Самсонова, В. В. Гурский, К. Н. Козлов, A.M. Самсонов, Системный подход к исследованию развития организмов, Научно-технические ведомости СПбГТУ, 2:222-234, 2006

10. А. Н. Тихонов, А. А. Самарский. Уравнения математической физики, М., Наука, 1966

11. Р. Н. Чураев, А. В. Галимзянов, Моделирование реальных эукариотических управляющих генных подсетей на основе метода обобщенных пороговых моделей, Молекулярная биология 35(6): 1088-1094, 2001

12. Л. П. Шилышков, А. Л. Шильников, Д. В. Тураев, Л. Чуа. Методы качественной теории в нелинейной динамике. Часть I, Москва-Ижевск, Изд. Института компьютерных исследований, 2004

13. U. Alon. An Introduction to Systems Biology: Design Principles of Biological Circuits, Boston, Chapman&Hall/CRC, 2006

14. G. Bronner, H. Jackie, Control and function of terminal gap gene activity in the posterior pole region of the Drosophila embryo, Mechanisms of Development 35:205-211, 1991

15. J. Carr, R. L. Pego, Metastable patterns in solutions of ut = e2uxx 4- f(u), Comm. Pure Appl. Math. 42:523-576, 1989

16. K. W. Chu, Y. Deng, J. Reinitz, Parallel simulated annealing by mixing of states, J. Comput. Phys. 148: 646-662, 1999

17. H. de Jong, Modeling and simulation of genetic regulatory systems: a literature review, Journal of Computational Biology 9(1): 67-103, 2002

18. P. C. Fife, Mathematical Aspects of Reacting and Diffusing Systems, Lecture Notes in Biomathematics (Vol. 28), Springer Verlag, New York, 1979

19. Р. С. Fife, Diffusive waves in inhomogeneous media, Proceedings of Edinburg Mathematical Society 32: 291-315, 1989

20. FlyEx. Электронная база данных по экспрессии сегментационных генов в дрозофиле (http://urchin.spbcas.ru/flyex; http://flyex.ams.sunysb.edu/FlyEx)

21. М. Frasch, М. Levine, Complementary patterns of even-skipped and fushi tarazu expression involve their differential regulation by a common set of segmentation genes in Drosophila, Genes Sz Development 1:981-995, 1987

22. D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, Journal of Physical Chemistry 81(25): 2340-2361, 1977

23. L. Glass, S. A. Kauffman, The logical analysis of continuous, nonlinear biochemical control netwoiks, Journal of Theoretical Biology 39:103-129, 1973

24. T. Gregor, W. Bialek, R. R. de Ruyter van Stcveninck, D.W. Tank, E.F. Wieschaus, Diffusion and scaling during early embryonic pattern formation, Proceedings of the National Academy of Sciences 102(51): 18403-18407, 2005

25. T. Gregor, E.F. Wieschaus, A. P. McGregor, W. Bialek, D.W. Tank, Stability and nuclear dynamics of the Bicoid morphogen gradient, Cell 130(1): 141-152, 2007

26. T. Gregor, D.W. Tank, E.F. Wieschaus, W. Bialek, Probing the limits to positional information, Cell 130(1): 153-164, 2007

27. V. V. Gursky, J. Jaeger, K. N. Kozlov, J. Reinitz, A. M. Samsonov, Pattern formation and nuclear divisions arc uncoupled in Drosophila segmentation: Comparison of spatially discrete and continuous models, Physica D 197: 286-302, 2004

28. V. V. Gursky, K.N. Kozlov, A.M. Samsonov, J. Reinitz, Cell divisions as a mechanism for selection in stable steady states of multi-stationary gene circuits, Physica D 218(1): 70-76, 2006

29. V. V. Gursky, J. Reinitz, A.M. Samsonov, How gap genes make their domains: An analytical study based on data driven approximations, Chaos 11(1): 132-141, 2001

30. L. H. Harlwell, J. J. Hopfield, S. Leibler, A. W. Murray, From molecular to modular cell biology, Nature 402: C47-C52, 1999

31. B. Houchmandzadeh, E. Wiescliaus, S. Leibler, Establishment of developmental precision and proportions in the early Drosophila embryo, Nature 415: 798-802, 2002

32. B. Houchmandzadeh, E. Wieschaus, S. Leibler, Precise domain specification in the developing Drosophila embryo, Physical Review E 72: 061920, 2005

33. M. Howard, P. R. ten Wolde, Finding the center reliably: Robust patterns of developmental gene expression, Physical Review Letters 95: 208103, 2005

34. J. Jaeger, S. Surkova, M. Blagov, H. Janssens, D. Kosman, K. N. Kozlov, Manu, E. Myasnikova, С. E. Vanario-Alonso, M. Samsonova, D. H. Sharp, J. Reinitz, Dynamic control of positional information in the early Drosophila embryo, Nature 430:368-371, 2004

35. M. Kaern, Т. С. Elston , W. J. Blake, J. J. Collins, Stochasticity in gene expression: from theories to phenotypes, Nature Rev. Genet. 6:451-464, 2005

36. S.A. Kauffman. The Origins of Order: Self-Organization and Selection in Evolution, Oxford, Oxford University Press, 1993

37. Т. B. Kepler, Т. C. Elslon, Slochasticity in transcriptional regulation: origins, consequences, and mathematical representations, Biophysical Journal 81(6): 3116-3136, 2001

38. H. Kitano, Biological robustness, Nature Rev. Genet. 5:826-837, 2004

39. P. A. Lawrence. The Making of a Fly: The Genetics of Animal Design, Oxford, Blackwell Sci. Publ., 1992

40. R. E. Baker, E. A. Gaffney, P. K. Maini, Partial differential equations for self-organization in cellular and developmental biology, Nonlinearity 21: R251-R290, 2008

41. Mann, S. Surkova, A. V. Spirov, V. V. Gursky, H. Janssens, A.-R. Kim, O. Radulescu, С. E. Vanario-Alonso, D. H. Sharp, M. Samsonova, J. Reinitz, Canalization of gene expression in the Drosophila blastoderm by dynamical attractors, submitted

42. P. McHale, W.-J. Rappel, H. Levine, Embryonic pattern scaling achieved by oppositely directed morphogen gradients, Physical Biology 3:107-120, 2006

43. H. Mcinhardt. Models of Biological Pattern Formation, London, Academic Press, 1982

44. E. Mjolsness, Statistical mechanics of equilibrium molecular complex composition models, submitted

45. E. Mjolsness, D.H. Sharp, J. Reinitz, A connectionist model of development, J. Theor. Biol. 152:429-453, 1991

46. M. Mlodzik, W. J. Gehring, Hierarchy of the genetic interactions that specify the anteroposterior segmentation pattern of the Drosophila embryo as monitored by caudal protein expression, Development 101:421-435, 1987

47. E. Myasnikova, A. Samsonova, К. Kozlov, М. Samsonova, J. Reinitz, Registration of the expression patterns of Drosophila segmentation genes by two independent methods, Bioinformatics 17(1): 3-12, 2001

48. C. Nusslein-Volhard, H. G. Frohnhofer, R. Lehmann, Determination of anteroposterior polarity in Drosophila, Science 238:1675-1681, 1987

49. J. Paulsson, Models of stochastic gene expression, Physics of Life Reviews 2:157-175, 2005

50. E. Poustelnikova, A. Pisarev, M. Blagov, M. Samsonova, J. Reinitz, A database for management of gene expression data in situ, Bioinformatics 20:2212-2221, 2004

51. M. H. Plotter, H. F. Weinberger. Maximum Principles in Differential Equations, Prentice-Hall. Englewood Cliffs, NJ, 1967

52. J. Reinitz, D. Kosman, С. E. Vanario-Alonso, D. Sharp, Stripe forming architecture of the gap gene system, Developmental Genetics 23:11-27, 1998

53. J. Reinitz, E. Mjolsness, D.H. Sharp, Model for cooperative control of positional information in Drosophila by bicoid and hunchback, J. Exp. Zool. 271:47-56, 1995

54. J. Reinitz, D.H. Sharp, Mechanism of eve stripe formation, Mech. Dev. 49:133-158, 1995

55. A. M. Samsonov, Travelling wave solutions for nonlinear dispersive equations with dissipation, Applic. Analysis 57:85-100, 1995

56. A. M. Samsonov, On exact quasistationary solutions to nonlinear reaction-diffusion equation, Phys. Lett. A 245(6): 527-536, 1998

57. A. M. Samsonov, V. V. Gursky, Exact solutions to a nonlinear reaction-diffusion equation and hyperelliptic integrals inversion, J. Phys. A: Math. Gen. 32: 6573-6588, 1999

58. P. Smolen, D.A. Baxter, J.H. Byrne, Modeling transcriptional control in gene networks — methods, recent results, and future directions, Bulletin of Mathematical Biology 62: 247-292, 2000

59. G. Struhl, K. Struhl, P. M. Macdonald, The gradient morphogen Bicoid is a concentration-dependent transcriptional activator, Cell 57:1259-1273, 1989

60. S. Surkova, M. Samsonova, D. Kosman, K. N. Kozlov, Manu, E. Myasnikova, A. Samsonova, A. Spirov, С. E. Vanario-Alonso, J. Reinitz, Characterization of the Drosophila segment determination morphome, Developmental Biology 313(2): 844-862, 2008

61. L. Sanchez, D. Thieffry, A logical analysis of the Drosophila gap-gene system, J. theor. Biol. 211:115-141, 2001

62. R. Thomas, R. d'Ari. Biological feedback, Boca Raton-Florida, CRC Press, 1990

63. R. Thomas, M. Kaufman, Multistationarity, the basis of cell differentiation and memory. I. Structural conditions of multistationarity and other nontrivial behavior, Chaos 11(1): 170-179, 2001

64. S.A. Vakulenko, O. Radulescu, V.V. Gursky, K.N. Kozlov, An effective algorithm for parameter fitting in gene networks with layered pattern formation, in preparation

65. Wolfram Research, Inc., Mathematica, Version 5.0, Champaign, IL, 2003

66. L. Wolpert. Positional information and the spatial pattern of cellular differentiation, Journal of Theoretical Biology 25:1-47. 1969

67. M. Zvelebil, J. Baum. Understanding Bioinformatics, Oxford-New York, Taylor&Francis Group, 2007