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

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

/¿/{<чо1Р<с4

Макаров Сергей Сергеевич

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

ДИНАМИЧЕСКИХ ПРОЦЕССОВ В ПРОСТРАНСТВЕННО-НЕОДНОРОДНЫХ БИОЛОГИЧЕСКИХ СИСТЕМАХ

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

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

3 О МАЙ 2013

Москва-2013

005060535

005060535

Работа выполнена на кафедре компьютерных методов физики Физического факультета Московского государственного университета имени М. В. Ломоносова.

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

Грачёв Евгений Александрович

Официальные оппоненты: Полежаев Андрей Александрович,

доктор физико-математических наук, заведующий сектором теоретических проблем биофизики Отделения теоретической физики Физического института имени П. Н. Лебедева Российской академии наук

Тихонов Николай Андреевич, доктор физико-математических наук, профессор кафедры математики физического факультета Московского государственного университета имени М. В. Ломоносова

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

учреждение науки Институт математических проблем биологии Российской академии наук, г. Пущино

Защита диссертации состоится 14 июня 2013 года в 16 часов на заседании диссертационного совета Д 501.002.09 при Московском государственном университете имени М.ВЛомоносова по адресу: 119991, Москва, Ленинские горы, МГУ^ НИВЦ, Большой конференц-зал.

С диссертацией можно ознакомиться в Научной библиотеке МГУ имени М. В. Ломоносова НИВЦ МГУ (Ломоносовский проспект, д. 27).

Автореферат разослан «30»£Ц«^<и5 2013 года. • ; .

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

Суворов В. В.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

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

Первый подход заключается в усреднении переменных по пространству и рассмотрению системы как гомогенной. Такие модели называют кинетическими. Математически они представляют собой совокупность систем обыкновенных дифференциальных уравнений относительно усреднённых значений исследуемых переменных. В частности, в недавних работах М.-А. Dronne, G. Chapuisat, Е. Grenier была представлена кинетическая модель поведения клеток головного мозга в условиях инсульта. Модель хорошо описывала изменение состояния клеток в эпицентре зоны поражения, в условиях сниженного кровотока. Однако медицинский интерес не ограничивается этими клетками. Важным является и описание пространственного развития поражения на соседние области, в том числе и на те, кровоток в которых не нарушен. Зачастую необходимо определить пространственные характеристики зоны поражения, и в этом случае кинетические модели являются недостаточно информативными.

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

з

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

Второй подход применим непосредственно для гетерогенных моделей. В нем динамика переменных системы описывается совокупностью уравнений в частных производных с граничными и начальными условиями. Такие уравнения, как правило, решаются численно с помощью методов, основанных на конечно-разностных аппроксимациях уравнений в частных производных. Так, в работе & СЬаршэа! рассматривалась конечно-разностная модель развития инсульта Однако эта модель имеет проблемы с заданием граничных условий, и кора головного мозга вынужденно рассматривается как область простой геометрии (квадрат). Существуют конечно-разностные модели, описывающие электронный транспорт на биологической мембране. Одним из самых значимых успехов такого моделирования являются работы А. В. Вершубского, А. И. Приклонского, А. Н. Тихонова. Вместе с тем, общими недостатками конечно-разностных моделей являются их общая трудоёмкость и сложность задания граничных условий, дополненная также тем обстоятельством, что исследуемый объект в природе может иметь различную форму (как, например, клетка или её органеллы).

Третий подход основан на дискретизации системы и переходе от переменных, описывающих концентрации, к отдельным дискретным элементам. Методы, используемые при таком моделировании, включают в себя методы Монте-Карло, методы молекулярной динамики, методы частиц и т.д. На настоящий момент не существует дискретных моделей развития ишемического инсульта. Дискретные стохастические модели, описывающие динамику отдельных переносчиков электронов в липидной мембране, появились относительно недавно и активно развивались Г. Ю. Ризниченко, А. Б. Рубиным, Е. А Грачёвым и другими. Достоинством дискретных моделей является удобство их реализации в виде параллельных программ, однако большое общее количество элементов, тем не менее, может сделать трудоёмкими и такие модели.

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

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

Предлагаемые модели сочетают в себе лучшие черты всех перечисленных выше подходов и в значительно степени свободны от

4

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

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

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

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

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

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

- модели электронного транспорта между белковыми комплексами в липидной мембране тилакоида во время фотосинтеза.

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

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

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

5

проведенных с использованием этих комплексов, доказана адекватность предложенных моделей.

Соответствие диссертации паспорту научной специальности. Работа соответствует паспорту специальности 05.13.18 - «Математическое моделирование, численные методы и комплексы программ», так как включает в себя разработку, исследование и редукцию новых математических моделей биологических систем, а также разработку численных методов решения сформулированных математических задач и создание новых комплексов программ для нахождения численного решения.

Основные результаты работы, выносимые на защиту

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

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

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

4) Получены результаты расчетов созданных программ, по существующим экспериментальным данным доказана адекватность математических моделей.

