CC BY
37УДК 517.958:57
Брацун Дмитрий Анатольевич
доктор физико-математических наук, заведующий кафедрой теоретической физики и компьютерного моделирования
Захаров Андрей Павлович
кандидат физико-математических наук, заведующий лабораторией кафедры теоретической физики и компьютерного моделирования
ФГБОУ ВПО «Пермский государственный гуманитарно-педагогический университет», Пермь, Россия 614990, Пермь, Сибирская, 24, (342) 238-63-64,
e-mail: bratsun @pspu. ru
ПРОСТРАНСТВЕННО-РАСПРЕДЕЛЕННЫЕ СТОХАСТИЧЕСКИЕ КОЛЕБАНИЯ С ЗАПАЗДЫВАНИЕМ В МНОГОКЛЕТОЧНОЙ
СИСТЕМЕ*
Dmitry A. Bratsun
DS, Head of the Theoretical Physics Department Andrej P. Zakharov
PhD, Head of the Laboratory at the Theoretical Physics Department
Federal State Budget Educational Institution of Higher Professional Education «Perm State Humanitarian Pedagogical University» 24, Sibirskaja, 614990, Perm, Russia, e-mail: bratsun@pspu. ru
SPATIALLY DISTRIBUTED DELAY-INDUCED STOCHASTIC OSCILLATIONS IN A MULTI-CELLULAR SYSTEM*
Аннотация: исследуется влияние внутреннего шума и запаздывания на структурообразование поля концентрации протеина в модели квазидвумерной ткани эпителия. Колебания вызываются процессами транскрипции/трансляции в эпителиальных клетках, обменивающихся химическими и механическими сигналами. Колебания описываются на основе одногенной модели авторепрессора с постоянным значением запаздывания по времени. В работе дается как детерминистское, так и стохастическое описание системы. Обнаружено, что совместное действие запаздывания, шума и пространственного обмена сигналами между клетками может приводить к образованию паттернов даже в том случае, когда детерминистское описание предсказывает устойчивое состояние.
© Брацун Д.А., Захаров А.П., 2015
* Работа выполнена при финансовой поддержке Министерства образования Пермского края (грант ^26/00.4) и гранта РФФИ (14-01-96022р_урал_а).
Ключевые слова: генная регуляция, стохастические колебания, запаздывание по времени, синхронизация.
Abstract: We explore the joint effect of the intrinsic noise and time delay on the spatial pattern formation within a multi-scale mobile lattice model of the epithelium. The protein fluctuations are driven by transcription/translation processes in epithelial cells exchanging chemical and mechanical signals and are described by a single-gene auto-repressor model with constant delay. Both deterministic and stochastic descriptions are given. We found that time delay, noise and spatial signaling can result in the protein pattern formation even when deterministic description exhibits no patterns.
Key words: gene regulation, stochastic oscillations, time delay, synchronization.
В последнее время всё больший интерес вызывают динамические системы с запаздывающим аргументом. Диапазон приложений включает в себя популяционную динамику и социальные процессы [13], нелинейные химические реакции и процессы генной регуляции [8,16], поведение систем с автоматическим управлением [12,14], механику жидкости [15] и т.д. Запаздывание может быть обусловлено самыми различными причинами, например, ограниченностью скорости распространения взаимодействия (электрического или светового сигнала), наличием инерционности некоторых элементов (при управлении с обратной связью) или существованием цепочек многоэтапных последовательных реакций с известным результатом в конце (при транскрипции генов).
Известно, что для изучения свойств динамической системы можно использовать как детерминистское, так и стохастическое описание. Так как большинство реальных динамических процессов являются стохастическими, понятен интерес исследователей к изучению стохастических дифференциальных уравнений с запаздыванием [19,2,3]. Были разработаны специальные методы решения таких уравнений. Например, в работе [2] система запаздывающих дифференциальных уравнений была сведена к двумерному отображению, а затем возмущена белым шумом. Техника применения мастер-уравнения для распределения плотности вероятности была использована в работе [3] для дискретной задачи о «движении пьяницы». Применение приближения «слабого» запаздывания к исходным системам дифференциальных уравнений было продемонстрировано в [19]. Отметим серию работ [20,4,7]. В работе [7] рассматривалась простая модель с двумя равнозначными состояниями равновесия, между которыми система может совершать запаздывающие переходы с определенной вероятностью. Было обнаружено, что корреляционная функция сигнала имеет ярко выраженную резонансную структуру. Авторы подчеркивают, что природа этого явления отлична от хорошо известного стохастического резонанса, который возникает, например, в той же бистабильной зашумленной системе с периодическим внешним воздействием. В данном случае система является автономной
и ведущий резонансный пик у корреляционной функции определяется временем запаздывания. В дальнейшем авторы рассматривали различные аспекты этого и примыкающего к нему явлений. В работе [20] рассматривалась система, обобщающая результаты [7] на случай связанных между собой сети бистабильных элементов с глобальным процессом запаздывания. Вопросы управления системой при помощи обнаруженного эффекта изучались в [4]. В целом можно заключить, что взаимодействие шума и запаздывания изучалось на примере искусственных или упрощенных механических систем. При этом шум, как правило, учитывался в рамках аддитивного слагаемого с наперед заданными свойствами в правой части дифференциальных уравнений (метод Ланжевена). Вопрос о том, как проявляет себя эффект в реальных, системах оставался открытым.
В работах автора [8-12,16, 17] вопросы взаимодействия запаздывания и стохастических флуктуаций впервые были рассмотрены для процессов транскрипции/трансляции генов. Обнаружен эффект возбуждения в подкритической области колебаний за счет взаимодействия шума и запаздывания обратной связи [1,9]. Экспериментально эффект был обнаружен в работах [11,12]. Приложения эффекта для конкретных биологических систем были рассмотрены в [9,10]. Как оказалось, обнаруженный эффект имеет особенно сильные проявления в генетических системах, где небольшое число вовлеченных в процесс молекул может приводить к сильным флуктуациям поля концентрации протеина. В работах [16,17] впервые была предложена модификация алгоритма Гиллеспи [18], используемого для прямого численного моделирования малоразмерных стохастических дискретных процессов, на случай систем с запаздыванием. За время появления работ [16,17] предложенный алгоритм приобрел популярность среди исследователей: вышло большое количество работ (более 300), развивающих приложения алгоритма к разным системам.
Важнейшим дальнейшим развитием алгоритма является построение его пространственно-распределенной версии. Несколько лет назад в работе [21] была высказан мысль о том, что пространство остается последней крепостью в организации стохастических вычислений в биологических системах. В последнее время в математической биологии появилось большое количество экспериментальных данных о пространственно-временной динамике концентраций протеинов, хотя абсолютное большинство теоретических работ по инерции заняты изучением только временной составляющей этой динамики. В последние несколько лет ситуация стала меняться [1,6,22], однако редкие примеры рассмотрения стохастических процессов в пространстве посвящены Марковским системам. Причина, в общем-то, понятна - теоретически обоснованного алгоритма расчета стохастических немарковских систем, обоснованного в стиле работы [16], до сих пор не предложено.
В данной работе мы рассматриваем эффект взаимодействия внутреннего шума и запаздывания по времени в рамках многоклеточной хемомеханической
модели эпителия, предложенной в работах [5,10]. Использование латтис-модели пространственной среды позволяет, в какой-то мере, обойти вопрос о последовательной реализации стохастического расчета пространственно -диффузионного канала и использовать хорошо отработанный алгоритм, предложенный в работах [16,17]. Конечно, вопрос о строгом рассмотрении алгоритма стохастического расчета просачивания молекул через мембраны клеток остается открытым - здесь же мы используем простейшую диффузионную модель для межклеточных сигналов. Главная же цель работы -рассмотрение эффекта стохастических подкритических колебаний, обнаруженных в [16], в случае пространственной эволюции системы.
Опишем кратко основные особенности модели эпителия, предложенной в [5] и развитой для моделирования раковых образований в [10]. Эта модель включает расчет динамики отдельных клеток, представленных в виде многоугольников с разным числом вершин (система откалибрована так, чтобы наиболее вероятной формой клетки являлась гексагональная ячейка, хотя появление других видов многоугольников также возможно). Клетки плотно примыкают друг к другу, образуя сплошную двумерную поверхность эпителия. Модель обладает целым набором свойств, позволяющих удачно имитировать поведение реальной эпителиальной ткани:
- возможность изменения размеров клеток в процессе эволюции ткани и изменения локальных механических свойств среды;
- возможность роста общего количества клеток в системе посредством их деления в определённых условиях эволюции;
- возможность перемещения клеток в общей массе эпителия посредством механизма интеркаляции;
- расчет для каждой клетки динамики концентрации веществ, участвующих в регуляции жизнедеятельности ткани;
- обмен химическими сигналами, осуществляемый между соседними клетками эпителия через общую границу.
Таким образом, каждая клетка в модели испытывает ряд хемомеханических воздействий, под влиянием которых она эволюционирует вместе со всей системой. Так как внутри клетки поля зависят только от времени, то мелкая структура пространственно-распределенных эффектов, связанная с неоднородностью полей внутри клеток, в модели не определяется. Однако на расстояниях больших, по сравнению с размерами одной клетки, пространственное структурообразование проявляет себя в полной мере. Разработанную модель можно классифицировать как дискретную сложную систему, демонстрирующую коллективные эффекты, с индивидуальным поведением отдельных элементов.
Опишем теперь механику самих клеток. Потенциальная энергия всей эпителиальной ткани может быть записана следующим образом:
E = 11 (кL + n(A - A0)2) 2
(1)
' cells
где ь и А - соответственно периметр и площадь клетки. Первое слагаемой в выражении (1) описывает действие сил, стремящихся сократить периметр каждой клетки, а второе выражает сопротивление клетки действию сил растяжения и сжатия, и стремлению её сохранить площадь А ■ Соответствующие коэффициенты кип являются важными параметрами задачи. Эпителиальная ткань может эволюционировать при перемещении узлов клеток. Механическая сила, действующая на ] -й узел, определяется стандартным образом как производная от потенциальной энергии:
F = --SE- + F st F SR. j
(2)
где и. - радиус вектор узла, ^ - белый шум, имеющий нулевое среднее значение.
Таким образом, уравнение движения клеток эпителия можно записать:
^ = h (frh\-f0) f
mech
(3)
где Н - функция Хэвисайда, ^ - параметр, определяющий пороговую силу, ниже которой узел остается неподвижным. Как видно из уравнения (3), для ансамбля клеток, составляющего ткань эпителия, лучше подходит механика Аристотеля. Это происходит по той же причине, почему для пористой среды лучше подходит закон Дарси - сопротивление среды по отношению к движению элементов настолько велико, что силы определяют скорости, а не ускорения [20].
Рис. 1. Схематическое описание процесса генной транскрипции. Красным цветом отмечены мономеры протеина, ответственного за внутриклеточные колебания (в черно-белом варианте эти молекулы более темные). Синим цветом отмечены мономеры транспортного протеина. Транспортный протеин диффундирует из клетки в клетку, осуществляя обмен сигналами.
Рассмотрим теперь модель процесса генной транскрипции, проходящего внутри каждой клетки (рис.1). Для описания внутриклеточных колебаний мы используем простейшую модель транскрипции одного гена с запаздывающей отрицательной обратной связью [16]. Это весьма популярный мотив при изучении регуляторных генных цепочек. Например, аналогичная модель без запаздывания была рассмотрена подробно как в рамках детерминистского, так и стохастического описания в работе [21]. Пусть протеин в клетках может существовать только в мономерной и димерной формах и переменные X, У обозначают концентрации мономеров и димеров соответственно. Тогда переходы между мономерами и димерами можно записать как
х + х > у , (4)
у —> х + х, (5)
где к±а, - скорости реакций перехода из одной формы в другую. Предположим далее, что протеин может синтезироваться и распадаться со скоростями А и В соответственно:
Д —^ д+т + х'+т, (6)
х' —^0'. (7)
Работу промотора гена можно учесть введением оператор-сайта, который может останавливать или возобновлять транскрипцию, - это достигается введением специальной функции д е{ Д, Д}, принимающей значение Д в случае открытия оператор-сайта ид в случае его закрытия. С формальной точки зрения состояние оператор-сайта можно рассматривать как дополнительный реагент в системе. Синтез белка начинается при условии Д(') = 1. Если же он закрыт, то Д(') = 0 и производство белка невозможно. Отметим, что реакция (6) является запаздывающей, и белок синтезируется только через время т после начала транскрипции. Кроме того, состояние оператор-сайта зависит от концентрации димеров У:
Д0 + У д, (8)
Д —До + У, (9)
где к - скорости переключений состояний оператор-сайта. Таким образом, в системе реализуется отрицательная обратная связь: при увеличении концентрации синтезированного протеина вероятность перехода оператор-сайта в состояние «закрыто» (8) возрастает, что влечет за собой прекращение транскрипция гена, а значит, и синтеза протеина (6).
Рассмотрим теперь часть модели, отвечающей за межклеточные обмены (см. рис.1). Для описания взаимодействия клеток мы вводим транспортный белок, концентрация Т которого регулируется белком х посредством положительной обратной связи:
0 >г, (10)
где А - скорость реакции. При такой записи 7-протеин синтезируется только в случае «закрытия» оператор-сайта Т-гена, причем синтез происходит немедленно, без всякого запаздывания:
В% + X —Б? , (11)
БТ —Б0 + X. (12)
Здесь к±2 - скорости переключений состояний оператор-сайта Т-гена. Для простоты регулирующих фактором транскрипции был выбран Х-протеин в мономерной форме.
В конце концов мы предполагаем, что при попадании транспортного протеина в определенную клетку он превращается в Х-мономер со скоростью В :
Т—^Х, (13)
причем, происходит это с числом копирования N (это означает, что синтезируется сразу N копий мономера). Все вместе реакции (4-13) определяют кинетику процессов генной регуляции в клетках.
Рассмотрим свойства модели в рамках детерминистского описания. Основное допущение, которое здесь обычно делается, состоит в предположении медленности процессов синтеза и деградации протеина (6, 7) по сравнению с прямой и обратной реакцией димеризации (4, 5), а также реакциями закрытия/открытия оператор-сайта (8, 9). Это означает, что динамика системы быстро приходит к состоянию локального равновесия [8,9,17]:
У = вX2, Б. = гЪБо1Х12, БТ = сБТХ , (14)
где в = к+а /к_а, 5 = к+1 / , а = к+2 / & 2, а индекс соответствует номеру клетки. С учетом выражений (14) после несложных преобразований получим модельные уравнения:
(1 + 4 в X;) ^ = А + ВТ (г) - ВХ; (0, (15)
ш 1 + в5Xi (г - т)
^ = ШгсХ^ - ВТ (г) + £ а ь (Т (г) - Т (г)), (16)
ёг 1 + аГ;(г) )
где а - коэффициент диффузии, а суммирование в правой части (16) ведется по соседним клеткам. Предполагается, что транспортный протеин диффундирует из клетки в клетку, при этом поток вещества зависит от периметра клетки Ь (см. рис.1). Таким образом, поле транспортного протеина является глобальным для всей ткани.
Рис. 2. Нейтральная кривая для колебательных возмущений на плоскости безразмерных параметров производства и деградации протеина: е = 0,1, 5 = 0,2 (а). Спектр Фурье стохастического сигнала, полученного ниже нейтральной кривой (нижняя точка на рис. 2, а)
при А = 20, В = 5 , т = 6, к+1 = 100, к_1 = 1000, = 200, к_л = 1000 (б)
Для системы, состоящей из одной клетки, уравнение (15) примет вид:
(1 + 4еX)^ = , "ВХ(г) . (17)
йг 1 + е5Х (г - т)
Впервые уравнение (17) рассмотрено в работе [16]. У этого уравнения имеется единственное положение равновесие, которое при определенных значениях параметров может терять устойчивость по отношению к колебательным возмущениям. На рис. 2 (а) приведена нейтральная кривая бифуркации Хопфа в координатах безразмерных скоростей медленных реакций. Колебательная неустойчивость развивается выше нейтральной кривой.
Для рассмотрения пространственных эффектов, возникающих за счет обмена сигналами между клетками, система запаздывающих дифференциальных уравнений (15,16) решалась численно с помощью явного метода Эйлера, чья устойчивость гарантировалась выбором достаточно мелкого шага по времени. Эта процедура была синхронизирована с расчетом механической эволюции ткани в рамках системы уравнений (1-3). В качестве начального состояния системы была выбрана гексагональная форма из 1560 клеток со случайным распределением фазы колебаний и фиксированным значением констант, определяющих свойства механической подмодели: к = 1,0, П = 1,0, = 3л/з/2, ^ = 0,02 .
На рис.3(а) представлены два кадра эволюции поля Х-протеина в рамках детерминистского описания для значений параметров выше нейтральной кривой (верхняя точка на параметрической плоскости рис. 2, а). Нелинейная динамика системы включает спиральные бегущие волны, которые возникают на фоне осциллирующего синхронизированного поля с частотой колебаний приблизительно равной утроенному времени запаздывания.
а
б
Рис. 3. Два последовательных кадра эволюции поля концентрации Х-протеина для моментов времени г = 260 и г = 340 в эпителии, состоящем из 1560 клеток, в рамках детерминистского (а) и стохастического (б) описания выше нейтральной кривой бифуркации Хопфа (верхняя точка, см. рис.2, а). Фиксированные значения параметров при: А = 500, Ат = 500, В = 5 , т = 6,
N = 1, а = 0.05, к+1 = 100, к-1 = 1000, к+2 = 100 , к-2 = 1000, к+ы = 200, к-а = 1000
Обсудим теперь стохастическое описание поведения системы. Одним из самых известных методов численного исследования стохастических малоразмерных реакций является метод Гиллеспи [18], который за время своего появления в 1977 г. стал поистине классическим из -за своей простоты и надежности. Математически метод Гиллеспи воспроизводит решение для распределения вероятности, задаваемого мастер-уравнением. Опишем кратко главную идею метода Гиллеспи и практический алгоритм его использования в численных исследованиях. Предположим, имеется К химических соединений Хух, вступающих в многоэтапную реакцию между собой. Каждый тип реактивного контакта, случающийся между молекулами, будем называть вслед за Гиллеспи «химическим каналом» Я. Пусть таких каналов имеется М штук.
Согласно алгоритму Гиллеспи, многосоставная реакция распадается на элементарные акты. При совершении каждого акта, нужно вычислить время, через которое произойдет следующая элементарная реакция, а также определить канал этой реакции, т.е. выяснить, какие именно компоненты
81
из вектора ^ должны вступить в контакт. Для Марковских процессов, для
которых изначально сформулирован метод Гиллеспи, распределение промежутков времени Лг, проходящих между реакциями, подчиняется распределению Пуассона:
V
Р(Лг) = £«VехР -Лг^«V • (18)
Здесь ^ = еД - вероятность того, что в данный момент реализуется V -й химический канал. Значение формируется благодаря двум факторам. Не углубляясь в тонкости кинетики реакций можно сказать, что е - это скорость реакции. Величина к определяется для каждого химического канала и равняется количеству различных молекулярных соединений из набора Xv. Это число легко определяется кинетикой заданной реакции [18]. Отбор следующей реакции зависит от дискретного распределения
Р^ = V *) = ^ • (19)
V
Классический алгоритм Гиллеспи [11] очень прост. Он включает в себя вычисление по заранее известным кинетическим уравнениям химических реакций (каналов) значений « . После этого генерируются два случайных числа
из отрезка [0,1]. Первое, согласно распределению (18), определяет время наступления ближайшей реакции Лг, а второе, согласно распределению (19), -химический канал V*, который реализуется через время Лг. Система делает прыжок г = г + Лг, и реакция совершается. Обновляется вектор состояния системы Xv, и вся процедура повторяется.
Модификация алгоритма на случай немарковских процессов была предложена в (17,4,7). Главная идея обновления алгоритма заключается в следующем. Предположим, один из химических каналов является
немарковским. В соответствии с распределениями (18) и (19) определяется канал V*, который должен реализоваться, и промежуток времени Лг, через который это должно произойти. Если выпадает реакция V* с запаздыванием, то она не выполняется, а помещается в стек памяти с тем, чтобы быть выполненной в момент времени Лг + т, где т - время запаздывания. Если же выпавший химический канал является Марковским, то прыжок во времени Лг , который должна совершить система, сравнивается со временами наступления отложенных в стек памяти реакций. Если в диапазоне между г и г + Лг никаких реакций не запланировано, то система совершает прыжок на Лг . Если в этот диапазон времени как раз попадает одна из запаздывающих реакций, то последняя процедура отбора игнорируется, а система на самом деле прыгает к моменту времени, которое ранее было уже запланировано,
V
но отложено из-за запаздывания. В данной работе мы расширяем этот алгоритм на случай пространственных реакций в рамках предложенной клеточной модели. Опишем предлагаемый алгоритм численного решения более подробно. В стохастических вычислениях величины Хг. (число мономеров), У{ (число
димеров), т (число молекул транспортного протеина), а также состояния оператор-сайтов считаются целыми. Индекс «г» относится здесь к клеткам эпителия. Тогда в каждой клетке производятся следующие вычисления:
1. Задание начальных значений для Х1, Гг, Тг, Вг, Бь, вТ, вТ. Количество химических соединений К = 7.
Количество химических каналов М = 9. Установка начального времени гг = 0.
Установка начального значения счетчика шагов по времени = 1.
2. Вычисление вероятности реализации для каждого канала а . , р = 1..М. Вычисление значение полной вероятности для каналов а0г = ^ацг.
3. Генерирование двух случайных чисел ги, ^ е[0,1].
4. Вычисление промежутка времени до следующей реакции Аг']осН = - 1п
а0г
5. Проверка на наличие отложенной (запаздывающей) реакции или эволюционного шаг процесса диффузии, которые должны произойти внутри рассматриваемого промежутка времени: г^е1ау е[г, г + Аг^осН] либо
г^ е [г, г + Аг}госк]:
А. Если в указанном промежутке [г, г + Аг^осН] в момент га1е1ау происходит
отложенная реакция, то
а) пункты алгоритма 2-4 игнорируются;
б) система совершает прыжок до момента времени гае1ау;
в) запаздывающая реакция совершается;
г) счетчик реакций обновляется ^ = ^ +1;
д) переход к пункту 2.
Б. Если же в указанном промежутке [г,г + Аг^°с1г] первым происходит акт
диффузии в момент времени г^ , то
а) значение для транспортного протеина обновляется согласно формуле:
аАг^ £ЦТ -Т)
гг
(20)
где квадратные скобки обозначают операцию взятия целого значения. б) механическое состояние обновляется в соответствии с (1-3); д) переход к пункту 6. В. Если же в указанном диапазоне времени события А или Б не происходят, то переходим к пункту 6.
ц*-1 ц*
6. Нахождение канала для следующей реакции ц* согласно ^ а^ < гъаы < ^ а^ .
ц=1 ц=1
7. Если выпавший канал ц* соответствует реакции без запаздывания, то
а) система совершает прыжок до момента времени г + ;
б) вектор решения обновляется в соответствии с выпавшим каналом;
г) счетчик реакций обновляется У. = У. +1;
д) переход к пункту 2.
8. Если выпавший канал ц* отвечает реакции с запаздыванием, то
а) выполнение реакции откладывается до момента времени г + т;
б) переход к пункту 2.
Необходимо отметить, что шаг по времени в уравнении (20) синхронизирован с шагом по времени механической эволюции системы (1-3), который обычно устанавливался равным Аг = 0,05. Так как типичный шаг по времени стохастической подсистемы гораздо меньше (обычно Аг « 0,00001 - 0,001), то возникает проблема состыковки двух численных схем для совместной работы. Один из вариантов решения этой проблемы приведен в описании алгоритма.
Рис. 4. Четыре последовательных кадра эволюции поля концентрации Х-протеина в рамках одного периода колебаний 2 т в эпителии, состоящем из 1560 клеток, в рамках стохастического описания ниже нейтральной кривой бифуркации Хопфа (нижняя точка, см. рис. 2, а). Фиксированные значения параметров те же, что и на рис.3, за исключением А = 20, N = 8
Обсудим теперь результаты расчетов системы в рамках стохастического описания. На рис. 2 (б) приведен спектр Фурье стохастического сигнала ниже бифуркационной кривой (точка отмечена на рис. 2 (а) в его нижней части) для случая одной замкнутой от внешних сигналов клетки. Примечательной особенностью спектра является явное присутствие в сигнале периодичности: первый пик соответствует частоте ю* = п/т « 0,524 . Остальные пики являются гармониками первого пика и удовлетворяют формуле
юп = (2п -1) ю*. (21)
Таким образом, запаздывание и шум могут генерировать периодические колебания даже в той области параметров, где детерминистское описание предсказывает абсолютную устойчивость системы. Это был важный результат работы [16]. Рис. 3 (б) и 4 дают два примера того, каким образом указанный эффект раскрывается в пространственном измерении. Кадры эволюции поля Х-протеина, приведенные на рис. 3 (б), получены выше бифуркационной кривой. В этом случае поведение системы весьма схоже с пространственно-временной динамикой системы в рамках детерминистского описания (см. рис. 3, а). В случае же расчета ниже нейтральной кривой (рис. 4) результат зависит от числа копирования N. При малых значениях этого параметра конечным состоянием эволюции является полностью синхронизированное поле, осциллирующее с частотой ю*. Увеличение числа копирования ведет к возникновению паттерна - на рис.4 (N = 8) хорошо видны два кластера клеток, которые осциллируют в противофазе. Этот же эффект иллюстрирует рис.5, где приведена эволюция распределения клеток по фазам колебаний. Структурообразование в клеточных системах в виде кластеров активно изучается в литературе, так как кластеры рассматриваются как причина для дальнейшей дифференциации клеток.
monomers 3 10 963
Рис. 5. Временная эволюция распределения клеток по фазам стохастических колебаний для значений параметров ниже нейтральной кривой. Хорошо видно, что распределение осциллирует с периодом 2 т . Значения параметров точно такие же, как на рис. 4
85
В заключение отметим следующее. В данной работе рассмотрено влияние внутреннего шума самих генетических реакций (называемых в зарубежной литературе «intrinsic noise») на структурообразование поля протеина, синтезируемое в клеточной ткани с некоторым запаздыванием. Так как пространственно-распределенная версия алгоритма Гиллеспи до сих пор в литературе не предложена, мы использовали клеточную модель, где обмен сигналами между клетками происходит на основании простых диффузионных соотношений. Важным новым результатом является демонстрация эффекта возбуждения подкритических колебаний взаимным действием шума и запаздывания в пространственном измерении. Показано, что в подкритической области в этом случае возможны не просто колебания, но и структурообразование в виде кластеризации групп клеток.
Список литературы
1. Брацун Д.А. Эффект возбуждения подкритических колебаний в стохастических системах с запаздыванием. Ч.1. Регуляция экспрессии генов // Компьютерные исследования и моделирование. - 2011. - Т. 3, № 4. - С. 431438.
2. Брацун Д.А., Захаров А.П. Моделирование пространственно-временной динамики циркадианных ритмов Neurospora crassa // Компьютерные исследования и моделирование. - 2011. - Т. 3, № 2. - С. 191-213.
3. Брацун Д.А., Захаров А.П., Письмен Л.М. Многоуровневое математическое моделирование возникновения и роста опухоли в ткани эпителия // Компьютерные исследования и моделирование. - 2014. - Т. 6, № 4. - С.585-604.
4. Брацун Д.А., Зюзгин А.В. Эффект возбуждения подкритических колебаний в стохастических системах с запаздыванием. 4.II. Управление равновесием жидкости // Компьютерные исследования и моделирование. -2012. - Т. 4, № 2. - С. 375-396.
5. Мюррей Дж. Математическая биология. Т. 1. - М: ИКИ-РХД, 2009. -774 с.
6. Об активном управлении равновесием жидкости в термосифоне // Д.А. Брацун, А.В. Зюзгин, К.В. Половинкин, Г.Ф. Путин. ПЖТФ. - 2008. - Т.34, вып.15. - С. 36-42.
7. Янушевский Р.Т. Управление объектами с запаздыванием. - М.: Наука, 1978. - 416 с.
8. Bratsun D.A. Effect of unsteady forces on the stability of non-isothermal particulate flow under finite-frequency vibrations // Microgravity Sci. Technol. -2009. - Vol. 21. - P. 153-158.
9. Delay-induced stochastic oscillations in gene regulation // D. Bratsun, D. Volfson, J. Hasty, L. Tsimring. PNAS. - 2005. - Vol. 102, № 41. - P. 1459314598.
10. Gillespie D. T. Exact stochastic simulation of coupled chemical reactions // J. Phys. Chem. - 1977. - Vol. 81. - P. 2340-2361.
11. Huber D., Tsimring L.S. Cooperative dynamics in a network of stochastic elements with delayed feedback // Phys. Rev. E. - 2005. - Vol. 71. - P. 036150036165.
12. Kepler T.B., Elston T.C. Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations // Biophys. J. - 2001. - Vol. 81, № 6. - P. 3116-3136.
13. Lemerle C., Di Ventura B., Serrano L. Space as the Final Frontier in Stochastic Simulations of Biological Systems // FEBS Letters. - 2005. - Vol. 579. - P. 1789-1794.
14. Li C.-W., Chen B.-S. Stochastic Spatio-Temporal Dynamic Model for Gene/Protein Interaction Network in Early Drosophila Development // Gene Regulation and Systems Biology. - 2009. - Vol. 3. - P. 191-210.
15. Losson J., Mackey M.C. Evolution of probability densities in stochastic coupled map lattices // Phys. Rev. E. - 1995. - Vol. 52. - P. 1403-1417.
16. Non-Markovian processes in Gene Regulation // D. Bratsun, D. Volfson, J. Hasty, L. Tsimring. Noise in Complex Systems and Stochastic Dynamics III / edited by Laszlo B. Kish, Katja Lindenberg, Zoltan Gingl, Proceeding of SPIE. -2005. - Vol. 5845. - P. 210-219.
17. Ohira T., Yamane T. Delayed stochastic systems // Phys. Rev. E. - 2000. -Vol. 61. - P. 1247-1257.
18. Pawlik A.H., Pikovsky A. Control of oscillators coherence by multiple delayed feedback // Physics Letters A. - 2006. - Vol. 358. - P. 181-185.
19. Salm M., Pismen L.M. Chemical and mechanical signaling in epithelial spreading // Phys. Biol. - 2012. - Vol. 9, №. 2. - P. 026009-026023.
20. Small delay approximation of stochastic delay differential equations // S. Guillouzic, I. L'Heureux, A. Longtin. Phys. Rev. E. - 1999. - Vol. 59. - P. 39703982.
21. Tsimring L.S. Noise in biology // Rep. Prog. Phys. - 2014. - Vol. 77. -P. 026601.
22. Tsimring L.S., Pikovsky A. Noise-induced dynamics in bistable systems with delay // Phys. Rev. Lett. - 2001. - Vol. 87. - P. 2506021-2506024.
