Моделирование биоритмов в масштабах гена, клетки и всего организма Текст научной статьи по специальности «Математика»

2013

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

Математика. Механика. Информатика

Вып. 3(22)

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

УДК 517.958.57

Моделирование биоритмов в масштабах гена, клетки и всего организма

А. П. Захаров, Д. А. Брацун

Пермский государственный гуманитарно-педагогический университет Россия, 614600, Пермь, ул. Сибирская, 24 az1211@mail.ru; +79222465522

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

Ключевые слова: биоритмы; синхронизация; системы с запаздыванием.

Введение

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

© Захаров А. П., Брацун Д. А., 2013

Работа выполнена при финансовой поддержке Министерства образования Пермского края (грант С-26/244), Министерства науки и образования РФ и ФЦП (ГК № 14.740.11.0352), проекта №031-Ф Программы стратегического развития ПГГПУ.

бильность периодичности изменения освещенности, температуры, влажности, геомагнитного поля и других параметров окружающей среды, обусловленных движением Земли и Луны вокруг Солнца, позволила живым системам в процессе эволюции выработать устойчивые к внешним воздействиям временные программы. Как сейчас уже стало понятно, такие ритмы закреплены в генетической структуре. Поэтому в отсутствие внешнего воздействия периоды таких ритмов отличаются от периодов соответствующих ритмов окружающей среды [2]. Хотя классификация биоритмов достаточно разнообразна [1; 3], в данной работе нас будут интересовать ритмы, синхронные с суточными колебаниями внешней среды, называемые циркадианными.

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

Отметим, что проблема пространственной синхронизации множества взаимодействующих осцилляторов весьма популярна у физиков [5], но в работах биологов встречается нечасто. Большинство источников по цир-кадианным колебаниям концентрируется на вопросах устройства механизма временных колебаний.

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

1. Уровни описания биологической системы

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

Как известно, клеточное ядро, в котором хранится генетическая информация, у эу-кариот имеет размер от 1 до 10 мкм, если не брать в расчет крупные ядра половых клеток, достигающих гигантских размеров (до 500 мкм). В ядре происходит процесс транскрипции, включающий в себя считывание инфор -мации с конкретных генов с помощью РНК-полимеразы и синтез мРНК. Затем молекулы мРНК покидают ядро и поступают в рибосому, где синтезируются белки. Рибосома имеет размер, на два порядка уступающий ядру, -около 0.03 мкм. Белки могут возвращаться в ядро и оказывать непосредственное влияние на процессы транскрипции, взаимодействуя с промоутерами генов. Следовательно, характерный масштаб пространственных процессов определяется здесь размерами ядра. С другой стороны, количество молекул, вступающих в реакции во время генных процессов транскрипции/трансляции, весьма невелико - скорость работы РНК-полимеразы, например, всего около 50 нуклеотидов в секунду. Таким образом, даже небольшие флуктуации концентраций рибону-

Рис. 1. Характерные пространственные масштабы биологической системы

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

На следующем уровне описания биологической системы, который мы определим как мезоскопический, выделяется клетка или несколько клеток, обменивающихся сигналами. Типичный размер клетки эукариота порядка 10-100 мкм (рис. 1). Этот размер задает пространственный масштаб описания системы. Сколько молекул протеина того или иного типа имеется в клетке? Исследования, дающие точный ответ на этот вопрос, немногочисленны. В фундаментальной и единственной в своем роде работе [8] довольно подробно были исследованы дрожжи. В среднем для разных типов белков были получены четырехзначные или пятизначные числа. Например, ген frq (внутренняя классификация YDR373W), отвечающий за поддержание механизма циркадианных ритмов в дрожжах и некоторых других грибах, обнаружен в количестве 7160 штук, приходящихся на одну клетку. Это означает, что можно пренебречь стохастическими свойствами системы, так ярко проявляющими себя на микроскопическом уровне описания, поскольку роль флук-туаций здесь невелика. С другой стороны, учитывая, что размер одной молекулы белка составляет порядка 0.001-0.01 мкм, можно ввести пространственно-распределенную модель описания динамики протеинового облака. Известно, что коэффициент диффузии разных белковых молекул в форме мономеров при невысокой их концентрации в водном растворе оценочно равен 10-7 см2/с [9]. Это значение несколько ниже стандартного значения для диффузии веществ в воде, так как мо-