Научно-практическая значимость работы. Предложенные в диссертации математические модели биологических систем могут быть использованы для решения прикладных задач биологии и медицины (диагностика ишемического инсульта и др.).

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

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

Личное участие автора в выполнении работы. Постановка биофизической задачи о моделировании поведения группы нейронов и астроцитов в условиях ишемического инсульта принадлежит автору совместно с научным руководителем, а также докторами медицинских наук Кочетовым А. Г. и Губским Л. В. Формулировка математической модели проведена автором лично, с консультациями Грачёва Е. А.

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

6

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

Апробация результатов. Материалы диссертации докладывались на международной конференции «Математика. Компьютер. Образование» в Объединенном Институте Ядерных Исследований, г. Дубна, январь 2010 года и январь 2012 года; на международной конференции студентов, аспирантов и молодых учёных по фундаментальным наукам "Ломоносов-2011" на физическом факультете МГУ, апрель 2011 года; на рабочем семинаре сектора информатики и биофизики сложных систем кафедры биофизики биологического факультета МГУ под руководством д.ф.-м.н. Г. Ю. Ризниченко, декабрь 2011 года и январь 2012 года; на научном семинаре кафедры высшей нервной деятельности биологического факультета МГУ, март 2012 года; на научном семинаре кафедры биофизики физического факультета МГУ, апрель 2012 года; на научном семинаре кафедры вычислительной математики Московского Физико-Технического Института (МФТИ), июнь 2012 года; на научном семинаре "Обратные задачи математической физики" под руководством д.ф.-м.н. А. Г. Яголы, д.ф.-м.н. А. Б. Бакушинского и д.ф.-м.н. А. В. Тихонравова, проводящемся в НИВЦ МГУ, октябрь 2012 года; на научном семинаре кафедры компьютерных методов физики физического факультета МГУ, ноябрь 2012 года; на научно-методологическом семинаре НИВЦ МГУ под руководством д.ф.-м.н. А. В. Тихонравова, декабрь 2012 г. на научном семинаре института математических проблем биологии РАН под руководством д.ф.-м.н. М. Н. У станина, декабрь 2012 года; на научном семинаре в государственной клинической больнице №31 под руководством д.м.н. Л. В. Губского, март 2013 года

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

Структура п объем диссертации. Диссертация состоит из титульного листа, оглавления, введения, четырех глав, заключения, двух приложений и списка литературы, содержащего 93 наименования. Общий объем диссертации составляет 117 страниц, в том числе 29 рисунков и 10 таблиц.

ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ

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

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

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

Рассматривается система из N нейронов, М астроцитов и межклеточного пространства. Система схематически изображена на Рис. 1 (для наглядности на Рис. 1 N = 3 ,М = 3 ).

(Й) нейроны

ястроцнты

межклеточное пространство

синпптнческяя связь

щелевые контакты

ионные к*йн;|.'1м

Рис. 1. Схематическое представление модели

Вершины nu...,nN соответствуют нейронам, вершины а1,....ам -астроцитам, а область Е изображает межклеточное пространство. Положениям клеток в межклеточном пространстве соответствуют точки ri >•••>»V+M- В разных областях и вершинах данной совокупности объектов зададим следующие переменные:

- и,.(г,г), геЕ, г' = 1,...,5, определяющие концентрации пяти типов ионов (соответственно, Са2+, Na+, К+, СГ, glu) в межклеточном пространстве;

- = j = l...,N, и w.bj.t),/=1,...,5, j = \,...,M, определяющие внутриклеточные концентрации пяти типов ионов для N нейронов и M астроцитов соответственно;

v{p},t), j =l,...,N, и v(aj,t), j = l,...,M, определяющие электрический потенциал на мембранах нейронов и астроцитов;

