CC BY
90УДК 517.958.57
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЦИРКАДИАННЫХ РИТМОВ НА РАЗНЫХ МАСШТАБНЫХ УРОВНЯХ ОПИСАНИЯ БИОСИСТЕМЫ А.П. Захаров, Д.А. Брацун
Пермский государственный гуманитарно-педагогический университет, г.Пермь
MATHEMATICAL MODELING OF CIRCADIAN RHYTHMS AT DIFFERENT SCALES OF BIOSYSTEMS DESCRIPTION A. P. Zakharov, D. A. Bratsun
Perm State Humanitarian Pedagogical University, Perm
В работе выделяются три характерных масштаба описания биосистемы: микроскопический (размер гена), мезоскопический (размер клетки) и макроскопический (размер организма). Для каждого случая обсуждается подход к моделированию циркадианных ритмов на примере предложенной ранее модели с запаздыванием. На уровне гена использовалось стохастическое описание. Показана устойчивость механизма ритмов по отношению к флюктуациям. На мезоскопическом уровне предложено детерминистское описание в рамках пространственно-распределенной модели. Обнаружен эффект групповой синхронизации колебаний в клетках. Макроскопические эффекты исследованы в рамках дискретной модели, описывающей коллективное поведение большого числа клеток. Обсуждается вопрос о соединении результатов, полученных на разных уровнях описания.
Ключевые слова: биоритмы, синхронизация, системы с запаздыванием.
In the paper three characteristic scales of a biological system are proposed: microscopic (gene's size), mesoscopic (cell’s size) and macroscopic level (organism’s size). For each case the approach to modeling of circadian rhythms is discussed on the base of a time-delay model. At gene’s scale the stochastic description has been used. The robustness of rhythms mechanism to the fluctuations has been demonstrated. At the mesoscopic scale we propose the deterministic description within the spatially extended model. It was found the effect of collective synchronization of rhythms in cells. Macroscopic effects have been studied within the discrete model describing the collective behavior of large amount of cells. The problem of cross-linking of results obtained at different scales is discussed.
© Захаров А.П., Брацун Д.А., 2013
Key words: biorhythms, synchronization, time-delay system
Биологические ритмы - периодически повторяющиеся изменения интенсивности биологических процессов и явлений, которые свойственны живой материи на всех уровнях ее организации - от молекулярных и субклеточных до биосферы в целом [1]. Возникновение ритмов в живых организмах связано с процессами их адаптации к окружающей среде в процессе эволюционного развития. Стабильность периодичности изменения освещенности, температуры, влажности, геомагнитного поля и других параметров окружающей среды, обусловленных движением Земли и Луны вокруг Солнца, позволила живым системам в процессе эволюции выработать устойчивые к внешним воздействиям временные программы. Как сейчас уже стало понятно, такие ритмы закреплены в генетической структуре. Поэтому при отсутствии внешнего воздействия периоды таких ритмов отличаются от соответствующих в окружающей среде [21]. Хотя классификация биоритмов достаточно разнообразна [1, 5], в данной работе нас будут интересовать ритмы, синхронные с суточными колебаниями внешней среды, называемые циркадианными.
Во времена, когда генетический характер закрепления механизма циркадианных колебаний еще не был ясен, их приписывали определенному уровню существования организма. Например, в монографии [7] предложена классификация ритмов, в которой клеткам и тканям приписываются ультрадианные ритмы (короткие околочасовые колебания). Циркадианные ритмы формируются на уровне всего организма в целом. Согласно этой гипотезе, последние служат ритмами-водителями: они лабильны по отношению к действию внешних факторов, могут синхронизироваться с ними, а также воздействовать на подчиненные им ритмы-ведомые. Однако обнаружение генетического механизма воспроизводства циркадианных колебаний привело к пониманию, что эта картина неверна. Механизм околосуточных ритмов функционирует даже в масштабах одного или нескольких генов, проявляясь в колебаниях уровня мРНК в процессах транскрипции. Он же проявляет себя в процессах изменения концентраций белка в клетке после трансляции. Наконец, циркадианные колебания неизбежно должны возникать на клеточном уровне после того, как транспортные белки просачиваются через клеточные мембраны и вступают в межклеточное взаимодействие. На уровне органов сигналы от отдельных клеток должны синхронизироваться, создавая единый ритм для всего организма.
Отметим, что проблема пространственной синхронизации множества взаимодействующих осцилляторов весьма популярна у физиков [20], но в работах биологов рассматривается нечасто. Большинство источников по циркадианным
колебаниям концентрируется на вопросах устройства механизма временных колебаний.
Данная работа посвящена обсуждению подходов к моделированию циркадианных ритмов на разных пространственных масштабах живого организма, а также изучению форм синхронизации колебаний, формирующихся на разных уровнях функционирования системы.
Уровни описания биологической системы. Выделим три основных уровня описания биологической системы, основываясь на характерных масштабах структурных единиц этого организма (рис.1).
Как известно, клеточное ядро, в котором хранится генетическая информация, у эукариот имеет размер от 1 до 10 мкм, если не брать в расчет крупные ядра половых клеток, достигающих гигантских размеров (до 500 мкм). В ядре происходит процесс транскрипции, включающий в себя считывание информации с конкретных генов с помощью РНК-полимеразы и синтез мРНК. Затем молекулы мРНК покидают ядро и поступают в рибосому, где синтезируются белки. Рибосома имеет размер, на два порядка уступающий ядру, - около 0,03 мкм. Белки могут возвращаться в ядро и оказывать непосредственное влияние на процессы транскрипции, взаимодействуя с промоутерами генов. Следовательно, характерный масштаб пространственных процессов определяется здесь размерами ядра. С другой стороны, количество молекул, вступающих в реакции во время генных процессов транскрипции/трансляции, весьма невелико - скорость работы РНК-полимеразы, например, всего около 50 нуклеотидов в секунду. Таким образом, даже небольшие флюктуации концентраций рибонуклеиновой кислоты (РНК) и белка могут иметь значительное влияние на общую динамику системы. В этом состоит принципиальное отличие этих систем от, скажем, гидродинамических систем с шумом, в эволюции которых принимает участие настолько огромное количество молекул, что флюктуации могут проявить себя только вблизи точки бифуркации [2;18]. Таким образом, описание системы в масштабах одного или нескольких генов должно быть принципиально стохастическим. При этом пространственно-распределенной динамикой молекулярного облака, очевидно, здесь можно пренебречь, так как речь идет о весьма небольшом количестве элементов (десятки или сотни молекул), и приближение сплошной среды здесь не имеет смысла. Этот уровень описания назовем микроскопическим (см. рис.1).
/ /
Микроскопический / Мезоскопический / Макроскопический
уровень I уровень I уровень
1 I
>> V \ \ \ \ \ Длина (мкм)^>
1 10 1 I \ ю ! со _ о см _ о 104
V У V У У
Ядро клетки / Клетка / Организм
' /
Рис. 1. Характерные пространственные масштабы биологической системы
На следующем уровне описания биологической системы, который мы определим как мезоскопический, выделяется клетка или несколько клеток, обменивающихся сигналами. Типичный размер клетки эукариота порядка 10-100 мкм. Этот размер задает пространственный масштаб описания системы. Сколько молекул протеина того или иного типа имеется в клетке? Исследования, дающие точный ответ на этот вопрос, немногочисленны. В фундаментальной и единственной в своем роде работе [17] довольно подробно были исследованы дрожжи. В среднем для разных типов белков были получены четырехзначные или пятизначные числа. Например, ген /щ (внутренняя классификация УВК373Ш), отвечающий за поддержание механизма циркадианных ритмов в дрожжах и некоторых других грибах, обнаружен в количестве 7160 штук, приходящихся на одну клетку. Это означает, что можно пренебречь стохастическими свойствами системы, так ярко проявляющими себя на микроскопическом уровне описания, поскольку роль флюктуаций здесь не велика. Но, учитывая, что размер одной молекулы белка составляет порядка 0,001-0,01 мкм, можно ввести пространственно-распределенную модель описания динамики протеинового облака. Известно, что коэффициент диффузии разных белковых молекул в форме мономеров при невысокой их концентрации в водном растворе оценочно равен
7 2
10 см /с [8]. Это значение несколько ниже стандартного для диффузии веществ в воде, так как мономеры белка представляют собой крупные тяжелые молекулы, состоящие из сотен и тысяч аминокислотных остатков. Однако диффузия белка в цитоплазме живой клетки ограничена строением внутриклеточного пространства. Как сейчас уже стало понятно, цитоскелет клетки, состоящий из актиновых филоментов, представляет собой сложную вязкоупругую среду, допускающую при определенных условиях разжижение [10; 19]. Значение коэффициента диффузии белковых клеток в таком внутриклеточном каркасе в данный момент неизвестно, но оно точно меньше вышеприведенного значения
для воды. Вероятно, хорошей оценкой «снизу» для этой величины является
12 2
значение для латеральной диффузии белка в клеточных мембранах - 10 см /с
[6]. Таким образом, на мезоскопическом уровне описания мы предполагаем, что пространственно-распределенная детерминистская модель процессов реакции-диффузии является наиболее адекватной строению системы.
Наконец, макроскопический уровень предполагает описание ансамбля большого числа взаимодействующих клеток (см. рис.1). Это уровень органа или даже всего организма в целом. В отличие от мезоскопических масштабов, на которых внутриклеточное пространство, вообще говоря, открыто для перемещения молекул белка, здесь система опять становится дискретной, так как клетка является элементарной единицей строения организма, обладающей способностью к самостоятельному существованию, воспроизведению и развитию. Из-за действия межклеточных мембран характерная величина диффузии химических сигналов, которыми обмениваются клетки, достаточно мала. Это значит, что обмен между клетками происходит сравнительно медленно, и изучение динамики системы на этом уровне допускает использование такой дискретной модели, которая предполагает уровень концентрации белка в отдельных клетках, не зависящий от пространства. Данный подход позволяет абстрагироваться от мезоскопической динамики внутри клеток и учитывать только межклеточные различия.
Модель циркадианных колебаний. На всех трех уровнях описания будем использовать динамическую модель циркадианных колебаний, предложенную авторами ранее [4]. Несмотря на то что модель предлагалась для описания суточных ритмов организма Neurospora crassa, она имеет достаточно общий характер. К примеру, в работе [23] модель подобного класса легко обобщается на случай циркадианных колебаний у мушки Drosophila. Главным элементом механизма колебаний является эффект запаздывания реакций синтеза белков в процессах транскрипции и трансляции генов (рис. 2). Эти процессы не просто медленные, но еще и состоят из многоэтапных биохимических реакций, в ходе которых последовательно образуются сложные органические соединения. Таким образом, эти процессы растянуты во времени, а значит идут с некоторым характерным запаздыванием. Подробнее вопросы, связанные с функционированием механизма запаздывающих колебаний, обсуждаются в работах [2, 4].
На рис.2 представлена схема взаимодействия двух белков, формально обозначенных как F и W, кодируемых двумя генами f и w. У ряда организмов были выделены такие пары принципиальных генов, отвечающих за работу механизма циркадианных колебаний. Например, в случае нейроспоры это гены frq и w^, у дрозофилы - per и dclock [14]. Это не значит, что механизм колебаний поддерживается работой исключительно этой парой генов, так как циркадианные ритмы обладают широким набором свойств (автономность, способность компенсировать температурные изменения среды, а также изменять фазы
колебаний под действием внешнего освещения или температуры и т.д.), задаваемых десятками генов. Однако указанные пары являются принципиальными для поддержания циркадианных ритмов.
Рис. 2. Схема взаимодействий белков в модели циркадианных колебаний.
Здесь: £ w - гены; F, Ж - соответствующие белки
Как видно из рис. 2, обобщенная модель включает в себя как положительную, так и отрицательную петлю обратной связи и является симметричной по отношению к обеим переменным. Кроме того, она учитывает процессы димеризации и деградации белков. Полный список реакций, происходящих при транскрипции генов, приведен в табл. 1. Концентрации белков-мономеров в таблице обозначены F и ж, димеров - F2 иж2. Динамика оператор-сайта промотора описывается бинарной функцией в е {В0, В1}, которая принимает значение В0 в случае его открытия и В1 в случае закрытия. Наиболее важными для
поддержания колебаний являются реакции синтеза белков, протекающие с характерным временем запаздывания т. Именно этот механизм определяет динамику системы и её устойчивость к внешним и внутренним возмущениям.
Таблица 1
Список реакций транскрипции генов frq и wc
Процесс Реакция
1 Димеризации р + р кр > р2 , ж + ж кЖ > ж2
2 Редимеризации р2 к-1 > Р + Р, ж2 кЖ >ж + ж
3 Динамика оператор-сайтов Вр + ж2 Вр , В^ + Р2 В? Вр—к=2> В0р + Ж2 , Ц'г—кЖ> В^ + Р2
4 Синтеза белков Вр (0 —Вр + Р‘+Тр, В? ^) кж > в? + ж‘+тж
5 Распада белков Р 7р >0, ж 7ж >0
6 Образования гетеродимера р + ж—^> 0
Примечание: здесь, и далее к, к], к-1, к2, к-2, ур, уж, Ьр, Ьж - скорости соответствующих реакций.
На основании цепочки связанных биохимических реакций, представленных в табл. 1, может быть получена динамическая модель циркадианных ритмов (1-2) в детерминистском описании [4]:
(1 + 4 КрГ) § = V
1
1
1 + К. К Рж 2(t - т)
Ж ёЖ
(1 + 4К| Ж)
= кт
Ж
1 -
1
1 + Кр К Ж Г 2^ - т)
у Г - кГЖ, (1)
Г
- у.Ж - кГЖ, (2)
где К = к / кч, Кг. = к / к-г.. В численных расчетах были использованы значения параметров, приведенные в табл. 2, которые позволяют получить колебания с периодом 22,65 часа.
Таблица 2
Параметры модели
т к кр кж К К К. К Ж УР Уж
6 h 30 пМ"1^1 8 пМЖ 4 nM/h 5 пМ"1 5 пМ"1 5 пМ"1 5 пМ"1 0.3 К1 0.4 К1
Отметим, что математическая модель (1-2) не равнозначна совокупности реакций, представленных в табл. 1, так как она была выведена из предположения, что часть реакций являются быстрыми, а часть медленными. Таким образом, на фоне медленно меняющихся величин (например, общего количества молекул белка) реагенты, участвующие в быстрых реакциях, стремительно достигают состояния локального статистического равновесия. Таким образом, в зависимости от уровня описания системы и совокупности сделанных допущений необходимо пользоваться либо исходной системой кинетических уравнений химических реакций (табл.1), либо упрощенной детерминистской моделью (1-2).
Микроскопический уровень. Как отмечалось выше, на микроскопическом уровне (1 - 10 мкм, см. рис.1) мы можем пренебречь пространственно-распределенными эффектами и рассматривать процессы транскрипции-трансляции зависящими лишь от времени. Особенностью описания динамических явлений в масштабах гена является небольшое количество взаимодействующих молекул белка, поэтому многократно возрастает влияние флюктуаций на поведение системы. Следовательно, математическая модель на этом уровне описания может быть только стохастической.
В стохастических исследованиях генных процессов принято выделять два вида шума: внешний шум, который генерируется за пределами клетки и связан с межклеточными различиями, и внутренний шум, связанный со случайной
природой протекающих химических реакций, флюктуациями температуры и т.д. Описания динамики системы в масштабах гена предполагает влияние лишь внутреннего шума.
Удобным инструментом исследования стохастических генетических моделей является метод Гиллеспи [16], который представляет собой разновидность методов Монте-Карло. Особенно интенсивно этот метод используется в вычислительной системной биологии при изучении внутриклеточной динамики, в которой участвует сравнительно небольшое количество молекул. Как известно, численное решение, полученное методом Гиллеспи, статистически воспроизводит точное решение основного кинетического уравнения (мастер-уравнения). Важной особенностью метода является использование непосредственно кинетических уравнений химических реакций. Отметим, что классическая версия алгоритма Гиллеспи была сформулирована только для марковских систем. Поэтому в работе [11] была предложена модификация алгоритма на случай запаздывающих химических реакций. Там же на примерах показано, что новый алгоритм корректно воспроизводит их решения численно, а также приведены примеры расчета для тестовых задач. В работе [4] даны результаты расчетов стохастической модели циркадианных колебаний с запаздыванием для нейроспоры. Приведем здесь принципиальный результат этих расчетов.
то
60
О » О
£ 40
п 30
О
£ м
и,
‘—' ю
о
а
Йте (Поигз)
б
Рис.3. Временная зависимость полного числа мономеров белка Г (сплошная линия) и Ж (штриховая линия), полученная при (а) численном интегрировании системы (1-2).
(б) Временная эволюция соответствующей стохастической системы, полученная с помощью алгоритма Гиллеспи, модифицированного для немарковских процессов.
Во врезках приведены фурье-спектры сигналов. Значения параметров системы в обоих
случаях одинаковы (см. табл.2)
На рис.За представлены развертывания по времени для полного числа мономеров обоего сорта, полученные при численном расчете детерминистской динамической системы (1-2) при отсутствии базальных транскрипционных факторов. На рис.Зб представлены результаты стохастических расчетов для тех же значений параметров, что и в детерминистском случае. Хорошо видно, что,
несмотря на небольшое количество вовлеченных в динамику молекул (в любой момент времени полное число мономеров не превышает 50 штук) и значительные флюктуации, базисный механизм колебаний уверенно работает. Он основывается на синхронизации колебаний концентраций белков Г и Ж в противофазе. Так как концентрации белков соединены друг с другом посредством положительной обратной связи, то последний вывод представляется вполне естественным. Можно заметить, что полной симметрии между колебаниями концентраций не наблюдается, так как несимметричен набор значений параметров (см. табл. 2). Стохастичность системы проявляет себя во флюктуациях, которые достигают 20-40 % от амплитуд колебаний в детерминистском случае. Во врезках к рисункам приведены фурье-спектры сигналов. Можно заметить, что стохастический сигнал имеет устойчивый максимум, примерно соответствующий периоду колебаний 2З часа. Подобное поведение системы еще раз подчеркивает устойчивость базисного механизма генерирования циркадианных колебаний, основанного на запаздывании в процессах транскрипции и трансляции.
Мезоскопический уровень. Моделирование синхронизации биоритмов в масштабах одной или нескольких клеток (10-100 мкм, см. рис.1) предполагает построение пространственно-распределенной модели. Для этого модель циркадианных колебаний (1-2) была дополнена диффузионным слагаемым, где в - коэффициент диффузии молекул белка внутри клетки.
а Г _ 1
^7 _ (1 + 4 КрГ)
дЖ 1
дг (1 + 4 к'Жж)
kг
\
k-----------Г---------у Г - kFW
Г 1+к'жк:ж'-(г - т) Ж
+ в
К
У V
Л
Ад2Г д2ГЛ дх^+ ду2
1 + ЦК ЖГ 2(г - т)
+ в
^д2Г д2ГЛ +
дх2 ду2
(3)
(4)
Уравнения (3-4) сформулированы для переменных Г (г, х, у) и Ж (г, х, у), которые теперь являются функциями времени и двух пространственных координат х, у. В данной статье мы ограничиваемся рассмотрением пространственно-временной динамики системы в двух пространственных измерениях.
Численное исследование системы (3-4) осуществлялось с помощью метода конечных разностей, который является одним из наиболее распространенных численных методов решения уравнений в частных производных. Однако использование метода сеток, как наиболее простого и надежного, для решения задач с запаздыванием влечет за собой определенные вычислительные трудности. На каждом шаге по времени требуется информация о значениях
функций в прошлом, приводит к необходимости хранения большого объема данных внутри диапазона времени запаздывания. Ситуация усугубляется тем, что для явной схемы очень важно контролировать шаг по времени - для устойчивости расчета он должен быть меньше определенного значения (критерий Куранта):
Для преодоления этих трудностей в работе [3] был предложен специальный алгоритм оптимального хранения полей, при котором запоминание данных осуществляется не на каждом шаге по времени, а только на некоторых из них, называемых опорными. При этом сами значения полей в эти моменты времени назовем опорными временными слоями.
Значения переменных для промежуточных временных слоев восстанавливаются путем интерполяции данных опорных временных слоёв, хранимых в памяти (рис. 4). Предполагается, что в зависимости от параметров конкретной модели, гладкости функций и размерности накладываемой сетки могут быть использованы различные виды интерполяции. Для повышения эффективности работы алгоритма предусмотрена адаптивная схема использования количества опорных слоев, то есть возможно, что в процессе расчета плотность расположения опорных слоев в пределах диапазона запаздывания может меняться. Частота выделения слоёв определяется гладкостью функции. Во-первых, это позволяет фиксировать точность вычислений, а во-вторых, снизить методические погрешности, связанные с динамически меняющимся шагом по времени. Подробнее численный метод изложен авторами в [3]. В работе показано, что адаптивная схема для рассматриваемой задачи позволяет примерно на порядок ускорить время расчета по сравнению с таковой, где запоминаются данные всех полей. При этом расчет динамики режима бегущих волн, весьма чувствительных к погрешности схемы, не выявил существенных отклонений от результатов тестовой задачи.
Приведем некоторые результаты численных расчетов, касающиеся пространственно-временной динамики белковых полей. При случайном распределении начальных возмущений в отдельных центрах зарождаются спиральные волны. Эти волны продвигаются во всех направлениях пространства, пока не захватывают всю область интегрирования
Рис. 4. Зависимость количества опорных слоев (обозначены как BS)
от гладкости функции
Внутри области неустойчивости возникают многочисленные центры возбуждения вторичных спиральных волн, нелинейное взаимодействие которых приводит к хаотическому поведению среды за гребнем первичной волны (рис. 5). Расчёт взаимно-корреляционной функции между парами точек среды, принадлежащих областям до и после прохождения волны, показывает, что корреляции между точками в волновой зоне практически отсутствуют [4]. Таким образом, обнаруженные волны можно трактовать как волны рассинхронизации, которые сбрасывают все локальные установки фаз часов в среде. Тем не менее поведение системы нельзя назвать пассивным - она, скорее, находится в «ждущем» режиме по отношению к внешней среде. Хотя поле концентрации белка ¥ (с белком ж происходит то же самое) хаотически осциллирует как во времени, так и в пространстве, но при определенной команде извне система готова быстро синхронизировать свои колебания в соответствии с внешним сигналом.
т * < Г7> ->
■ V -Л о
•л,], ч-<^* 4 Г -V Ч'^ \
1|
Рис. 5. Характерный вид волновых структур, образуемых белком ¥ во внутриклеточном пространстве из случайных начальных условий в условиях полной темноты. .0=0,01
Отметим, что многие биологические процессы, начиная с генной экспрессии, клеточной пролиферации и вплоть до высокоуровневых процессов, таких как зрение, память и обучение, предполагают наличие обмена информацией с окружающей средой. Эти процессы осуществляются, в том числе, посредством передачи сигнала через клеточную мембрану или протеиновые рецепторы (каналы), которые соединяют клетку с окружающим пространством. Мембраны не только выделяют различные реакционные объемы живых клеток, отделяя живое от среды. Они играют ключевую роль в метаболизме, селективно пропуская потоки неорганических ионов и органических молекул. Известно, что изменение уровня этих потоков является одним из механизмов влияния на химический состав и химические структуры внутри клетки.
Рис.6. Схема расчетной области
Моделирование данного процесса выполнено следующим образом: в расчетной области (рис.6) выделяются определенные участки, которые мы будем условно называть «клетками», в которых протекают реакции (3-4), задающие циркадианные ритмы. Так как механика самих клеток в данный момент нас не интересует, то модельные клетки имеют форму квадрата и фиксированный размер. Граница клеток не однородна: большая её часть непроницаема для реагента (отмечена на рис.6 черной линией), а оставшаяся часть имеет диффузионную проницаемость, вычисляемую по конечноразностной формуле (6) [12,13],
¥,( х,0) = ¥2( х,0) =
0(4¥1 (х,-И) - ¥ (х,-2И)) + 00 (4¥2 (х, И) - ¥2 (х,2И))
(6)
3( О + А)
где h - шаг по сетке, О0 - коэффициент диффузии белка во внеклеточном пространстве. Таким образом, каждый участок, в котором протекает реакция, является моделью живой клетки с собственным циркадианным ритмом внутри. Каждая клетка имеет четыре мембраны фиксированной ширины, через которые
она может обмениваться с межклеточным пространством. Для простоты мы не стали вводить в модель отдельный вид белка, ответственный за транспортировку циркадианного сигнала наружу. Транспортную функцию в предлагаемой модели выполняет всё тот же белок ¥, важный для поддержания циркадианных ритмов. Его партнер, белок Ж, может перемещаться только внутри пространства клеток - для него клеточные мембраны закрыты.
Предполагается, что в межклеточном пространстве транспортный белок не вступает в реакцию: он может лишь распространяться там благодаря диффузии:
дР = D дt 0
ґд2F д2р' д^ + ду2
(7)
у
Попадая внутрь другой клетки извне, белок ¥ вновь может участвовать в процессах транскрипции и трансляции. В этом случае он выступает как базальный фактор транскрипции по отношению к новой клетке.
Перейдем к результатам моделирования системы с клетками. На рис. 7 представлены четыре кадра динамики изменения концентрации белка ¥ в полной темноте и в отсутствие факторов базальной транскрипции, когда только одна из клеток (#3) генерирует циркадианный ритм. В остальных трёх клетках механизм транскрипции выключен: k¥ = 0. Как было установлено
выше (см. рис. 5), система в этом случае находится в состоянии фазовой турбулентности, которое в отсутствии внешнего сигнала может сохраняться сколь угодно долго. Несмотря на то что существует поток белка ¥ через мембраны клетки #3 в межклеточное пространство, клетка продолжает функционировать в прежнем режиме хаотических пространственных волн. Убыль белка компенсируется его синтезом внутри клетки. Попадая в межклеточное пространство, выполняющий транспортные функции белок ¥ либо диффундирует там, либо проникает в клетки #1, #2 и #4. Механизм циркадианных колебаний в этих клетках парализован, однако процессы линейной и нелинейной деградации белка действуют. Поэтому концентрация ¥ внутри выключенных клеток быстро падает до нуля (рис. 7). Представленный на рисунке сценарий не зависит от размеров клеток и ширины межклеточного зазора между ними. Синхронизация колебаний в клетке #3 на рис.7 не наступает из-за отсутствия внешнего сигнала.
Рис. 7. Динамика концентрации белка ¥ в полной темноте. Только в одной клетке (#3) происходит реакция синтезирования белка. В остальных клетках (#1, #2 и #4) белок только деградирует. Последовательно представлены моменты времени /=200,
5000, 10000, 20000. Коэффициент диффузии в клетках равен 0=0,01, в межклеточном пространстве 0=0,01.
Качественное изменение в системе возникает, когда при тех же условиях начинают функционировать все четыре клетки (рис.8). Выявлено, что спонтанная пространственная синхронизация возникает, когда концентрация белка в межклеточном пространстве достигает определенного уровня концентрации ¥ = 0,8, примерно соответствующей среднему значению, около которого осциллирует поле концентрации в клетках. Именно общее для всех клеток межклеточное диффузионное пространство настраивает систему таким образом, что колебания приобретают упорядоченный синхронизированный вид. Отметим, что процесс синхронизации занимает значительное время / * 10000, особенно в сравнении с периодом базового колебания.
Рис. 8. Динамика концентрации белка ¥ в полной темноте. Все четыре клетки функционируют, вырабатывая белок ¥. Последовательно представлены моменты времени /=200, 5000, 10000, 20000. Коэффициент диффузии в клетках равен 0=0,01,
в межклеточном пространстве 0=0,01
Рис. 9. Динамика концентрации белка ¥ в системе из 16 плотно примыкающих друг к другу клеток с разными фазами циркадианного ритма в полной темноте.
Последовательно представлены моменты времени /=0, 2000, 4000, 6000, 8000, 10000.
Коэффициент диффузии в клетках равен 0=0,01.
Численные эксперименты показали, что важную роль в синхронизации играет ширина межклеточного зазора. Чем шире зазор, тем больше времени требуется для достижения состояния синхронных колебаний. Межклеточное диффузионное пространство играет здесь точно такую же роль, какую выполняют пружинки связи в системе связанных механических осцилляторов [20]. Чем сильнее связь между осцилляторами, тем легче достигается синхронизация. На рис.9 представлены кадры эволюции 16 клеток, которые настолько тесно расположены друг к другу, что межклеточное пространство здесь практически отсутствует. На самом деле такая конфигурация системы более реалистична, так как клетки тканей органов обычно образуют плотную сплошную среду без зазоров. На рис.9 клетки продолжают обмениваться друг с другом сигналами посредством того же самого диффузионного механизма, что и раньше. В начальный момент эволюции в каждой клетке задавалось постоянное для этой клетки значение концентрации обоих взаимодействующих белков, однако сами значения в разных клетках выбирались случайным образом. Так как уровень концентрации белков задаёт фазу циркадианного ритма (максимальный уровень ¥ определяет локально «полдень»), то система интегрировалась со случайным значением фазы колебаний в каждой клетке. При этих начальных условиях спиральные волновые структуры, приведенные на рис.5, в системе не возникают. Для этого нужны резкие локальные флюктуации концентрации белка.
В случае, представленном на рис. 9, реализуется решение в виде стоячей волны. Хорошо видно, что сравнительно быстро фаза всей группы клеток выравнивается, - происходит фазовая синхронизация. Начиная примерно с момента времени / = 4000, система осциллирует с общей фазой. Отметим, что синхронизация здесь наступает гораздо раньше, чем в случае, показанном на рис.8.
Для изучения динамики системы было рассмотрено несколько вариантов начальных условий. Например, устанавливались разными не только фазы, но и частоты колебаний в клетках. Это достигается варьированием параметров к¥ и кж, ответственных за скорость протекания реакций синтезирования
белков. Результаты моделирования системы показывают, что в конечном итоге в пространстве клеток в этом случае также устанавливается единая фаза и единая частота колебаний.
Макроскопический уровень. Исследование синхронизации биоритмов в масштабах органа или даже целого организма предполагает описание функционирования и взаимодействия большого числа клеток. В недавней работе [22] была предложена единственная в своем роде хемомеханическая модель растущей ткани эпителия. Эпителиальная ткань представляет собой слой клеток, выстилающий поверхность (полость) органа или тела. Таким образом, при моделировании поведения эпителия можно ограничиться квазидвумерной системой, что существенно упрощает вычисления.
Рис.10. Элементы хемомеханической модели эпителия
Модель, предложенная в статье [22], включает расчет динамики отдельных клеток, представленных в виде многоугольников с разным числом вершин (система откалибрована так, что наиболее вероятной формой клетки является гексагональная ячейка, хотя появление других видов многоугольников также возможно). Клетки плотно примыкают друг к другу, образуя сплошную двумерную поверхность эпителия (рис.10). Модель обладает целым набором свойств, позволяющих удачно имитировать поведение реальной эпителиальной ткани:
- возможность изменения размеров клеток в процессе эволюции ткани (например, затягивание раны) и изменения локальных механических свойств среды;
- возможность роста общего количества клеток в системе посредством их деления в определённых условиях эволюции;
- возможность перемещения клеток в общей массе эпителия посредством механизма интеркаляции;
- расчет динамики концентрации веществ, участвующих в регуляции жизнедеятельности ткани, для каждой клетки сообщества;
- обмен химическими сигналами, осуществляемый между соседними клетками эпителия через общую границу (см. рис.10);
- учет эффекта поляризации клеток, которая происходит спонтанно или под воздействием внешних условий.
Таким образом, каждая клетка в модели испытывает ряд хемомеханических воздействий, под влиянием которых она эволюционирует вместе со всей системой. Так как внутри каждой клетки все поля зависят только от времени, то мелкая структура пространственно-распределенных эффектов, связанная с неоднородностью полей внутри клеток, в модели не определяется. Однако на расстояниях, больших по сравнению с размерами одной клетки, пространственное структурообразование проявляет себя в полной мере. Разработанную модель можно классифицировать как дискретную сложную систему, демонстрирующую коллективные эффекты, с индивидуальным поведением отдельных элементов системы. Познакомиться подробнее с деталями реализации хемомеханической модели растущего эпителия можно в соответствующей публикации [22] .
Модель эпителия мы дополнили расчетом циркадианных ритмов в каждой клетке сообщества согласно уравнениям (8, 9),
dFi 1
dt (1 + 4K[Fl)
kF-----------^----------yFF. — kF Ш
F 1 + К?К%Ж, 2(і — т) F 1 11
k
V " 1 ^1^2" і V- V у
dW 1 (, К ,
kW-------------,------У ш Ш. — kF Ш .
V ш 1 + к;к2^ (t — т) г гу
dt (1 + 4 К'шшг)
+ “Е (Fj — F,), (8)
і
і = , (9)
где а = 0,1 - коэффициент переноса белка ¥ через клеточную мембрану.
Слагаемое, ответственное за перенос белка из клетки в клетку, в уравнении (8) записано в виде простого разностного соотношения, в котором Ь.. имеет
смысл длины границы между .-й и .-й клетками (рис. 10). Суммирование в формуле (8) ведется только по клеткам, которые являются соседними для -й клетки. Если уровень концентрации белка ¥ в рассматриваемой клетке выше, чем у её окружения, то формируется поток молекул белка наружу. В обратном случае наблюдается приток белка извне. Транспортную функцию, как и прежде, выполняет только белок ¥, а для Ж границы клеток остаются
** Мы благодарны разработчикам модели за то, что они согласились передать программный комплекс авторам данной работы для исследований в области сложных живых систем.
непроницаемыми. В результате деления какой-то клетки на две части считается, что новые клетки наследуют фазу циркадианного ритма материнской клетки. Так как в данной статье нас не интересует влияние циркадианных ритмов на хемомеханические свойства среды, обратная связь уравнений (8, 9) с общим функционированием клеточной модели отсутствует.
Несмотря на то что уравнения (8, 9), в отличие от (3, 4), представляют собой систему обыкновенных дифференциальных уравнений, общее их количество достаточно велико, так как N (количество клеток) может достигать нескольких тысяч. Кроме того, при расчете уравнений необходимо проводить процедуру запоминания полей концентраций в пределах диапазона запаздывания и делать это для каждой клетки отдельно.
Перейдём к изложению результатов моделирования. На рис. 11 представлены несколько кадров эволюции эпителия, состоящего в начальный момент из 1600 клеток. В качестве начального условия в клетках случайным образом задавалась фаза циркадианного ритма, а внешние воздействия на систему исключались. Такой подход позволяет определить возможные собственные формы коллективного поведения клеток по синхронизации колебаний. Результаты расчетов показывают, что в пространстве большого количества клеток полной синхронизации ритмов в смысле абсолютного выравнивания фазы колебаний по всем клеткам, как это происходило в случае небольшого сообщества (см. рис.9), не получается. Здесь проявляется макроскопический эффект кластеризации колебаний - клетки формируют два примерно равных сообщества, которые коллективно осциллируют в противофазе (см. рис.11). Между двумя группами есть небольшая прослойка клеток, в которых реализуются колебания с промежуточными значениями фазы.
Кластеризация в системах с большим количеством элементов, обменивающихся химическими сигналами, с некоторых пор привлекает к себе внимание исследователей. Например, в недавней работе [14] подробно изучена растущая группа взаимодействующих друг с другом синтетических генетических осцилляторов. Обнаружено, что с течением времени происходит кластеризация ткани на два типа осциллирующих клеток. Авторы связывают это явление с двумя возможными устойчивыми состояниями равновесия у системы. Отмечается, что кластеризация, по-видимому, является важнейшей характерной особенностью больших сообществ и может служить причиной дальнейшей дифференциации клеток в органах.
Рис. 11. Эффект кластеризации циркадианных ритмов в эпителиальной ткани: колебания концентрации белка ¥ в условиях полной темноты в ткани эпителия, состоящего более чем из 1600 клеток. В качестве начальных условий использовано случайное распределение фазы колебаний в клетках. Последовательно представлены
моменты времени 1=0, 25, 50, 75, 100, 125
Заметим, что в нашей модели клетки могут перемещаться по пространству эпителия. В связи с этим интересно отметить, что переход клетки за границу кластера не сопровождается сохранением у неё прежней фазы колебаний: попадая в новое для себя окружение, клетка подстраивается под общую для этого кластера фазу колебаний. Таким образом, наблюдается осциллирующая на месте стоячая волна поля концентрации белка ¥ (см. рис. 11).
Был рассмотрен также случай специальных начальных условий, представляющий собой небольшое малое возмущение концентрации протеина ¥ на фоне нулевых значений полей в остальных точках области. Физически это означает, что циркадианный ритм всей ткани полностью синхронизирован по фазе, за исключением одной клетке в центре сообщества. Такие начальные условия являются гораздо более реалистичными, так как в тканях отдельные клетки постоянно сигнализируют сообществу о происходящих в них изменениях. Как эволюционирует система в этом случае изображено на рис. 12.
Рис. 12. Эффект бегущей волны циркадианных ритмов в эпителиальной ткани: колебания концентрации белка ¥ в условиях полной темноты в ткани эпителия, состоящего более чем из 1600 клеток. В качестве начальных условий использовано возмущение концентрации белка в одной клетке в центре сообщества. Последовательно представлены моменты времени /=0, 20, 40, 60, 80, 100
Из рис. 12 видно, что через пространство эпителия распространяются концентрические бегущие волны, имеющие центр в точке возмущения.
В системе устанавливается состояние динамической синхронизации процессов регуляции циркадианных ритмов.
Вопрос о соединении результатов моделирования, полученных на разных уровнях описания системы, является важным для понимания того, как интерпретировать результаты. Мы рассмотрели динамику циркадианных колебаний в масштабах гена (стохастическая модель), клетки (детерминистская пространственно-распределённая модель с диффузией) и организма (дискретная хемомеханическая модель клеток эпителия). На всех уровнях описания в модели был заложен генный механизм генерирования ритмов на основе бинарного взаимодействия двух белков (см. рис. 2). Соответственно базовые колебания с периодом примерно 23 часа проявляют себя во всех моделях. Было показано, что механизм колебаний устойчив по отношению к внутреннему шуму, связанному с кинетикой реакций (см. рис. 3). При отсутствии сильных возмущений в масштабах клетки колебания быстро синхронизируют свою фазу по всему внутриклеточному пространству. В случае взаимодействия нескольких клеток обмен диффузионными потоками приводит к синхронизации ритмов уже у всей группы клеток (см. рис. 9).
Если внутри клетки есть значительные локальные отклонения концентрации, то система может перейти к новому режиму осцилляций, представляющему собой хаотическое поле бегущих волн (см. рис.5). Характерный размер пространственных структур, возникающих в этом случае, примерно на порядок меньше размеров самой клетки. Сценарий синхронизации ритмов таких возбужденных клеток, взаимодействующих друг с другом, несколько отличается от прямолинейного выравнивания фазы. За счет транспортного белка в межклеточном пространстве устанавливается общее для всех клеток среднее поле, которое приводит к регуляризации бегущих волн и пространственно-временной синхронизации их распространения (см. рис. 8). Важность взаимодействия между клетками легко установить, если отключить процессы синтеза циркадианных белков во всех клетках, за исключением одной (см. рис. 7). В этом случае состояние хаотически распространяющихся бегущих волн может поддерживаться сколь угодно долго.
Любопытная картина складывается, когда мы поднимаемся к уровню описания групповой динамики большого числа клеток. Во-первых, здесь сразу же пропадают мелкие волновые структуры, возникающие в одиночных клетках (приведенные, например, на рис. 8). На этом уровне описания, минимальное пространственное разрешение системы определяется размером одной клетки. Не удивительно поэтому, что более мелкие структуры (а волны на рис. 8 на порядок мельче) сразу пропадают. Во-вторых, здесь возникают крупномасштабные структуры, состоящие из сотен и тысяч клеток.
Мы обнаружили эффект разбиения пространства эпителия по крайней мере на два домена сложной формы, внутри каждого из которых происходит выравнивание фаз колебаний (см. рис.11). В некоторых работах кластеризация рассматривается как важнейший начальный этап дифференциации клеток.
Наконец, при локальном выведении из равновесия фазы циркадианного ритма эпителия на общем синхронизированном поле обнаружено распространение глобальных концентрических волн, информирующих сообщество клеток о внезапной активизации одного из своих членов (см. рис. 12). Это позволяет клеткам всё время подстраивать свою фазу колебаний, а всему сообществу - быть готовым коллективным образом реагировать на изменение внешних условий.
В научной литературе существует ряд экспериментальных свидетельств, что процессы транскрипции/трансляции белков происходят с запаздыванием. Например, в работе [15] для организма Neurospora crassa было установлено, что такая запаздывающая обратная связь формируется у пары белков FRQ и WCC, синтезируемых генами frq и wcc. Что особенно важно, именно эта пара белков ответственна за поддержание циркадианных ритмов у нейроспоры. В данный момент трудно сказать, является ли запаздывание фундаментальным механизмом, присущим транскрипции любых генов. Однако то, что это важное свойство хотя бы некоторых из них, установлено достоверно.
Современные технологии позволяют экспериментаторам вместо изучения уже существующих организмов работать с искусственно созданными цепочками ДНК. Именно таким путем пошли авторы исследования [9]. Ими была создана синтетическая генная цепочка, которая включала в себя осциллирующую часть, основанную на паре генов LuxI и LuxR. Колебания уровня соответствующих белков достигались за счет запаздывания обратной связи. Кроме того, цепочка включала в себя ген, синтезирующий транспортный белок AHL (N-ацил-гомосерин лактон/ Для искусственной визуализации групповой динамики в каждую цепочку был добавлен ген yemGFP, вызывающий свечение под действием транспортного белка AHL. Синтетическая цепочка была интегрирована в каждую бактерию кишечной палочки, а вся колония помещена в специальную среду с питательным раствором. Таким образом, внедренные осцилляторы имитировали в каждом организме работу его биологических часов, а белок AHL осуществлял глобальную координацию поведения всей колонии. Численность колонии достигала несколько тысяч особей. В результате наблюдения за колонией бактерий было обнаружено спонтанное возникновение бегущих концентрических волн свечения. В сопроводительных файлах на сайте журнала Nature можно ознакомиться с различными видеозаписями экспериментов.
Отметим, что концентрические волны, распространяющиеся в хемомеханической модели эпителия (см. рис. 12), в деталях воспроизводят пространственно-временную динамику в колонии сообщающихся друг с другом бактерий.
Выводы. На основе предложенной ранее модели циркадианных колебаний с запаздыванием проведено исследование синхронизации циркадианных ритмов на различных уровнях описания биологической системы. Изучение взаимодействия в масштабах нескольких генов осуществлено с помощью стохастического моделирования. Пространственно-временная динамика поля концентрации белков в масштабах одной или нескольких клеток исследованы в рамках детерминистской модели с диффузией. Наконец, макроскопические эффекты, возникающие в сообществе большого числа клеток, рассмотрены на базе дискретной хемомеханической модели растущего эпителия, клетки которого могут обмениваться химическими сигналами. Выявлены формы синхронизации циркадианных ритмов на каждом уровне описания системы.
Список литературы
1. Ашофф Ю. Биологические ритмы. - М.: Мир, 1984. - С. 412.
2. Брацун Д.А. Эффект возбуждения подкритических колебаний в стохастических системах с запаздыванием. Ч. I. Регуляция экспрессии генов // Компьютерные исследования и моделирование. - 2011. - Т.3, № 4. -С.431-438.
3. Брацун Д.А., Захаров А.П. К вопросу о численном расчете пространственно-распределенных динамических систем с запаздыванием по времени // Вестн. Перм. ун-та. Сер. Математика. Механика. Информатика. - 2012. - Вып. 4(12). - С. 32-41.
4. Брацун Д.А., Захаров А.П. Моделирование пространственно-временной динамики циркадианных ритмов Neurospora crassa // Компьютерные исследования и моделирование. - 2011. - Т.3, № 2. - С.191-213.
5. Моисеева Н. И., Сысуев В. М. Временная среда и биологические ритмы. -Л.: Наука, 1981. - С. 127.
6. СеверинаЕ.С. Биохимия: учеб. для вузов. - М.: ГЭОТАР-Медиа, 2003. - С. 784.
7. Степанова С.И. Биоритмологические аспекты проблем адаптации. - М.: Наука, 1986. - С. 244.
8. Якубке Х.-Д., Ешкайт Х. Аминокислоты, пептиды, белки. - М.: Мир, 1985. -С. 457.
9. A synchronized quorum of genetic clocks / Danino T., Mondragon-Palomino O., Tsimring L., Hasty J. // Nature. - 2010. - №.463. - P.326-330.
10. An active biopolymer network controlled by molecular motors / Koenderink et al. // PNAS. - 2009. - Vol.106, №36. - P.15192-15197.
11. Bratsun D., Volfson D., Hasty J., Delay-induced stochastic oscillations in gene regulation // PNAS. - 2005. - Vol.102, №41. - P.14593-14598.
12. Bratsun D.A., De Wit A. Buoyancy-driven pattern formation in reactive immiscible two-layer systems. // Chem. Eng. Sci. - 2011. - Vol.66, №22. -P.5723-5734.
13. Bratsun D.A., De Wit A. Control of chemoconvective structures in a slab reactor // Techn. Phys. - 2008. - Vol.53. - P.146-153.
14. Cooperative differentiation through clustering in multicellular populations / A. Koseska, E. Ullner, E. Volkov, J. Kurths, J. Garcia-Ojalvo // J. Theor. Biol. -2010. -Vol. 263. - P.189-202.
15. Denault D.L., Loros J.J., Dunlap J.C. WC-2 mediates WC-1-FRQ interaction within the PAS protein-linked circadian feedback loop of Neurospora // EMBO J.
- 2001. - Vol.20. - P. 109-117.
16. Gillespie D.T. Exact stochastic simulation of coupled chemical reactions // J. Phys. Chem. - 1977. - Vol.81. - P.2340-2361.
17. Global analysis of protein expression in yeast / S. Ghaemmaghami, H. Won-Ki, K. Bower, R.W. Howson, A. Belle, N. Dephoure, E.K. O'Shea, J.S. Weissman // Nature. - 2003. - Vol.425. - P. 737-741.
18. Hasty J., Collins J. J. Translating the noise // Nature Genetics. - 2002. - Vol. 31.
- P. 13-14.
19. Morozov K.I., Pismen L.M. Cytoskeleton fluidization versus resolidification: Prestress effect // Phys. Rev. E. - 2011. - Vol.83. - P.051920-051928.
20. Pikovsky A., Rosenblum M., Kurths J. Synchronization - a universal concept in nonlinear sciences. Cambridge University Press, 2001. - P.432.
21. Pittendrigh C.S. Temporal organization: reflections of a Darwinian clock-watcher // Annu. Rev. Physiol. - 1993. - Vol.55. - P. 16-54.
22. Salm M., Pismen L.M. Chemical and mechanical signaling in epithelial spreading // Phys. Biol. - 2012. - Vol.9. - P.026009-026023.
23. Smolen P., Baxter D. A., Byrne J. H. Reduced models of the circadian oscillators in Neurospora crassa and Drosophila melanogaster illustrate mechanistic similarities // OMICS. - 2003. - Vol.7, №4. - P.337-354.
