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

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

Акинышш Андрей Александрович

Математическое и численное моделирование искусственных регуляторных контуров генных

сетей

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

Автореферат диссертации на соискание учёной степени кандидата физико-математических наук

28 НОЯ 2013

005540613

Барнаул — 2013

005540613

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

Научный руководитель: доктор физико-математических наук, профессор

Голубятников Владимир Петрович

Официальные оппоненты: Перцев Николай Викторович,

доктор физико-математических наук, профессор, Омский филиал Федерального государственного бюджетного учреждения науки Института математики им. С.Л. Соболева СО РАН, ведущий научный сотрудник Медведев Сергей Борисович,

доктор физико-математических наук, Федеральное государственное бюджетное учреждение науки Институт вычислительных технологий СО РАН, старший научный сотрудник

Ведущая организация: Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования «Санкт-Петербургский государственный политехнический университет»

Защита состоится 17 декабря 2013 г. в 15 часов на заседании диссертационного совета Д 003.061.02 на базе Федерального государственного бюджетного учреждения науки Института вычислительной математики и математической геофизики Сибирского отделения Российской академии наук по адресу: 630090, г. Новосибирск, пр. академика Лаврентьева, 6.

С диссертацией можно ознакомиться в библиотеке Федерального государственного бюджетного учреждения науки Института вычислительной математики и математической геофизики СО РАН.

Автореферат разослан 15 ноября 2013 года.

Ученый секретарь диссертационного совета Д 003.061.02, д.ф.-м.н.

Сорокин Сергей Борисович

Общая характеристика работы

Актуальность темы. Начало ХХ1-го века связано с резким ростом интереса к такой области науки, как системная компьютерная биология (она же постгеномная биоинформатика). За последние два десятилетия накопились огромные объёмы экспериментальных биологических данных. Увеличение вычислительных мощностей позволило вывести обработку и анализ этих данных на новый уровень. Развитие средств разработки обеспечило возможность писать высокопроизводительные программные продукты для автоматизации анализа сложных биологических систем. К сожалению, математический аппарат оказался не готов предоставить основу для подобных расчётов. Поэтому на сегодняшний день необходимо параллельно развивать соответствующее программное обеспечение и сопутствующие математические методы.

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

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

онные явления (реорганизация механизмов функционирования биологической системы).

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

Целью диссертации является изучение систем следующего вида:

±г

Х2

<

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

В некоторых случаях анализ системы (1) удаётся свести к анализу единственного уравнения с запаздывающим аргументом:

¿(г) = /(х(г - г)) - (2)

которое также само по себе используется при рассмотрении многих биологических и физических моделей.

Для достижения поставленной цели решались следующие задачи:

1. Исследовать системы (1), (2) на предмет стационарных точек.

2. Исследовать системы (1), (2) на предмет циклических траекторий.

3. Вывести условия, определяющие качественные и количественные характеристики систем (1), (2).

4. Разработать программный комплекс, способный выполнять моделирование систем (1), (2).

Основные положения, выносимые на защиту:

1. Проведён анализ нелинейного уравнения с запаздывающим аргументом (2), изучены происходящие бифуркации Андронова-Хопфа (теорема 1, лемма 1), рассмотрены стационарные точки и периодические траектории, исследованы вопросы устойчивости, выведена формула первого коэффициента Ляпунова (теорема 2).

= А(хп) - Р&Ъ

= /2(^1) - /З2Х2, = /п{%п—1)

2. Построена и изучена (леммы 3, 4, 5, 6, 7, 8) дискретная структура (граф кластеров), которая сужает область поиска циклов системы (1).

3. Проведён анализ многомерной циклической динамической системы (1), выведены условия устойчивости стационарных точек (лемма 2), рассмотрен метод 2 сведения этой системы к уравнению (2) (лемма 9).

4. Для рассматриваемых систем описано 4 метода поиска циклов: метод кластеров 1 (таблица 1), метод свёртки 2 (теорема 3), метод подобия 3 (лемма 10), метод развёртки 4 (теорема 4).

5. Разработан программный комплекс Phase Portrait Analyzer, реализующий получение качественных и количественных характеристик систем (1), (2).

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

Личный вклад автора. Результаты работы базируются на ряде фундаментальных исследований. В работе [1] В.А. Лихошвай описывает метод 2 сведения многомерной системы к уравнению с запаздывающим аргументом для поиска симметричных циклов. В настоящей работе была описана численная реализация данного метода, включённая в разработанный программный продукт. В работе [3] Тереза Фарья описывает схему вычисления первого коэффициента Ляпунова для уравнений с запаздывающим аргументом. Данная схема использована при исследовании уравнения (2). Для анализа трёхмерных и пятимерных систем вида (1) В.П. Голубятников предложил (см. [4,5]) разбить окрестность стационарной точки на 8 частей и 32 части соответственно, анализ поверхностей которых использовался при поиске периодических траекторий. В настоящей работе данные рассуждения были обобщены для многомерного случая, дискретные построения получили название «граф кластеров», был установлен ряд полезных свойств для данного графа.

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