номеры белка представляют собой крупные тяжелые молекулы, состоящие из сотен и тысяч аминокислотных остатков. Однако диффузия белка в цитоплазме живой клетки ограничена строением внутриклеточного пространства. Как сейчас уже стало понятно, ци-тоскелет клетки, состоящий из актиновых фи-ломентов, представляет собой сложную вязко-упругую среду, допускающую при определенных условиях ожижение [10, 11]. Значение коэффициента диффузии белковых клеток в таком внутриклеточном каркасе в данный момент неизвестно, но оно точно меньше вышеприведенного значения для воды. Вероятно, хорошей оценкой снизу для этой величины является значение для латеральной диффузии белка в клеточных мембранах - 10-12 см2/с [12].

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

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

2. Модель циркадианных колебаний

Рис. 2. Схема взаимодействий белков в модели циркадианных колебаний. Здесь: f, V -гены; F, Ж - соответствующие белки

На всех трех уровнях описания будем использовать динамическую модель циркадианных колебаний, предложенную авторами ранее в [13]. Несмотря на то, что модель предлагалась для описания суточных ритмов организма №и1герога crassa, она имеет достаточно общий характер. К примеру, в работе [14] модель подобного класса легко обобщается на случай циркадианных колебаний у мушки Drosophila. Главным элементом механизма колебаний является эффект запаздывания реакций синтеза белков в процессах транскрипции и трансляции генов (рис. 2). Эти процессы не просто медленные, но еще и состоят из многоэтапных биохимических реакций, в ходе которых последовательно образуются сложные органические соединения. Таким образом, эти процессы растянуты во времени, а значит идут с некоторым характерным временем запаздывания. Подробнее вопросы, связанные с функционированием механизма запазды-

вающих колебаний, обсуждаются в работах [7, 13].

На рис. 2 представлена схема взаимодействия двух белков, формально обозначенных как F и W, кодируемых двумя генами f и w. У ряда организмов были выделены такие пары принципиальных генов, отвечающих за работу механизма циркадианных колебаний. Например, в случае нейроспоры это гены frq и wcc, у дрозофилы - per и dclock [14]. Это не значит, что механизм колебаний поддерживается работой исключительно этой пары генов, так как циркадианные ритмы обладают широким набором свойств (автономность, способность компенсировать температурные изменения среды, а также изменять фазы колебаний под действием внешнего освещения или температуры и т.д.), задаваемых десятками генов. Однако, указанные пары являются принципиальными для поддержания циркадианных ритмов.

Как видно из рис. 2, обобщенная модель включает в себя как положительную, так и отрицательную петлю обратной связи и является симметричной по отношению к обеим переменным. Кроме того, она учитывает процессы димеризации и деградации белков. Полный список реакций, происходящих при транскрипции генов, приведен в табл. 1. Концентрации белков-мономеров в таблице обозначены F и W, димеров - F2 и W2. Динамика оператор-сайта промотора описывается бинарной функцией B X {B0, B1}, которая принимает значение B0 в случае его открытия и Bl в случае закрытия. Наиболее важными для поддержания колебаний являются реакции синтеза белков, протекающие с характерным временем запаздывания Ш. Именно этот механизм определяет динамику системы и её устойчивость к внешним и внутренним возмущениям.

Таблица 1. Список реакций транскрипции генов frq и м>с. Здесь и далее ^ k^2, уг, уЖ,

Процесс димеризации F G F I kF I? f2 , W g W I F?? W2

Процесс дедимеризации F2 I kF'ii F g F, W2 I kW? W g W

Динамика оператор-сайтов DF g w2 I kF ?? DF , d0w G f21 kw ? dww, DF I k2 ?? DF G w2 , DF I f ?? DW g f2