- ?.(°j>f)> ' = 1> -;3, 7 = 1,...,М, определяющие ряд биофизических величин, влияющих на изменение концентрации ионов Са2+ в астроцитах.

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

8

л, - г,,...,дм — Глг+м соответствует модель транспорта ионов через расположенные на мембране ионные каналы. Отрезки, соединяющие вершины, обозначенные буквой п, отражают синаптическую передачу мембранного потенциала между нейронами. Отрезки, соединяющие вершины с буквой а, отражают транспорт молекул инозитол-1,4,5-трифосфата (1Р3) между астроцитами, стимулирующий увеличение концентрации ионов кальция в этих клетках. Наконец, в центре схемы располагается область £, которая представляет собой трёхмерную область и предполагается бесконечной ввиду того, что рассматриваемая система составляет лишь малую часть от всего головного мозга.

Разработана математическая постановка задачи, соответствующая процессам, протекающим в биологической системе, и представляющая собой следующую систему уравнений:

(г,/) = О, Д и, (г,/) + g¡ (и,,..., ы5, -и- 5, V, г) / = 1,...,5,

от

"/(«■> 01 г-н. =Иом / = 1,...Д

А'+Л/ ^ (г)

ы5,IV,,...,К,г) = £ /у (и, те, —, 1 = 1,...,5,

¡Л V

у; (и,,...,и,, »с,,..., и>5, К) = /я, (и, (г )...,и5(Гу)V,{л,)...,(и,), )) +

+ 2>Д И"*)) Л», (и, (г,),..., И, (гу } V, (и,(И/ ), У(пк к

I = 1,... ,5, У = 1,. ..,ЛГ,

/.,■ (и,,..., «5, м-,,..., И',,К) = /а,, (и, (г,),..., И5 (г,), (ау_„),..., (а;.л,), )),

е)=-/0-(и„...,и„ IV,, / = 1,...Ду=1,...,Л^,

^ .')="/, (»„-."з».».,^5, К) + Зп<р{^ ,д„...,д,\

/ = 1,...,5,7 = 1,..„Л/,

<=1

.м

дополненную начальными условиями. Здесь Ц - коэффициенты диффузии ¡-х ионов в межклеточном пространстве, «0, - значения

концентраций ионов в физиологических условиях, Xjir) ~ характеристическая функция области coj вокруг j-й клетки, в которой функция источника предполагается отличной от нуля, v - объем областей со j (считается одинаковым для всех j), In,, la, - нелинейные функции, описывающие ионные токи через мембрану нейронов и астроцитов в соответствии с моделью М.-А. Dronne, J. P. Boissel, E. Grenier1, ajk -коэффициент синаптической связи между j-м и к-м нейроном, z, -постоянный заряд /'-го иона, Q„, Qa - постоянные величины, А, и ç> -нелинейные функции, описывающие изменение переменных <7, и w, в соответствии с моделью T. Hofer2.

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

Заметим, что уравнения в частных производных относительно м,(г,/) в системе (1) с учетом начальных и граничных условий могут быть заменены на интегральные уравнения

"Хг>')-"о/ = J^s Jûfrgj(M|,...)t<5,w1,...,w5,F)G(r>s,i-r) ,

Я) Lo

где C(r,s,i-r) - функция Грина. Воспользовавшись видом функции g из (1) и предполагая, что внутри каждой области ш, /(r,f)« f(r,,t\ Vr e получим

Î2)

Для решения (2) и для целей моделирования нам достаточно определить значения и,(гу,<), _/' = 1,...,Л' + Л/. Кроме того, решая уравнение относительно значений ыДг^,?) около у'-й клетки, можно считать все остальные клетки точечными источниками. Тогда

j !

v », Lo

X»M ' *=I. k*j о

i = l,...,5,y = \,...,N + M.

(3)

Dronne M.-A., Boissel J. P., Grenier E. A mathematical model of ion movements in grey matter during a stroke. Journal of Theoretical Biology, 240, 599-615 (2006)

2 Hofer T., Venance L, Glaume C. Control and plasticity of interceiiufar calcium waves in astrocytes: a modeling approach. Journal of Neuroscience, 22,4850-4859 (2002)