Апробация работы. Основные результаты работы докладывались на 23 российских и международных конференциях: международная конференция «The Bioinformatics Research and Education Workshop» (2013, Германия, г. Берлин), международная конференция «Federation of European Biochemical Societies CONGRESS Mechanisms in Biology» (2013, г. Санкт-Петербург), международная

5

конференция «VI-th international conference Solitons, Collapses and Turbulence: Achievements, Developments and Perspectives» (2012, г. Новосибирск), международная конференция «The eighth international conference on bioinformatics of genome regulation and structure systems biology» (2012, г. Новосибирск), международная конференция «Математика. Компьютер. Образование» (2013, г. Пу-щино), международная научная конференция студентов, аспирантов и молодых учёных «Ломоносов» (2013, г. Москва), международная конференция «Дифференциальные уравнения. Функциональные пространства. Теория приближений» (2013, г. Новосибирск), международная научно-практическая конференция студентов и молодых учёных «Современные техника и технологии» (2011, 2012, г. Томск), международная научная студенческая конференция «Студент и научно-технический прогресс» (2010, 2011, 2012, г. Новосибирск), конференция «Дни геометрии в Новосибирске» (2011, г. Новосибирск), всероссийская научно-техническая конференция студентов, аспирантов и молодых учёных «Наука и молодёжь» (2010, 2011, 2012, 2013, г .Барнаул), международная школа-конференция «Биофизика сложных систем. Анализ и моделирование» (2013, г. Пущино), международная молодёжная школа-конференция «Теория и численные методы решения обратных и некорректных задач» (2013, г. Новосибирск), международная школа молодых учёных «Systems Biology and Bioinformatics» (2012, г. Новосибирск), международная школа-семинар «Ломоносовские чтения на Алтае» (2012, г. Барнаул), всероссийская молодёжная школа-семинар «Анализ, геометрия и топология» (2013, г. Барнаул).

На 8 семинарах: семинар профессора Р. Хофештадта (2013, Германия, г. Билефельд, Университет Билефельда), Санкт-Петербургский семинар по компьютерной биологии СПбГПУ (2013, г. Санкт-Петербург), объединённый семинар ИВМиМГ и кафедры вычислительной математики НГУ (2013, г. Новосибирск), семинар «Избранные вопросы математического анализа» ИМ СО РАН (2013, г. Новосибирск), семинар «Законы сохранения и инварианты» ИВТ СО РАН (2013, г. Новосибирск), семинар Лаборатории обратных задач математической физики ИМ СО РАН (2013, г. Новосибирск), краевой семинар по геометрии и математическому моделированию АГУ (2013, г. Барнаул), расширенный семинар кафедры прикладной математики АлтГТУ (2013, г. Барнаул).

Диссертационная работа была выполнена при поддержке РФФИ, грант 12-01-00074; стипендии Президента РФ для молодых учёных и аспирантов, СП-561-2012.5.

Публикации. Основные результаты по теме диссертации изложены в 27 печатных изданиях [7-33], 4 из которых изданы в рецензируемых журналах из списка ВАК [7-10], 8 — в прочих изданиях [11-18], 15 — в тезисах докладов [19-33].

Содержание работы

Диссертация состоит из введения, четырёх глав, заключения и двух приложений. Полный объем диссертации составляет 150 страниц с 51 рисунком и 3 таблицами. Список литературы содержит 116 наименований.

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

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

±(t) = f(x(t-r))+g(x(t)). (3)

Уравнения данного типа широко используются при анализе ряда биологических и физических моделей. Интересным феноменом в теории динамических систем является явление бифуркации — момент качественного изменения фазового портрета системы в зависимости от изменения одного или нескольких параметров (см. [3,6]). Одной из наиболее часто встречающихся бифуркаций является бифуркация Андронова-Хопфа. Она происходит, когда пара комплексно-сопряжённых собственных значений матрицы линеаризации некоторой стационарной точки проходит через мнимую ось — в этот момент рождается или исчезает некоторый предельный цикл.

Для уравнений без запаздываний бифуркация Андронова-Хопфа детально изучена. Однако, в реальных системах обычно имеют место некоторые задержки — физические и биологические процессы не могут происходить мгновенно. Теория уравнения с запаздывающим аргументом для бифуркации Андронова-Хопфа заметно сложнее классического варианта и требует отдельного рассмотрения.

В данной главе исследуется появление бифуркационных циклов Андронова-Хопфа в окрестности стационарных точек уравнений вида (3): формулируются условия существования бифуркации, выводятся бифуркационные значения, анализируются бифуркационные циклы на устойчивость. Также рассматриваются несколько полезных приёмов изучения уравнений с запаздыванием.

Для некоторой найденной стационарной точки х функции / и д считаются достаточно гладкими в её окрестности, и для их производных вводятся обозначения: fj = /(j'(i), <jj = д,3Цх), j 6 N. Формулируются и доказываются следующие утверждения:

Теорема 1. В уравнении (3) возникает серия бифуркаций Андронова-Хопфа при г е К}, где

Tt*=T*+fcr, т* =

лги

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

Т= Иш

Лемма 1. Если и д\ имеют одинаковый знак, то 2т* ^ Г* ^ 4т*, а если разный, то 2т* ^ Т*.