Процесс синтеза белков DF(t) I f* ?? DF g FtGШ , D? (t) I \w ? D? G W'Ш

Процесс распада белков F I I" ? й , W I Г? ? й

Процесс образования гетеродимера F G W I Ik ? й

Таблица 2. Параметры модели

100 k kF kw kF K 2F K к w V °W

6 h 30 nM-1h-1 8 nM/h 4 nM/h 5 nM-1 5 nM-1 5 nM-1 5 nM-1 0.3 h-1 0.4 h-1

На основании цепочки связанных биохимических реакций, представленных в табл. 1, может быть получена динамическая модель циркадианных ритмов (1-2) в детерминистском описании [13]:

1

ff)dF т %i1 Я

(1G 4K F) ^ T kp +1

— ¡Я

1G KKW2(t Я 1) | (1) Я °fF Я kFW,

dW ^ (1 G 4KWW)— T kW +1Я dt |

1

1G KKF2(tii)Iя

Я 0WW Я kFW,

(2)

F a F F

где K. s k. / k я

K

W 4

1 W I i W D

k. / k,. . В численных

Н и ^г * ^г ' г

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

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

ческих реакций (табл.1), либо упрощенной детерминистской моделью (1-2).

3. Микроскопический уровень

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

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

Удобным инструментом исследования стохастических генетических моделей является метод Гиллеспи [15], который представ-

time (hours) time (hours)

Рис. 3. (а) Временная зависимость полного числа мономеров белка F (сплошная линия) и W (штриховая линия), полученная при численном интегрировании системы (1-2). (b) Временная эволюция соответствующей стохастической системы, полученная с помощью алгоритма Гиллеспи, модифицированного для немарковских процессов. Во врезках приведены фурье-спектры сигналов. Значения параметров системы в обоих случаях одинаковы (табл. 2)

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

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

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

На рис. 3 а представлены развертки по времени для полного числа мономеров обоего сорта, полученные при численном расчете де -терминистской динамической системы (1-2) в отсутствие базальных транскрипционных факторов. На рис. 3Ь приведены результаты стохастических расчетов для тех же значений параметров, что и в детерминистском случае. Хорошо видно, что, несмотря на небольшое

Можно заметить, что полной симметрии между колебаниями концентраций не наблюдается, так как несимметричен набор значений параметров, приведенный в табл. 2. Стохастич-ность системы проявляет себя во флуктуаци-ях, которые достигают 20-40% от амплитуд колебаний в детерминистском случае. Во врезках к рисункам приведены фурье-спектры сигналов. Можно заметить, что стохастический сигнал имеет устойчивый максимум, примерно соответствующий периоду колебаний 23 часа. Подобное поведение системы еще раз подчеркивает устойчивость базисного механизма генерирования циркадианных колебаний, основанного на запаздывании в процессах транскрипции и трансляции.

4. Мезоскопический уровень

Моделирование синхронизации биоритмов в масштабах одной или нескольких клеток (10-100 мкм, рис. 1) предполагает построение пространственно-распределенной модели. Для этого модель циркадианных колебаний (1-2) была дополнена диффузионным слагаемым, где D - коэффициент диффузии молекул белка внутри клетки. Уравнения (3-4) сформулированы для переменных F(X, х, у) и Ж(X, х, у), которые теперь являются функциями времени и двух пространственных координат х, у. В данной статье мы ограничиваемся рассмотрением пространственно-временной динамики системы в двух пространственных измерениях.

Численное исследование системы (3-4) осуществлялось с помощью метода конечных

Ж

т

1

в х (1 в 4к;ж)

й1

к , I, А в2р , 8 2рI

л ,а й ■) н н к¥ж х с я £4Л Ьр- «■

в К1 К 2 ¥ (X й ■) х 18 х 8 у X

(4)