Проведена оценка арифметической сложности численного решения уравнений (3). Показано, что количество необходимых арифметических операций составляет o[(N + М)2Т2), где Т — число шагов по времени, на которые надо разбить промежуток времени от 0 до t для нахождения численного решения. Предложен метод ускорения решения неоднородного уравнения диффузии в бесконечном пространстве, позволяющий значительно сократить это количество.

Рассмотрим следующую вспомогательную задачу в области Л3:

f^(r,r)=£Mr,i)+/(r,0, г е. Rd,t >0, at

•«M»-=0, ' > 0, (4)

и(г,0) = 0, гей,.

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

I

u{T,t)=\dzjds /(s,r)G(r,s,/-r)

о п

Пусть функция / отлична от нуля в конечной области П, и необходимо найти решение задачи (4) в области Л, размеры которой много меньше размеров области ii. Выделим некоторую область П0 (Лей,), размер которой много меньше размера области П. Внутри области П0 оставим функцию источника/в прежнем виде. Область ПШ„ разобьем на конечное число подобластей n„i = \,...,М, в каждой из которых заменим функцию источника на один точечный источник, расположенный в г, (центре области Q,), мощностью

/,(/) = Jrfs /(s,/), i = l„.„Л/.

п,

Тогда решение задачи с измененной функцией источника будет иметь вид

ff(r,i)=JrfrjAG(r,s,/-r)/(i,r)+|;}rfrG(r,s,l-r)7(r) (5)

О £!0 <=' О

В работе доказана Теорема. Погрешность

£= sup sup)/(r,/)-S(r,f)j

lc[/,,r,] геЛ

предложенного метода ускорения вычислений решения уравнения диффузии в области Л имеет оценку сверху

г