Устойчивость бифуркационного цикла анализируется с помощью первого коэффициента Ляпунова (см. [6]): отрицательное значение соответствует устойчивому циклу, а положительное — неустойчивому. Формулируется и доказывается следующая теорема для поиска первого коэффициента Ляпунова:

Теорема 2. Первый коэффициент Ляпунова уравнения (3) можно найти по формуле

к = ^Ке(гс2оСц +ШС21),

где

его = Фг(0)/г2о, сц = Ф1(0)/ць с21 = Ф1(0)/г2ь ^1(0) = ^^[у+^т2' где

/120 = /2е"2шТ + 32, ^11 = /2 + 92,

/121 = /2(2 ъии(-т)е-^т + №20 (-т)ешт) + 52(2№„(0) + ш20(0))+ + /зе-^ + аз,

где и)ц(0), Шц(-т), и;2о(0), и>2о(—т) находятся из:

(_е2^>20(_Г) + Ш2о(0) = ^(1 - О + - е31П,

/1Ш2о(-т) + (Э1 - 2ш)№2о(0) = С20 +_С02 ~ ^20,

-Шц(-т) + «п(0) = _ 1} _ !£П (е»г _ 1})

И] 10 /1Иц(-т) + 51№ц(0) = Си + Сц - /гп.

В конце главы формулируются некоторые подходы, которые могут быть полезны при анализе изучаемого уравнения; проводится систематизация полученных результатов (таблица 1); делается обзор различных исследований, содержащих в себе изучаемые дифференциальные уравнения. Большое значение для последующих рассуждений имеет «приведённый вид» уравнения (3) с выполненной нормализацией запаздывания:

Хг(и) = т(/{хТ{и - 1)) + з(а;г(*г)))- (4)

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