количество вовлеченных в динамику молекул разностей, который является одним из наибо-(в любой момент времени полное число моно- лее распространенных численных методов ре-

р

1

8 X (1 Л 4К

¥¥) £

кр й

1 с КЖК рж 2(х Й ■)

й орр й кРЖ{ Л я

А82р, 8 2р I

;8 х2

8 У 2 X

8 У X

(3)

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

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

----relative rate of change of function density of base slices--time steps

Рис. 4. Зависимость количества опорных слоев (обозначены как BS) от гладкости функции

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

1 t

h2

2 max F ,W

(5)

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

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

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

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

уrUf -'с,-а

./ г ^ -i

С

j ч \ ■- гч. ^ г . ,

таfrjil

Рис. 5. Характерный вид волновых структур, образуемых белком F во внутриклеточном пространстве из случайных начальных условий в условиях полной темноты. D=0.01

Fj(x,0) t F2(x,0) t

D(4^(x,Hh) H F1(x,H2h)) G D0(4F2(x,h) H F2(x,2h))

(6)

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

Внутри области неустойчивости возникают многочисленные центры возбуждения вторичных спиральных волн, нелинейное взаимодействие которых приводит к хаотическому поведению среды за гребнем первичной волны (рис. 5). Расчёт взаимно-корреляционной функции между парами точек среды, принадлежащих областям до и после прохождения волны, показывает, что корреляции между точками в волновой зоне практически отсутствуют [13]. Таким образом, обнаруженные волны можно трактовать как волны рас-синхронизации, которые сбрасывают все локальные установки фаз часов в среде. Тем не менее поведение системы нельзя назвать пассивным - она, скорее, находится в "ждущем" режиме по отношению к внешней среде. Хотя поле концентрации белка Р (с белком Ж происходит то же самое) хаотически осциллирует как во времени, так и по пространству, но при определенной команде извне система готова быстро синхронизировать свои колебания в соответствии с внешним сигналом.

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

3(D G Do)

Моделирование данного процесса выполнено следующим образом: в расчетной области (рис. 6) выделяются определенные участки, которые мы будем условно называть "клетками", в которых протекают реакции (3-4), задающие циркадианные ритмы. Так как механика самих клеток в данный момент нас не интересует, то модельные клетки имеют форму квадрата и фиксированный размер. Граница клеток не однородна: большая её часть непроницаема для реагента (отмечена на рис. 6 черной линией), а оставшаяся часть имеет диффузионную проницаемость, вычисляемую по конечно-разностной формуле (6) [18, 19], где h - шаг по сетке, D0 - коэффициент диффузии белка во внеклеточном пространстве.

Таким образом, каждый участок, в котором протекает реакция, является моделью живой клетки с собственным циркадианным ритмом внутри. Каждая клетка имеет четыре мембраны фиксированной ширины, через которые она может обмениваться с межклеточным пространством. Для простоты мы не стали вводить в модель отдельный вид белка, ответственный за транспортировку циркадиан-ного сигнала наружу. Транспортную функцию в предлагаемой модели выполняет всё тот же белок F, важный для поддержания циркадианных ритмов. Его партнер, белок W , может перемещаться только внутри пространства клеток - для него клеточные мембраны закрыты.

I II I

cell #1

I__I

extracel

Г

cell #2

I__I

ular space

Г и

cell #3

cell #4

I__II__I

Рис. 6. Схема расчетной области

Рис. 7. Динамика концентрации белка Р в полной темноте. Только в одной клетке (#3) происходит реакция синтезирования белка. В остальных клетках (#1, #2 и #4) белок только деградирует. Последовательно представлены моменты времени (=200, 5000, 10000, 20000. Коэффициент диффузии в клетках равен Я=0.01, в межклеточном пространстве Я=0.01

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

3/7 A3 2 г 3 2/7 1

8 F t D0 fe^F G ^ j.

8 t

:8 x2

2

8 У X

(7)

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

Перейдем к результатам моделирования системы с клетками. На рис. 7 представлены четыре кадра динамики изменения концентрации белка Р в полной темноте и в отсутствие факторов базальной транскрипции, когда только одна из клеток (#3) генерирует цирка-дианный ритм. В остальных трёх клетках механизм транскрипции выключен: кР Т 0. Как было установлено выше (рис. 5), система в этом случае находится в состоянии фазовой турбулентности, которое в отсутствие внешнего сигнала может сохраняться сколь угодно долго. Несмотря на то, что существует поток белка Р через мембраны клетки #3 в межкле-

ris>>- T^kîv т-ësï

точное пространство, клетка продолжает функционировать в прежнем режиме хаотических пространственных волн. Убыль белка компенсируется его синтезом внутри клетки. Попадая в межклеточное пространство, выполняющий транспортные функции белок F либо диффундирует там, либо проникает в клетки #1, #2 и #4. Механизм циркадианных колебаний в этих клетках парализован, однако процессы линейной и нелинейной деградации белка действуют. Поэтому концентрация F внутри выключенных клеток быстро падает до нуля (рис. 7). Представленный на рисунке сценарий не зависит от размеров клеток и ширины межклеточного зазора между ними. Синхронизация колебаний в клетке #3 на рис. 7 не наступает из-за отсутствия внешнего сигнала.

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

ES г ч

f.

m

ш ш

Рис. 8. Динамика концентрации белка Р в полной темноте. Все четыре клетки функционируют, вырабатывая белок Р. Последовательно представлены моменты времени (=200, 5000, 10000, 20000. Коэффициент диффузии в клетках равен Я=0.01, в межклеточном пространстве Я=0.01

Рис. 9. Динамика концентрации белка F в системе из 16 плотно примыкающих друг к другу клеток с разными фазами циркадианного ритма в полной темноте. Последовательно представлены моменты времени t=0, 2000, 4000, 6000, 8000, 10000. Коэффициент диффузии в клетках равен Б=0.01

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

Численные эксперименты показали, что важную роль в синхронизации играет ширина межклеточного зазора. Чем шире зазор, тем больше времени требуется для достижения состояния синхронных колебаний. Межклеточное диффузионное пространство играет здесь точно такую же роль, какую выполняют пружинки связи в системе связанных механических осцилляторов [5]. Чем сильнее связь между осцилляторами, тем легче достигается синхронизация. На рис. 9 представлены кадры эволюции 16 клеток, которые настолько тесно расположены друг к другу, что межклеточное пространство здесь практически отсутствует. На самом деле такая конфигурация системы более реалистична, так как клетки тканей органов обычно образуют плотную сплошную среду без зазоров. На рис. 9 клетки продолжают обмениваться друг с другом сигналами посредством того же самого диффузионного механизма, что и раньше. В начальный момент эволюции в каждой клетке задавалось постоянное для этой клетки значение концентрации обоих взаимодействующих белков, однако сами значения в разных клетках задавались случайным образом. Так как уровень концентрации белков задаёт фазу циркадианного ритма (максимальный уровень F определяет локально «полдень»), то система интегрировалась со случайным значением фазы колебаний в каждой клетке. При этих начальных условиях спиральные волновые структуры, приведенные на рис. 5, в системе не возни-

кают. Для этого нужны резкие локальные флуктуации концентрации белка.

В случае, представленном на рис. 9, реализуется решение в виде стоячей волны. Хорошо видно, что сравнительно быстро фаза всей группы клеток выравнивается, - происходит фазовая синхронизация. Начиная примерно с момента времени t Т 4000, система осциллирует с общей фазой. Отметим, что синхронизация здесь наступает гораздо раньше, чем в случае, показанном на рис. 8.

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

5. Макроскопический уровень

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

Модель, предложенная в статье [20], включает расчет динамики отдельных клеток, представленных в виде многоугольников с разным числом вершин (система откалиброва-на так, что наиболее вероятной формой клет-

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

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

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

dР4

Т

1

\кР Н

кР

dt (1 в 4Кр Fl) у 1 й К^ КFWl 2 ( Н

- возможность перемещения клеток в общей массе эпителия посредством механизма ин-теркаляции;

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

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

- учет эффекта поляризации клеток, которая происходит спонтанно или под воздействием внешних условий.

Рис. 10. Элементы хемомеханической модели эпителия

-Н °рР, Н kFlWl * й~ А (р Н Fl). (8)

V Ж ]

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

Модель эпителия мы дополнили расчетом циркадианных ритмов в каждой клетке сообщества согласно уравнениям (8-9), где ~ Т 0.1 - коэффициент переноса белка Р через клеточную мембрану.

dWi

1

dt (1 й 4KÍWW1)

кж Н

к

W

ж

1 й Кр К 2[р;

(t Н

0 Ж Н кРЖ *

Ж ГГ I и П,! ^Г I

Ж

I Т 1, N^). (9)

Рис. 11. Эффект кластеризации циркадианных ритмов в эпителиальной ткани: колебания концентрации белка Р в условиях полной темноты в ткани эпителия, состоящего из более 1600 клеток. В качестве начальных условий использовано случайное распределение фазы колебаний в клетках. Последовательно представлены моменты времени (=0, 25, 50, 75, 100, 125

Слагаемое, ответственное за перенос белка из клетки в клетку в уравнении (8) записано в виде простого разностного соотношения, в котором Lij имеет смысл длины границы между /-ой и у-ой клетками (рис. 10). Суммирование в формуле (8) ведется только по клеткам, которые являются соседними для /ой клетки. Если уровень концентрации белка Р в рассматриваемой клетке выше, чем у её окружения, то формируется поток молекул белка наружу. В обратном случае наблюдается приток белка извне. Транспортную функцию, как и прежде, выполняет только белок Р, а для Ж границы клеток остаются непроницаемыми. В результате деления какой-то клетки на две части считается, что новые клетки наследуют фазу циркадианного ритма материнской клетки. Так как в данной статье нас не интересует влияние циркадианных ритмов на хемомеханические свойства среды, обратная связь уравнений (8-9) с общим функционированием клеточной модели отсутствует.

Несмотря на то, что уравнения (8-9) в отличие от (3-4) представляют собой систему обыкновенных дифференциальных уравнений, общее их количество достаточно велико, так как N (количество клеток) может достигать нескольких тысяч. Кроме того, при рас-

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

Перейдём к изложению результатов моделирования. На рис. 11 представлены несколько кадров эволюции эпителия, состоящего в начальный момент из 1600 клеток. В качестве начального условия в клетках случайным образом задавалась фаза циркадиан-ного ритма, а внешние воздействия на систему исключены. Такой подход позволяет определить возможные собственные формы коллективного поведения клеток по синхронизации колебаний. Результаты расчетов показывают, что в пространстве большого количества клеток, полной синхронизации ритмов в смысле абсолютного выравнивания фазы колебаний по всем клеткам, как это происходило в случае небольшого сообщества (рис. 9), не получается. Здесь проявляется макроскопический эффект кластеризации колебаний -клетки формируют два примерно равных сообщества, которые коллективно осциллируют в противофазе (рис. 11). Между двумя группами есть небольшая прослойка клеток, в которых реализуются колебания с промежуточными значениями фазы.

Рис. 12. Эффект бегущей волны циркадианных ритмов в эпителиальной ткани: колебания концентрации белка Р в условиях полной темноты в ткани эпителия, состоящего из более 1600 клеток. В качестве начальных условий использовано возмущение концентрации белка в одной клетке в центре сообщества. Последовательно представлены моменты времени (=0, 20, 40, 60, 80, 100

Кластеризация в системах с большим6. количеством элементов, обменивающихся химическими сигналами, с некоторых пор привлекает к себе внимание исследователей. Например, в недавней работе [21] подробно изучена растущая группа взаимодействующих друг с другом синтетических генетических осцилляторов. Обнаружено, что с течением времени происходит кластеризация ткани на два типа осциллирующих клеток. Авторы связывают это явление с двумя возможными устойчивыми состояниями равновесия у системы. Отмечается, что кластеризация, по-видимому, является важнейшей характерной особенностью больших сообществ и может служить причиной дальнейшей дифференциации клеток в органах.

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

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

Обсуждение

Вопрос о сшивании результатов моделирования, полученных на разных уровнях описания системы, является важным для понимания того, как интерпретировать результаты. Мы рассмотрели динамику циркадианных колебаний в масштабах гена (стохастическая модель), клетки (детерминистская пространственно-распределённая модель с диффузией) и организма (дискретная хемомеха-ническая модель клеток эпителия). На всех уровнях описания в модели был заложен генный механизм генерирования ритмов на основе бинарного взаимодействия двух белков (рис. 2). Соответственно, базовые колебания с периодом примерно 23 часа проявляют себя во всех моделях. Было показано, что механизм колебаний устойчив к внутреннему шуму, связанному с кинетикой реакций (рис. 3). В отсутствие сильных возмущений в масштабах клетки колебания быстро синхронизируют фазу колебаний по всему внутриклеточному пространству. В случае взаимодействия нескольких клеток обмен диффузионными потоками приводит к синхронизации ритмов уже у всей группы клеток (рис. 9).

Если внутри клетки есть значительные локальные отклонения концентрации, то система может перейти к новому режиму осцил-ляций, представляющему собой хаотическое поле бегущих волн (рис. 5). Характерный размер пространственных структур, возникающих в этом случае, примерно на порядок меньше размеров самой клетки. Сценарий синхронизации ритмов таких возбужденных клеток, взаимодействующих друг с другом, несколько отличается от прямолинейного выравнивания фазы. За счет транспортного белка в межклеточном пространстве устанавливается общее для всех клеток среднее поле, которое приводит к регуляризации бегущих волн и пространственно-временной синхронизации их распространения (рис. 8). Важность взаимодействия между клетками легко установить, если отключить процессы синтеза циркадиан-ных белков во всех клетках, за исключением одной (рис. 7). В этом случае состояние хаотически распространяющихся бегущих волн может поддерживаться сколь угодно долго.

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

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

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

В литературе существует ряд экспериментальных свидетельств, что процессы транскрипции/трансляции белков происходят с запаздыванием. Например, в работе [22] для организма №и1герога crassa было установлено, что такая запаздывающая обратная связь формируется у пары белков FRQ и WCC, синтезируемых генами frq и м>сс. Что особенно важно, именно эта пара белков ответственна за поддержание циркадианных ритмов у нейроспо-ры. В данный момент трудно сказать, является ли запаздывание фундаментальным механизмом, присущим транскрипции любых генов. Однако то, что это важное свойство хотя бы некоторых из них, установлено достоверно.

Современные технологии позволяют экспериментаторам вместо изучения уже существующих организмов экспериментировать с искусственно созданными цепочками ДНК. Именно таким путем пошли авторы работы [23]. Ими была создана синтетическая генная цепочка, которая включала в себя осциллирующую часть, основанную на паре генов Ьих1 и LuxR. Колебания уровня соответствующих белков достигались за счет запаздывания

обратной связи. Кроме того, цепочка включала в себя ген, синтезирующий транспортный белок AHL (N-ацил-гомосерин лактон! Для искусственной визуализации групповой динамики в каждую цепочку был добавлен ген yemGFP, вызывающий свечение под действием транспортного белка AHL. Синтетическая цепочка была интегрирована в каждую бактерию кишечной палочки, а вся колония помещена в специальную среду с питательным раствором. Таким образом, внедренные осцилляторы имитировали в каждом организме работу его биологических часов, а белок AHL осуществлял глобальную координацию поведения всей колонии. Численность колонии достигала несколько тысяч особей. В результате наблюдения за колонией бактерий было обнаружено спонтанное возникновение бегущих концентрических волн свечения. В сопроводительных файлах на сайте журнала Nature можно ознакомиться с различными видеозаписями экспериментов. Отметим, что концентрические волны, распространяющиеся в хемомеханической модели эпителия (рис. 12), в деталях воспроизводят пространственно-временную динамику в колонии сообщающихся друг с другом бактерий.

Заключение

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

Авторы выражают свою признательность профессору Л.М.Писмену (L.M. Pismen) и доктору М. Салму (M. Salm) за разрешение

использования их модели растущего эпителия

в данном исследовании.

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

1. Ашофф Ю. Биологические ритмы. М.: Мир, 1984. С. 412.

2. Pittendrigh C.S. Temporal organization: reflections of a Darwinian clock-watcher // Annu. Rev. Physiol. 1993. V.55. P. 16-54.

3. Моисеева Н.И., Сысуев В.М. Временная среда и биологические ритмы. Л.: Наука, 1981. С. 127.

4. Степанова С.И. Биоритмологические аспекты проблем адаптации. М.: Наука, 1986. С. 244.

5. Pikovsky A., Rosenblum M., Kurths J. Synchronization - а universal concept in nonlinear sciences. Cambridge University Press, 2001. P. 432.

6. Hasty J., Collins J.J. Translating the noise // Nature Genetics. 2002. V. 31. P. 13-14.

7. Брацун Д.А. Эффект возбуждения подкри-тических колебаний в стохастических системах с запаздыванием. Ч. I. Регуляция экспрессии генов // Компьютерные исследования и моделирование. 2011. Т. 3, № 4. С.431-438.

8. Ghaemmaghami S., Won-Ki H., Bower K., Howson R. W., Belle A., Dephoure N., O'Shea E.K., Weissman J.S. Global analysis of protein expression in yeast // Nature. 2003. V. 425. P.737-741.

9. Якубке Х.-Д., Ешкайт Х. Аминокислоты, пептиды, белки. М.: Мир, 1985. С. 457.

10. Koenderink et al. An active biopolymer network controlled by molecular motors // PNAS. 2009. V. 106, № 36. P. 15192-15197.

11.Morozov K.I., Pismen L.M. Cytoskeleton fluidization versus resolidification: Prestress effect // Phys. Rev. E. 2011. V. 83. P. 051920-051928.

12. Северина Е.С. Биохимия: учебник для вузов. М.: ГЭОТАР-Медиа, 2003. С. 784.

13. Брацун Д.А., Захаров А.П. Моделирование пространственно-временной динамики цир-кадианных ритмов Neurospora crassa //

Компьютерные исследования и моделирование. 2011. Т. 3, № 2. С. 191-213.

14. 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. V. 7, № 4. P. 337-354.

15. Gillespie D.T. Exact stochastic simulation of coupled chemical reactions // J. Phys. Chem. 1977. V. 81. P. 2340-2361.

16. Bratsun D., Volfson D., Hasty J., Tsimring L. Delay-induced stochastic oscillations in gene regulation // PNAS. 2005. V. 102, № 41. P. 14593-14598.

17. Брацун Д.А., Захаров А.П. К вопросу о численном расчете пространственно-распределенных динамических систем с запаздыванием по времени // Вестник Пермского университета. Сер.: Математика. Механика. Информатика. 2012. Вып. 4(12). С. 3241.

18. Bratsun D.A., De Wit A. Buoyancy-driven pattern formation in reactive immiscible two-layer systems. // Chem. Eng. Sci. 2011. V. 66. № 22. P.5723-5734.

19. Bratsun D.A., De Wit A. Control of chemoconvective structures in a slab reactor // Techn. Phys. 2008. V. 53. P. 146-153.

20. Salm М., Pismen L.M. Chemical and mechanical signaling in epithelial spreading // Phys. Biol. 2012. V. 9. P. 026009-026023.

21.Koseska A., Ullner E., Volkov E., Kurths J., Garcia-Ojalvo J. Cooperative differentiation through clustering in multicellular populations // J. Theor. Biol. 2010. V. 263. P. 189-202.

22.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.

23. Danino T., Mondragón-Palomino O., Tsimring L., Hasty J. A synchronized quorum of genetic clocks // Nature. 2010. № 463. P. 326-330.

Modeling biorhythms in the scale of a gene, a cell and a whole organism

A. P. Zakharov, D. A. Bratsun

Perm State Humanitarian Pedagogical University, Russia, 614990, Perm, Sibirskaya st., 24

az1211@mail.ru; +79222465522

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.

Key words: biorhythms; synchronization; systems with time-delay.