е< F sup supraaxfdrmax|G(r,r*,/-r)-G(r,r,,/-rl|

"Л ' i «п.1

где /г= 5иР П/(м)<&, а г" - либо наиболее близкая, либо наиболее

удалённая от области Л точка области О,

Были проведены численные расчёты для модельной двумерной задачи. Предположили, что области Л, П0, С2 представляют собой концентрические окружности с центром в точке г. Области О,, г = 1 ,-.,М, определили как сектора колец с центром в г с разбиением на Nр областей по координате р и Л^ областей по координате <р (М = ИрИЯ1). Выбрали произвольные значения параметров £>, р0, Я, NЫ<р, Г. Получили зависимости погрешности метода е от параметров задачи. Показано, что е резко убывает, стремясь к нулю, при возрастании р0. Погрешность также убывает, стремясь к ненулевому значению при возрастании Nр или N9, и возрастает как функция от Т, но возрастание при больших Т резко замедляется.

Количество операций, требуемых для расчёта решения "огрубленной" задачи, может быть оценено как о(р„). В то же время, погрешность метода с ростом р02 убывает. На Рис. 2 приведена диаграмма, иллюстрирующая процесс выбора оптимального р„.

0.14 :

0.12; (ПО; 0.08; 0.(16- 1

0 .'2 20 44 60 ХО НЮ 0

1 О

Рис. 2. Зависимость погрешности метода (е) от времени численного счета (пропорционально

Численное решение по формуле (5) потребует хранения в памяти компьютера значений всех переменных в каждый момент времени от 0 до I. Чтобы избежать этой необходимости, во второй главе также предложено заменить в формуле (3) нижний предел интегрирования с 0 на /-/*, если 1>г . Предложена процедура поиска 1 при заданном уровне погрешности решения е.

Система уравнений (1) с учётом вышесказанного может быть переписана следу ющим образом:

М I

Л

1=1

; = 1,...,5, у = 1,...,А/, ЯДи„...,и5, тс5)91,...,9з,К); » = 1....Д у = 1.....А/,

Л

(-1

+ начальные условия.

где индекс / означает суммирование по А/, источникам, появившимся в результате применения описанного метода ускорения решения уравнения диффузии в точке гу.

Система (6) включает в себя интегральные уравнения относительно иДг./»0> обыкновенные дифференциальные уравнения относительно w|{nj>t), и'Да^г), и Я+М алгебраических уравнений,

позволяющих определить динамику, переменных , ?)

входящих в систему уравнений неявным образом.

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

воспользуемся методом Рунге-Кутга, неявным по переменным У(пу,(), и явным по остальным переменным. Таким образом, выразим значения <"(«,), ^,"'(«,), как функции Теперь из

алгебраических уравнений имеем

Решив систему методом Ньютона, получим значения V'*r(ny), Далее, по уже заданным функциональным зависимостям определим ^'(иД Аналогично методом Рунге-Кутта можно

определить Наконец, для определения m'+t(iv) воспользуемся

верхними уравнениями из (6), заменив интегрирование по времени на суммирование по дискретным шагам, соответствующим промежутку времени от t — t" до t, после чего все переменные на новом шаге по времени определены.

В третьей главе представлен программный комплекс, написанный на языке Pascal в среде Lazarus IDE v0.9.30, реализующий численное решение модели, приведена блок-схема модели, обсуждаются вопросы реализации программного комплекса на суперкомпьютерах. Приведены и обсуждены результаты компьютерного моделирования.

Подробно рассматривали модельную ситуацию, когда число нейронов N = 60, число астроцитов М = 60, и клетки расположены вдоль одной координаты. Неизвестные параметры системы определили, найдя стационарное состояние модели в физиологических условиях. Далее считали, что первый нейрон и первый астроцит имеют недостаток кровоснабжения (вид функций In,, Iaj изменен), а остальные клетки находятся в зоне достаточного кровоснабжения. Исследовали динамику изменения мембранных потенциалов нейронов V. Распространение зоны поражения во время инсульта представлено на Рис. 3.

ЫЛ1.1АП L'П РТ1/Ц

0 20 40 60 80 100

время, мин

Рис. 3. Динамика изменения мембранных потенциалов V нейронов, находящихся на разном расстоянии от эпицентра поражения. Результаты моделирования.

Исследовали также динамику изменения основных маркеров развития зоны поражения во время инсульта (концентрации К+, glu" в межклеточном пространстве, Na" в нейронах, Са2+ в астроцитах). Результаты (не показаны) подтвердили ряд существующих теоретических предположений, в соответствии с которыми необратимые процессы в эпицентре зоне поражения начинаются с увеличения внутриклеточной

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

В четвертой главе рассматривается задача о математическом моделировании электронного транспорта в мембране тилакоида во время световой фазы фотосинтеза.

Тилакоид - компартмент внутри хлоропласта, ограниченный мембраной. Форма тилакоида напоминает стопки плоских дисков, соединённые друг с другом через диски большего диаметра или трубки. На мембране тилакоида расположен ряд белковых комплексов (см. Рис. 4, а): фотосистема II (ФС2), цитохромный b6f комлекс (цит b6f), фотосистема 1 (ФС1).

Во время световой фазы фотосинтеза под действием света две компоненты ФС2 (акцепторы QA и QB) переходит в восстановленное состояние, характеризующееся наличием свободного электрона. Подвижный переносчик пластохинон (ПХ), диффундируя во внутримембранном пространстве, принимает электрон от ФС2 и доставляет его к цит b6f. Далее, подвижный переносчик пластоцианин (Пц), перемещаясь во внутритилакоидном пространстве, переносит электрон от 4HTb6f к ФС1. Наконец, переносчик ферредоксин (Фд) принимает электрон от ФС1 и перемещается во внешнем пространстве, отдавая электрон на белок ФХР, связанный с цит b6f или на белок ФНР, локализованный во внешнем пространстве (Рис. 4, б).

а б

Рис. 4. а) Схематическое представление пространственного строения тилакоида; б) Схематическое изображение цепи электронного транспорта в мембране тилакоида.

Важными свойствами системы являются подвижность белковых переносчиков, асимметрия расположения комплексов в мембране (ФС2 сосредоточены в «стопке дисков», ФС1 - в межгранальном диске), а также локализация процессов электронного транспорта. Минимальной пространственной структурой для построения модели, учитывающей эти свойства, является одномерный отрезок. Таким образом, структура тилакоида проецируется на чёрную ломаную кривую, изображённую на Рис. 4, а.

Математическая модель формулируется как начально-краевая задача диффузии на отрезке с граничными условиями второго рода:

*б(о,4 /=t>о,

at дх

^|„о,= 0, i = \,...,N, t> О, дп 1

, M,(x,0) = u0(4 X<E[Q,L) (7)

Л, п,

i = l,...,N,t> О,

где L — длина отрезка, N — число компонент модели, ut(x,t), / = \,...,N -функции, описывающие распределение концентрации каждой компоненты модели вдоль отрезка в момент времени /, к, — коэффициент диффузии /-й компоненты (дляусловно неподвижных компонент к, = 0, vJt j - \,...,М -скорости химических реакций на отрезке.

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

Численное решение системы разностных уравнений было реализовано с помощью программного комплекса, написанного на языке Pascal в среде Lazarus IDE v0.9.30. Адекватность модели проверяли, сравнивая результаты вычислительного эксперимента с данными натурного эксперимента, проведённого Т. К. Анталом на кафедре биофизики биологического факультета МГУ. В экспериментах измерялась интенсивность флуоресценции, наблюдаемой от листьев растений, освещаемых после длительного пребывания в темноте, которая может быть вычислена из данных нашего моделирования по формуле

а,

F{t) =-i-,

af+ad + a4Unx (t) + apUgA (t)

LI 2 LI 2

Unx{f) = J "i t)dx, UgA(t)= \u2{x,t)dx. о 0

где ad, aq, ap - константы, ut{x,t), w2(x,/) - переменные,

описывающие концентрацию восстановленного пластохинона и восстановленного QA.

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

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

I, мс

Рис. 5. Зависимость интенсивности флуоресценции Р от времени I по результатам модели и натурного эксперимента, проведённого Т. К. Анталом на кафедре биофизики Биологического факультета МГУ.

В заключении перечислены основные результаты и выводы, полученные в диссертации:

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

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

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

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

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

Публикации в изданиях из Перечня ВАК:

1. Макаров С. С., Исаева A.B., Грачёв Е. А., Сердобольская М. JI. Ускорение вычислений при решении неоднородного уравнения диффузии с помощью перенормировочного метода // Вычислительные методы и программирование, 2012, т. 1, с. 239-245.

2. Макаров С. С., Грачёв Е. А., Антал Т. К. Математическое моделирование электрон-транспортной цепи в тилакоидной мембране с учетом пространственной гетерогенности мембраны // Математическая биология и биоинформатика, 2012, т. 7, выпуск 2, с. 508-528.

3. Макаров С. С., Джебраилова Ю. Н., Грачёва М. Е., Грачёв Е. А., Кочетов А. Г., Губский JI. В. Моделирование группы нейронов и астроцитов в условиях ишемического инсульта // Журнал неврологии и психиатрии им. С. С. Корсакова, 2012, N 8, выпуск 2, с. 59-62.

Публикации в сборниках тезисов научных конференций:

4. Макаров С. С., Грачев Е. А. Моделирование межклеточной кальциевой волны в плоской ткани при непостоянной проводимости щелевых контактов // Тезисы докладов семнадцатой международной научной конференции "Математика. Компьютер. 0бразование"-2010. Секция S3. Анализ сложных биологических систем. Эксперимент и модели. Дубна, 2530 января 2010 г.

5. Макаров С. С. Моделирование группы нейронов и астроцитов в условиях ишемического инсульта // Тезисы докладов XVIII Международной научной конференции Ломоносов-2011. Секция "Физика". Подсекция "Биофизика". Москва, 11-15 апреля 2011 г.

6. Макаров С.С., Антал Т. К., Плюснина Т. Ю., Грачёв Е. А. Математическая модель индукции кривой флуоресценции хлорофилла а с учетом электронного транспорта в тилакоидной мембране // Тезисы докладов девятнадцатой международной научной конференции "Математика. Компьютер. 0бразование"-2012. Секция S3. Анализ сложных биологических систем. Эксперимент и модели. Подсекция "Молекулярная и клеточная биофизика". Дубна, 30 января - 4 февраля 2012 г.

Подписано в печать: 25.04.2013 Тираж: 70 экз. Заказ №981 Отпечатано в типографии «Реглет» . Москва, Ленинградский проспект д.74 (495)790-47-77 www.reglet.ru

Московский государственный университет имени М. В. Ломоносова

Физический факультет Кафедра компьютерных методов физики

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

Макаров Сергей Сергеевич

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ДИНАМИЧЕСКИХ ПРОЦЕССОВ В ПРОСТРАНСТВЕННО-НЕОДНОРОДНЫХ БИОЛОГИЧЕСКИХ СИСТЕМАХ

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

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель: к.т.н., доцент Грачёв Е. А.

Москва-2013

ОГЛАВЛЕНИЕ

стр.

ВВЕДЕНИЕ 4

ГЛАВА 1. РАЗРАБОТКА МАТЕМАТИЧЕСКОЙ МОДЕЛИ РАЗВИТИЯ

ИШЕМИЧЕСКОГО ИНСУЛЬТА 10

1.1. Моделируемая система 10

1.1.1. Система клеток головного мозга 10

1.1.2. Развитие ишемического инсульта 14

1.2. Обзор математических моделей 16

1.3. Принципы настоящей модели 20

1.3.1. Моделирование ионных токов через клеточную мембрану 22

1.3.2. Моделирование кальциевой сигнализации между астроцитами 23

1.3.3. Моделирование синаптической связи между нейронами 25

1.3.4. Моделирование диффузии ионов в межклеточном пространстве 27

1.4. Общая математическая постановка задачи 27

1.5. Заключение (по Главе 1) 30

ГЛАВА 2. ИССЛЕДОВАНИЕ И МЕТОДЫ ЧИСЛЕННОГО РЕШЕНИЯ

МОДЕЛИ РАЗВИТИЯ ИШЕМИЧЕСКОГО ИНСУЛЬТА 32

2.1. Исследование задачи о диффузии вещества в межклеточном пространстве 32

2.2. Численный метод решения системы уравнений математической модели 35

2.2.1. Формулировка метода 35

2.2.2. Оценка арифметической сложности метода 37

2.3. Метод ускорения решения неоднородного уравнения диффузии 38

2.3.1. Формулировка метода 39

2.3.2. Оценка погрешности метода для задачи произвольной размерности 42

2.3.3. Численная оценка погрешности метода для модельных задач 44

2.3.4. Применение метода для моделирования развития ишемического инсульта 47

2.4. Ускорение численного решения с помощью пренебрежения малыми значениями функции Грина 49

2.5. Заключение (по Главе 2) 52

ГЛАВА 3. КОМПЬЮТЕРНАЯ РЕАЛИЗАЦИЯ МОДЕЛИ РАЗВИТИЯ

ИШЕМИЧЕСКОГО ИНСУЛЬТА. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ 55

3.1. Принципы создания комплекса программ 55

3.1.1. Реализация модели на персональном компьютере 55

3.1.2. Принципы реализации модели на параллельных вычислительных машинах 60

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

3.2.1. Физиологические условия

3.2.2. Патологические условия 3.3. Заключение (по Главе 3)

61 63 72

ГЛАВА 4. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЭЛЕКТРОННОГО

ТРАНСПОРТА НА БИОЛОГИЧЕСКОЙ МЕМБРАНЕ 74

4.1. Описание моделируемой системы и обзор моделей 74

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

4.2.1. Основные положения 78

4.2.2. Компоненты модели 79

4.2.3. Сравнение с экспериментальными данными 87

4.2.4. Реализация модели 88

4.3. Результаты и обсуждение 90

4.4. Заключение (по Главе 4) 97

ЗАКЛЮЧЕНИЕ 98

ПРИЛОЖЕНИЕ 1. Описание ионных каналов, рассмотренных в модели

развития ишемического инсульта 101

ПРИЛОЖЕНИЕ 2. Уравнения и параметры модели кальциевой

сигнализации в астроцитах 107

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

ВВЕДЕНИЕ

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

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

Первый подход заключается в усреднении переменных по пространству и рассмотрению системы как гомогенной. Такие модели называют кинетическими. Математически они представляют собой совокупность систем обыкновенных дифференциальных уравнений относительно усредненных значений исследуемых переменных. В частности, в недавних работах М.-А. Dronne, G. Chapuisat, Е. Grenier была представлена кинетическая модель поведения клеток головного мозга в условиях инсульта. Модель хорошо описывала изменение состояния клеток в эпицентре зоны поражения, в условиях сниженного кровотока. Однако медицинский интерес не ограничивается этими клетками. Важным является и описание пространственного развития поражения на соседние области, в том числе и на те, кровоток в которых не нарушен. Зачастую необходимо определить пространственные характеристики зоны поражения, и в этом случае кинетические модели являются недостаточно информативными.

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

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

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

Третий подход основан на дискретизации системы и переходе от переменных, описывающих концентрации, к отдельным дискретным элементам. Методы, используемые при таком моделировании, включают в себя методы Монте-Карло, методы молекулярной динамики, методы частиц и т.д. На настоящий момент не существует дискретных моделей развития ишемического инсульта. Дискретные стохастические модели, описывающие динамику отдельных переносчиков электронов в липидной мембране, появились относительно недавно и активно развивались Г. Ю. Ризниченко, А. Б. Рубиным, Е. А. Грачёвым и другими. Достоинством дискретных моделей является удобство их реализации в виде параллельных программ, однако большое общее количество элементов, тем не менее, может сделать трудоемкими и такие модели.

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

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

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

—

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

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

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

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

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

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

- модели электронного транспорта между белковыми комплексами в липидной мембране тилакоида во время фотосинтеза.

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

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

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

Соответствие диссертации паспорту научной специальности. Работа соответствует паспорту специальности 05.13.18 - «Математическое моделирование, численные методы и комплексы программ», так как включает в себя разработку, исследование и редукцию новых

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

Основные результаты работы, выносимые на защиту

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

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

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

4) Получены результаты расчетов созданных программ, по существующим экспериментальным данным доказана адекватность математических моделей.

Научно-практическая значимость работы. Предложенные в диссертации математические модели биологических систем могут быть использованы для решения прикладных задач биологии и медицины (диагностика ишемического инсульта и др.).

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

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

Личное участие автора в выполнении работы. Постановка биофизической задачи о моделировании поведения группы нейронов и астроцитов в условиях ишемического инсульта принадлежит автору совместно с научным руководителем, а также докторами медицинских наук Кочетовым А. Г. и Губским Л. В. Формулировка математической модели проведена автором лично, с консультациями Грачёва Е. А.

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

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

Апробация результатов. Материалы диссертации докладывались на международной конференции «Математика. Компьютер. Образование» в Объединенном Институте Ядерных Исследований, г. Дубна, январь 2010 года и январь 2012 года; на международной конференции студентов, аспирантов и молодых ученых по фундаментальным наукам "Ломоносов-2011" на физическом факультете МГУ, апрель 2011 года; на рабочем семинаре сектора информатики и биофизики сложных систем кафедры биофизики биологического факультета МГУ под руководством д.ф.-м.н. Г. Ю. Ризниченко, декабрь 2011 года и январь 2012 года; на научном семинаре кафедры высшей нервной деятельности биологического факультета МГУ, март 2012 года; на научном семинаре кафедры биофизики физического факультета МГУ, апрель 2012 года; на научном семинаре кафедры вычислительной математики Московского Физико-Технического Института (МФТИ), июнь 2012 года; на научном семинаре "Обратные задачи математической физики" под руководством д.ф.-м.н. А. Г. Яголы, д.ф.-м.н. А. Б. Бакушинского и д.ф.-м.н. А. В. Тихонравова, проводящемся в НИВЦ МГУ, октябрь 2012 года; на научном семинаре кафедры компьютерных методов физики физического факультета МГУ, ноябрь 2012 года; на научно-методологическом семинаре НИВЦ МГУ под руководством д.ф.-м.н. А. В. Тихонравова, декабрь 2012 г. на научном семинаре института математических проблем биологии РАН под руководством д.ф.-м.н. М. Н. Устинина, декабрь 2012 года; на научном семинаре в государственной клинической больнице №31 под руководством д.м.н. Л. В. Губского, март 2013 года.

Структура и объем диссертации. Диссертация состоит из титульного листа, оглавления, введения, четырех глав, заключения, списка литературы, содержащего 93 наименования, и двух приложений. Общий объем диссертации составляет 117 страниц, в том числе 29 рисунков и 10 таблиц.

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

Публикации в журналах из перечня ВАК:

1. Макаров С. С., Исаева А. В., Грачёв Е. А., Сердобольская М. Л. Ускорение вычислений при решении неоднородного уравнения диффузии с помощью перенормировочного метода // Вычислительные методы и программирование, 2012, т. 1, с. 239-245.

2. Макаров С. С., Грачёв Е. А., Антал Т. К. Математическое моделирование электрон-транспортной цепи в тилакоидной мембране с учетом пространственной гетерогенности мембраны // Математическая биология и биоинформатика, 2012, т. 7, выпуск 2, с. 508-528.

3. Макаров С. С., Джебраилова Ю. Н., Грачёва М. Е., Грачёв Е. А., Кочетов А. Г., Губский Л. В. Моделирование группы нейронов и астроцитов в условиях ишемического �