¿1= fl(x„) - 01X1, ( ¿2 = /2(2:1) - Р2Х2, ^

Хп fn{Xn—\) PnXni

Для стационарных точек данной системы выводится условие устойчивости в зависимости от дополнительного вводимого свойства — чётности точки:

Лемма 2. Устойчивость стационарной точки определяется следующим выражением:

| Т < 1, для Ф > О (чётная точка),

[Т < 1/cos(7r/n), для Ф < О (нечётная точка),

ф = т = \/[ф|.

дхп дх\ дх2 dxn-i' v

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

Пусть для некоторой стационарной точки х имеется некоторая инвариантная окрестность X, которая отвечает следующим условиям:

Условие 1. Векторное поле, определяемое системой (5), по всей поверхности X направлено к стационарной точке, т.е. траектории могут только входить в инвариантную окрестность I и не могут из неё выходить.

Условие 2. Окрестность X не содержит точек экстремума функций fi, т.е. для любой точки х € Т по любой координате хг выполняется ¡¡J' i (x¡-\) ф 0.

Назовём кластером QSí;2...£„ (с — двоичный n-мерный вектор, £г G {0,1}) часть окрестности X в которой: хг ^ í¿, если с, — 0; хг > x¿, если £г = 1. Будем говорить, что кластеры граничат по i-ой грани, если их бинарные векторы различаются только в г-ой компоненте. Показывается, что внутри отдельно взятого кластера цикл появится не может, а на их границе между двумя кластерами векторное поле во всех точках направлено в одну и ту же сторону. Составим граф, вершины которого соответствуют кластерам, а ребра — переходам вдоль

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

9

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

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

Лемма 3. Переходы между кластерами описываются следующими схемами. Для^<0: 00-»01, 11-» 10. Для 0: 01 -» 00, 10 -» 11.

Здесь кластеры обозначаются парой бинарных индексов £ — 1 и Ь, стрелки показывают направление соответствующего ребра.

Лемма 4. При переходе из одного кластера в другой потенциал либо не меняется, либо уменьшается на 2.

Лемма 5. В рамках одной системы потенциалы всех кластеров обладают одинаковой чётностью.

Лемма 6. Чётность потенциалов кластеров системы совпадает с чётностью стационарной точки.

Лемма 7. В пределах фиксированной размерности пространства п графы кластеров стационарных точек одной чётности изоморфны друг другу.

Лемма 8. Количество вершин на р-ам потенциальном уровне определяется выражением 2СЦ.

Граф кластеров примечателен тем, что значительно сужает область поиска циклов. Заметим, что каждый существующий цикл принадлежит строго одному потенциальному уровню, так как по лемме 4 траектория может либо оставаться в текущем потенциальном уровне, либо опускаться на более низкие. Поэтому имеет смысл говорить о существовании циклов на том или ином уровне; при нахождении нового цикла для его качественного описания нам необходимо знать потенциальный уровень, на котором он находится. В работе [5] проводится подробное доказательство существования цикла на первом потенциальном уровне для нечётной гиперболической стационарной точки через теорему Брауэра. На основе этих рассуждений сформулируем метод поиска циклов в рассматриваемых динамических системах:

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

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

Рис. 1: Пример графа кластеров для 3-мерной системы

Также рассматривается симметричная система вида:

= ^(Ц, Хп, Хп-1, ..., х2), X2 = XI, Хп, ..., Ж3), ^^

Хп Р(хп, Хп_\, ■ ■ ■ > ^1)-

Формулируются и доказываются следующие утверждения:

Лемма 9. Для симметричного цикла системы (б) с периодом Т выполняется условие сдвига:

^Т) =!*_!(*), г = Т~гг, ре {0,1,...,п-1}. (7)

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

Таким образом, задача поиска симметричных циклов системы (6) сводится к поиску циклов уравнения с запаздыванием фиксированного периода. Сформулируем соответствующий метод (основа метода взята из [1]):

Метод поиска циклов 2 (Метод свёртки). Если для симметричной системы (б) удалось построить график зависимости соответствующего приведённого запаздывающего уравнения (4) от т, то условие Тг = ^ будет определять всё множество силшетричных циклов системы.

В главе 3 полученные результаты применяются к системам хилловского типа, которые активно используются при аппроксимации биологических процессов. В систему (5) подставляется функция Хилла /:

Хг = /(а*-!) - ¡Зхи /(*) = (8)

1 + х7

Рис. 2: Зависимость Тт от параметров а, /3, г

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

Формулируются и доказываются дополнительные свойства систем хил-ловского типа, приводятся два дополнительных метода поиска циклов (системы задаются четвёрткой параметров (п, а, /3,7) — размерностью пространства и параметрами из (8)):

Лемма 10. Множества симметричных циклов систем {п,а\, ¡З^^) и (п, 02,^2,7) для которых выполняется условие

Т = Т V

/?1 /32

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

Метод поиска циклов 3 (Метод подобия). Если для системы (п, а, /3,7) найден цикл, то для системы (п, а/с, ¡Зк, 7} (к е К+/) также можно получить соответствующий цикл по лемме 10.

Теорема 4. Симметричные циклы системы (п, а, /3,7) отображаются на циклы системы (пк, а, /3,7) (к 6 причём каждый цикл второй системы можно

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

12

Метод поиска циклов 4 (Метод развёртки). Если для системы (п, а, ¡5, у) найден цикл, то для системы (пк, а, ¡3,7) (к £ N) можно получить соответствующий цикл по теореме 4. Данную операцию будем называть развёртыванием цикла.

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

В главе 4 описывается программный комплекс Phase Portrait Analyzer, который был разработан в рамках настоящей работы для численного моделирования изучаемых процессов. Визуальная часть разработана на платформе .NET, она включает в себя гибкий интерфейс для формирования целевых систем, двумерные и трёхмерные интерактивные визуализации, графики различных параметров, форматированные отчёты, диспетчеризацию пользовательских команд. Математическое ядро программного продукта разработано на языке R, который представляет собой открытую кроссплатформенную среду для математических вычислений с открытым программным кодом. Для численного решения систем дифференциальных уравнений используется алгоритм Isoda, разработанный в Ливерморской национальной лаборатории, который способен автоматически адаптировать параметры алгоритма в ходе моделирования и переключаться между жёсткими и нежёсткими задачами.

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

В приложении А приводятся примеры изучаемых систем хилловского типа с подробным описанием качественных и количественных характеристик, включая расчётные значения, графики параметров и проекции найденных циклов на различные плоскости (пример приведён на рисунке 3).

Рис. 3: Конфигурации циклов 105-мерной системы, проекция на плоскость собственных

векторов

В приложении В приводятся листинги кода на языке Я для ряда проводимых численных экспериментов.

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

1. Лихошвай В.А., Голубятников В.П., Демиденко Г.В., Фадеев С.И., Евдокимов A.A. Теория генных сетей // Системная компьютерная биология. — Интеграционные проекты СО РАН / Рос. акад. наук, Сиб. отд-ние. Новосибирск: СО РАН, 2008. - С. 397^182.

2. Лихошвай В.А., Матушкин Ю.Г., Ратушный A.B., Ананько Е.А., Игнатьева Е.В., Подколодная О.В. Обобщенный химико-кинетический метод моделирования генных сетей // Молекулярная биология.

- 2001. - Т. 35, № 6. - С. 1072-1079.

3. Hale J., Magalhaes L.T., Oliva W. Dynamics in infinite dimensions H Applied Mathematical Sciences. — 2002. - Vol. 47.

4. Гайдов Ю.А., Голубятников В.П. О некоторых нелинейных динамических системах, моделирующих несимметричные генные сети // Вестник ИГУ. Серия «Математика, механика, информатика». - 2007. - Т. 7, № 2. - С. 8-17.

5. Голубятников В.П., Голубятников И.В., Лихошвай В.А. О существовании и устойчивости циклов в пятимерных моделях генных сетей // Сибирский журнал вычислительной математики. — 2010. — Т. 13, №4,- С. 403-411.

6. Kuznetsov Y. A. Elements of Applied Bifurcation Theory // Applied Mathematical Sciences. — 1998. — Vol. 112.- 591 pp.

Статьи автора в журналах из списка ВАК

7. Акиньшин A.A. Бифуркация Андронова-Хопфа для некоторых нелинейных уравнений с запаздыванием // Сиб. журн. индустр. матем.. — 2013. — Т. XVI, № 3 (55). — С. 3-15.

8. Акиньшин A.A. Изучение дискретных структур в некоторых циклических динамических системах И Ползуновский вестник. — 2012. — № 4. — С. 214-218.

9. Акиньшин A.A., Голубятников В.П., Голубятников И.В. О некоторых многомерных моделях функционирования генных сетей И Сиб. журн. индустр. матем.. — 2013. — Т. XVI, № 1 (53). — С. 3-9.

10. Акиньшин A.A., Голубятников В.П. Геометрические характеристики циклов в некоторых симметричных динамических системах // Вестник НГУ. Серия «Математика, механика, информатика».

- 2012. - Т. 12, № 2. - С. 3-12.

В прочих изданиях ^

11. Akinshin A.A., Golubyatnikov VP., Likhoshvai V.P. Mathematical and computational models of gene networks functioning. — The 11 th Bio informatics Research and Education Workshop. — Berlin, Germany.

- 2013. — 5 pp.

12. Акиньшин A.A. Поиск периодических траекторий в математических моделях генных сетей // Труды семинара по геаиетрии и математическому моделированию. — 2013. — С. 4-9.

13. Акиньшин А.А., Голубятников В.П. Вопросы единственности циклов у нелинейных динамических систем специального вида // Труды международной конференции «Дни геометрии в Новосибирске

- 2011». — 2012. — С. 7-15.

14. Акиньшин А.А. Применение численных методов для решения задач биоинформатики // Сборник научных статей шкалы-семинара «Ломоносовские чтения на Алтае». — 2012. — Т. 2. — С. 11-16.

15. Акиньшин А.А. Пример использования языка R для решения задач биоинформатики // Журнал «Горизонты образования». 10-ая Конференция «Наука и молодёжь». — 2013. — С. 6-8.

16. Акиньшин А.А., Голубятников В.П. Математическое и численное описание фазовых портретов некоторых нелинейных динамических систем // Журнал «Горизонты образования». 9-ая Всероссийская научно-техническая конференция «Наука и молодёжь». — 2012. — С. 6-9.

17. Акиньшин А.А., Голубятников В.П. Некоторые математические и вычислительные задачи биоинформатики // «Горизонты образования». 8-ая конференция «Наука и молодёжь». — 2011. — С. 6-9.

18. Акиньшин А.А., Голубятников В.П. Трёхмерные модели генных сетей // Журнал «Горизонты образования». 7-ая конференция «Наука и молодёжь». — 2010. — С. 8-10.

В тезисах докладов

19. Акиньшин А.А. Численный анализ некоторых моделей биоинформатики // ХХ-ая конференция серии «Математика. Компьютер. Образование». Тезисы. — Москва: 2013. — С. 203.

20. Акиньшин А.А. Численный анализ периодических режимов генных сетей // Материалы Международного молодежного научного форума «Ломоносов-2013» . — Москва. — 2013. — 1 с.

21. Akinshin А.А. Analysis of phase portraits in some gene networks models // 5th International Young Scientists School "System Biology and Bio informatics". — Novosibirsk, Institute of Cytology and Genetics.

- 2013,- P. 27.

22. Akinshin A.A. Numerical analysis of gene networks models // 8th FEBS Congress, Saint Petersburg, Russia, July 6-11, 2013. - Vol. 280. - 2013. - P. 547. - SW06.W29-32 - DOI: 10.11U/febs.l2396.

23. Akinshin A.A., Golubyatnikov V.P. Oscillating trajectories in some nonlinear dynamical systems // International Conference «Differential Equations. Function Spaces. Approximation Theory» / 2013. — P. 311.

24. Акиньшин A.A., Голубятников В.П. Компьютерные и математические модели функционирования генных сетей // 50-ая конференция «Студент и научно-технический прогресс» / 2012. — С. 308.

25. Акиньшин А.А., Гачубятников В.П. Математическое и компьютерное моделирование периодических режимов генных сетей // XVIII Международная научно-практическая конференция студентов и молодых учёных «Современные техника и технологии». / Т. 2. — 2012. — С. 255-256.

26. Akinshin А.А., Golubyatnikov V.P. Non-uniqueness of cycles in gene networks models // International conference on bioinformatics of genome regulation and structure \ systems biology. — Novosibirsk, Institute of Cytology and Genetics. — 2012. — P. 28.

27. Akinshin A.A., Golubyatnikov V.P., Gaidov Yu.A., Golubyatnikov I. V. Unstable cycles in gene networks models H intern, conf. on bioinformatics of genome regulation and structure. — Novosibirsk, Institute of Cytology and Genetics. - 2012. - P. 29.

28. Akinshin A.A. Computer analysis of phase portraits in gene networks models // Abstracts of Young scientist's school "Bioinformatics ans systems biology". — 2012. — P. 13.

29. Akinshin A.A., Golubyatnikov V.P. On Nonuniqueness of Cycles in Dissipative Dynamical Systems of Chemical Kinetics // Vl-th intern, conf. Solitons, Collapses and Turbulence. — 2012. — Pp. 71-72.

30. Акиньшин A.A., Голубятников В.П. Математическое и компьютерное моделирование периодических режимов функционирования генных сетей // XVII Международная научно-практическая конференция «Современные техника и технологии». / Т. 2. — 2011. — С. 283-284.

31. Акиньшин А.А., Голубятников В.П. Компьютерные и математические модели генных сетей // 49-ая Международная научная студенческая конференция «Студент и научно-технический прогресс». Секция «Математика» / НГУ. - Новосибирск: РИЦ НГУ, 2011. - С. 311.

32. Akinshin А.А., Golubyatnikov V.P., Golubyatnikov I.V. Mathematical and numerical modeling of gene network functioning // International Conference "Modern Problems of Mathematics, Informatics and Bioinformatics". — Novosibirsk: 2011. — P. 81.

33. Акиньшин A.A., Голубятников В.П. Трёхмерные модели генных сетей // XLVIII Конференция «Студент и научно-технический прогресс». Секция «Математика» / 2010. — С. 309-310.

Подписано в печать 13.11.2013 г. Печать цифровая. Бумага офсетная. Формат 60x84/16. Усл. печ. л. 2 Тираж 150 экз. Заказ № 189

Отпечатано в типографии «Срочная полиграфия» ИП Малыгин Алексей Михайлович 630090, Новосибирск, пр-т Академика Лаврентьева, 6/1, оф. 104 Тел. (383) 217-43-46, 8-913-922-19-07

Алтайский Государственный Технический Университет

им. И.И. Ползунова

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

04201 451 444 Акиньшин Андрей Александрович

Математическое и численное моделирование искусственных регуляторных контуров генных сетей

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

программ»

Диссертация на соискание учёной степени кандидата физико-математических наук

Научный руководитель: д.ф.-м.н, профессор Голубятников В.П.

Барнаул - 2013

Содержание

Введение..................................................................5

1 Дифференциальные уравнения с запаздыванием................19

1.1 Постановка задачи..................................................19

1.2 Бифуркация Андронова-Хопфа....................................20

1.3 Первый коэффициент Ляпунова ..................................24

1.4 Общие подходы к анализу..........................................30

1.5 Основные результаты ..............................................33

1.6 Обзор моделей процессов с запаздыванием......................34

2 Многомерные циклические системы...............38

2.1 Постановка задачи..................................................38

2.2 Стационарные точки................................................39

2.3 Граф кластеров......................................................42

2.3.1 Построение графа......................42

2.3.2 Потенциалы вершин........................................44

2.3.3 Пример графа........................47

2.3.4 Значение графа кластеров....... ...........48

2.4 Симметричные системы............................................51

3 Системы Хилловского типа.....................54

3.1 Постановка задачи..................................................54

3.2 Стационарные точки................................................56

3.3 Алгоритм поиска симметричных циклов..........................58

3.4 Подобие систем Хилловского типа................60

3.5 Биологическая интерпретация................... 64

4 Компьютерное моделирование...................68

4.1 Численный анализ..................................................68

4.1.1 Выбор математического пакета............................68

4.1.2 Моделирование динамических систем....................73

4.2 Функциональное описание программного комплекса......76

4.3 Техническое описание программного комплекса ........80

4.3.1 Модульная структура ......................................80

4.3.2 Формат хранения данных..................................81

4.3.3 Программная архитектура..................................82

Заключение ..............................................................85

Список рисунков ............................87

Список таблиц..............................90

Литература................................91

А Примеры циклических систем Хилловского типа........108

А.1 3-мерная система..........................110

А.2 4-мерная система..........................113

А.З 5-мерная система..........................116

А.4 6-мерная система..........................118

А.5 7-мерная система..........................120

А.6 8-мерная система..........................124

А.7 9-мерная система..........................128

А.8 12-мерная система.........................130

А.9 18-мерная система.........................134

АЛО 105-мерная система.........................137

В Листинги кода на К.........................145

Введение

Актуальность

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

Одним из наиболее актуальных направлений системной биологии в последние годы стала теория генных сетей (см. [1]). Генная сеть — это комплекс координировано функционирующих генов, обеспечивающих формирование различных фенотипических признаков организмов (молекулярных, биохимических, физиологических, морфологических, поведенческих и т.д.). К основным составляющим генной сети относятся ДНК, РНК и белки. Реальные генные сети являются очень сложными объектами, включающими в себя тысячи различных элементов, проводить их полноценный анализ весьма затруднительно. В определённых ситуациях при условии по-

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

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

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

Основные исследователи

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

Одними из первых исследователей, которые занялись изучением интересных для нас систем, являются Леон Гласс (Leon Glass1) и Михаил Маки (Michael С. Mackey2) из университет Макгилла. Одной из их самых основательных работ является монография [3], в которой подобно рассматриваются определённые классы динамических систем с их физическими и биологическими приложениями, особое внимание уделяется хаотическим явлениям. В 1977 году Гласс и Маки выписали одно из наиболее простых биологически значимых уравнений с запаздыванием (уравнение Гласса-Маки), генерирующее хаос (см. [4], [5]). В работах [6,7] Гласс рассматривает в окрестности стационарных точек дискретные построения подобные графу кластеров, построенному в настоящей работе. В работах [8,9] Маки описывает некоторые интересные особенности функций Хилла, которая также анализируется в настоящей работе.

Значительный вклад в развитие математической биологии внёс Джеймс Мюррей (James D. Murray3). В его работах [10, 11] выполняется подробный анализ разнообразных математических моделей с биологическим приложением. Также обсуждается важность и значимость поднимаемых проблем, некоторые из которых рассматриваются в настоящей работе.

'http ://www.medicine.mcgill.ca/physio/glasslab/

2http://www.medicine.mcgill.ca/physio/mackeylab/

3http://depts.Washington.edu/amath/people/faculty/murray/

Классические результаты по динамическим системам (в том числе с запаздывающим аргументом) можно найти в работах [12-17] Джона Маллет-Парета (John Mallet-Paret4) и Роджера Ньюсбаума (Roger D. Nussbaum5).

Кроме прочего, в работе используется метод анализа устойчивости бифуркационных циклов в системах запаздывающего типа с помощью первого коэффициента Ляпунова. Впервые основательно подобный метод описала Тереза Фарья (Teresa Faria) в работе [18], а в дальнейшем её результаты в более развёрнутом виде были представлены в главе 8.2 работы [19]. Конкретные вычисления по соответствующей схеме на примере модели лейкемии были рассмотрены в статьях [20-22] Анкой Вероникой Ион (Anca Veronica Ion6).

Ряд специфических особенностей уравнений с запаздывающим аргументом также рассматривался группой учёных из Гергели Роста, Тибора Кристина и Эдуардо Лица (Gergely Röst7, Tibor Krisztin8, Eduardo Liz9). Некоторые их идеи (см. [23-30]) успешно применяются для анализа изучаемых в рамках настоящей работы систем.

Некоторые полезные свойства уравнений с запаздыванием для описания генных сетей рассматривает Геннадий Владимирович Демиденко10 в работах [31-36].

Непосредственное влияние на данную работу оказали результаты Ли-хошвая Виталия Александровича11, который подробно описал (см. [{]) изучаемые ниже биологические системы. Некоторые методы анализа рассмат-

4http://research.brown.edu/research/profile.php?id=10399

5http://www.math.rutgers.edu/-nussbaum/

chttp://www.ima.ro/people/averion.html

7http://www.math.u-szeged.hu/~rost/

8http://www.math.u-szeged.hu/-krisztin/

9http://scholar.google.com/citations?user=gqDmmBcAAAAJ&hl=en

I0http://www.math.nsc.ru/LBRT/d5/ruswin/demidenko.htm

"http://www.bionet.nsc.ru/kib/?page_id=589

риваемых моделей базируются именно на его результатах. Важные характеристики целевых моделей также описываются в работах [37,38].

Цели работы

Основной целью работы является изучение систем следующего вида:

/

¿1 ¿2

<

•Еп

V

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

В некоторых случаях анализ системы (1) удаётся свести к анализу единственного уравнения с запаздывающим аргументом:

х(Ь) = /(х(г-т))-Рх®, (2)

которое также само по себе используется при рассмотрении многих биологических и физических моделей.

Для достижения поставленной цели необходимо было решить следующие задачи:

1. Исследовать системы (1), (2) на предмет стационарных точек.

2. Исследовать системы (1), (2) на предмет циклических траекторий.

3. Вывести условия, определяющие качественные и количественные характеристики систем (1), (2).

4. Разработать программный комплекс, способный выполнять моделирование систем (1), (2).

= Л Оп) -= /2(^1) ~ /02^2,

= /тгО^П—1)

Основные результаты

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

1. Проведён анализ нелинейного уравнения с запаздывающим аргументом (2), изучены происходящие бифуркации Андронова-Хопфа (теорема 1, лемма 1), рассмотрены стационарные точки и периодические траектории, исследованы вопросы устойчивости, выведена формула первого коэффициента Ляпунова (теорема 2).

2. Построена и изучена (леммы 3, 4, 5, 6, 7, 8) дискретная структура под названием граф кластеров, которая сужает область поиска циклов системы (1).

3. Проведён подробный анализ многомерной циклической динамической системы (1), выведены условия устойчивости стационарных точек (лемма 2), рассмотрен метод 2 сведения этой системы к уравнению (2) через условие сдвига (лемма 9).

4. Для рассматриваемых систем описано 4 метода поиска циклов: метод кластеров 1 (таблица 2), метод свёртки 2 (теорема 3), метод подобия 3 (лемма 10), метод развёртки 4 (теорема 4).

5. Разработан программный комплекс Phase Portrait Analyzer, автоматизирующий получение качественных и количественных характеристик систем (1), (2). Исследуемые объекты моделируются различными численными методами, результаты которых можно изучать с помощью двумерных и трёхмерных интерактивных визуализаций и в виде форматированных отчётов.

Научная новизна и практическая значимость

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

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

Результаты могут быть полезны при анализе разнообразных биологических и физических моделей.

Апробация работы

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

23 конференции

• Международная конференция «The Bioinformatics Research and Education Workshop»12, 2013, Германия, г. Берлин

• Международная конференция «Federation of European Biochemical Societies CONGRESS Mechanisms in Biology»13, 2013, г. Санкт-Петербург

• Международная конференция «VI-th international conference Solitons, Collapses and Turbulence: Achievements, Developments and Perspectives»14, 2012, г. Новосибирск

12http://cmb.molgen.mpg.de/brew2013/

13http://www.febs-2013.org/eng/default.aspx

14http://conf.nsc.ru/sct2012/en

• Международная конференция «The eighth international conference on bioinformatics of genome regulation and structure systems biology»15, 2012, г. Новосибирск

• Международная конференция «Математика. Компьютер. Образование»16, 2013, г. Пущино

• Международная научная конференция студентов, аспирантов и молодых учёных «Ломоносов»17, 2013, г. Москва

• Международная конференция «Дифференциальные уравнения. Функциональные пространства. Теория приближений»18, 2013, г. Новосибирск

• Международная научно-практическая конференция студентов и молодых учёных «Современные техника и технологии»19, 2011, 2012, г. Томск

• Международная научная студенческая конференция «Студент и научно-технический прогресс»20, 2010, 2011, 2012, г. Новосибирск

• Конференция «Дни геометрии в Новосибирске»21, 2011, г. Новосибирск

• Всероссийская научно-техническая конференция студентов, аспирантов и молодых учёных «Наука и молодёжь», 2010, 2011, 2012, 2013, г .Барнаул

15http: //соnf . rise. ru/BGRSSB2012

16http: / /www. nice. su/

,7http://lomonosov-msu.ru/rus/lomonosov.html

18http://math.nsc.ru/conference/sobolev/105/

l9http://portal.tpu.ru/science/копf/ctt

20http : //issc . nsu . ru/

21http://math.nsc.ru/conference/geomtop2011/

• Международная школа-конференция «Биофизика сложных систем. Анализ и моделирование»22, 2013, г. Пущино

• Международная молодёжная школа-конференция «Теория и численные методы решения обратных и некорректных задач»23, 2013, г. Новосибирск

• Международная школа молодых учёных «Systems Biology and Bioinformatics»24 25, 2012, 2013, г. Новосибирск

• Международная школа-семинар «Ломоносовские чтения на Алтае»26, 2012, г. Барнаул

• Всероссийская молодёжная школа-семинар «Анализ, геометрия и топология»2', 2013, г. Барнаул

8 семинаров

• Семинар профессора Р. Хофештадта28, 2013, Германия, г. Билефельд, Университет Билефельда (Universität Bielefeld)

• Санкт-Петербургский семинар по компьютерной биологии, 2013, г. Санкт-Петербург, Санкт-Петербургский государственный политехнический университет

• Объединённый семинар ИВМиМГ и кафедры вычислительной математики НГУ, 2013, г. Новосибирск, Институт вычислительной математики и математической геофизики СО РАН

22http://www.mce.su/biophys/

23http://conf.nsc.ru/tcmiip2013/ru/general_info

24http://conf.nsc.ru/SBB2013/

25http://conf.nsc.ru/BGRSSB2012/bgrssb2012_young_scientist_school

26https://sites.google.com/site/lomchten/home

27https://sites.google.com/site/geometryaltai/home

28http://www.techfak.uni-bielefeld.de/ags/bi/staff/headofdepartment/index.shtml

• Семинар «Избранные вопросы математического анализа»29, 2013, г. Новосибирск, Институт математики СО РАН

• Семинар «Законы сохранения и инварианты»30, 2013, г. Новосибирск, Институт вычислительных технологий СО РАН

• Семинар лаборатории обратных задач математической физики31, 2013, г. Новосибирск, Институт Математики СО РАН

• Краевой семинар по геометрии и математическому моделированию32, 2013, г. Барнаул, Алтайский Государственный Университет

• Расширенный семинар кафедры прикладной математики АлтГТУ33, 2013, г. Барнаул, Алтайский Государственный Технический Университет им. И.И. Ползунова

Публикации

Основные результаты по теме диссертации изложены в 27 печатных изданиях [39-65], 4 из которых изданы в рецензируемых журналах из списка ВАК [53,50,57,59], 8 — в прочих изданиях [39,42,47,51,54,60-62], 15 — в тезисах докладов [40,41,43-46,48-50,52,55,58,63-65]. Всего автору принадлежит 33 публикаций [39-7]].

29Ь1^р : //таЬЬ. пзс. ги/seITlinar/пlatarl/2013 . Ы;т1

301п^рз : //зл^е5 . доод1е . сот/з^е/2акапузоЬгапеп1а11гуу,аг1апЪу/

31111^р: //та 1:11 . пэс. ги/зетз.паг/а11. ЬШ!

32Ы:1:р: //№№». ип1-аД^а1. ги/1£то/кдтте/деот_Бет:1паг/

33Ь^р: //ют*. аИ;з1:и. ги/з1:гис1:иге/сЬа1г/рт/

Личный вклад автора

Результаты работы базируются на следующих исследованиях:

• В работе [1] В.А. Лихошвай описывает метод 2 сведения многомерной системы к уравнению с запаздывающим аргументом для поиска симметричных циклов. В настоящей работе была описана численная реализация данного метода, которая была включена в разработанный программный продукт.

• В работе [19] Тереза Фарья описывает схему вычисления первого коэффициента Ляпунова для уравнений с запаздывающим аргументом. Данная схема использована при исследовании уравнения (2).

• Для анализа трёхмерных и пятимерных систем вида (1) В.П. Голубятников предложил (см. [72,73]) разбить окрестность стационарной точки на 8 частей и 32 части соответственно, анализ поверхностей которых использовался при поиске периодических траекторий. В настоящей работе данные рассуждения были обобщены для многомерного случая, дискретные пос