О поиске универсальных параметров для моделей морфогенеза Текст научной статьи по специальности «Математика»

Владикавказский математический журнал 2013, Том 15, Выпуск 3, С. 58-66

УДК 517.958:57+519.710.22

О ПОИСКЕ УНИВЕРСАЛЬНЫХ ПАРАМЕТРОВ ДЛЯ МОДЕЛЕЙ МОРФОГЕНЕЗА

М. Н. Назаров

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

Ключевые слова: моделирование морфогенеза, межклеточный сигнальный обмен, согласование моделей разных классов.

Введение

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

© 2013 Назаров М. Н.

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

дискретного пространства и дискретного времени, а клетки в рамках имитационных моделей представляются как объекты фиксированного объема и формы. Все ключевые составляющие морфогенеза, такие как деление клеток, настройки типов и т. д. в подобных моделях воспроизводятся с помощью определенного набора правил. Классическим примером имитационной модели является модель клеточных автоматов (см. [8]).

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

1) наличие биологического смысла для каждого из параметров;

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

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

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

1. Выбор универсальных параметров

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

• T = {ti, ..., тп} — типы клеток;

• is (ti), ..., Is (тп) — интенсивности синтеза сигналов;

• Ре(т\ ,ti ),..., Ре(тп,тп) — вероятности передачи сигнального воздействия клетке при столкновении с сигналом ( 0 ^ Ре ^ 1 );

• Pr(Ti,Ti),..., Рг(тп,тп) — приоритеты сигналов ( Pr ^ 0 );

• D(t1,t1 ),..., В(тп, тп) — новые типы для дочерних клеток;

• K(т1,т1),..., K(тп, тп) — элементарные приращения скорости наступления деления или программируемой смерти;

• V(т1),..., V(тп) — объемы клеток ( Vmax = max V(т) ).

т

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

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

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

№ (т ),..., Бк (т )} .

Для определенности положим, что в единицу времени клетка выделяет в среднем: /j (т) сигналов Sj(т), где ] = 1 ,к. Тогда в качестве /^(т) выберем /^(т) = юш^^(г). В этом случае можно ввести усредненный сигнал Б(т) как набор реальных сигналов:

Вероятности Ре(т, т*) передачи сигнального воздействия со стороны сигналов Б(т*) клетке т напрямую зависят от способа передачи сигнального раздражения внутрь клетки. Для растворимых в жирах сигнальных молекул Sj(т*) вероятность их передачи через мембрану клетки определяется растворимостью. Если же сигналы Sj• (т*) для передачи воздействия связываются с рецепторами, то вероятность передачи зависит от количества данных рецепторов на мембране клеток типа т. Зная вероятности Ре(т, Sj• (т*)) для отдельных химических сигналов, можно приблизительно оценить вероятность Ре(т, т*) для усредненного по формуле (1) идеализированного сигнала Б(т*) как:

1 к

Ре(т,п) = ^ (2)

j(Т*) j=1

Приоритеты сигналов Рг(т,т*) тем больше, чем сильнее клетка реагирует на сигнальное раздражение. Для отдельного химического сигнала Sj (т*) приоритет можно по определению считать равным суммарной биологической активности вторичных посредников5, которые были задействованы при передаче раздражения внутрь клетки типа т. Зная приоритеты Рг (т, Sj(т*)) для отдельных химических сигналов, можно оценить приоритет Рг(т,т*) с помощью усреднения, полностью аналогичного формуле (2):

1к

Рг(т,п) = ^1,(п)-Рг(т,Б,(п)). (3)

lj (Т*) j=l

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

3 В простейшем случае состояние гена можно определить как количество актов копирования гена в единицу времени.

4с 8т

4Форма покоя — это форма клетки в отсутствии сокращения.

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

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

2. Описание базовой количественной модели

Переменными количественной модели будут:

• N(т, t) — численности клеток как функции от типов т;

• Vd(T,t) — интенсивности деления/программируемой смерти. Сигнальное воздействие одних клеток на другие можно посчитать как:

Inf(г f, N) > 0;

Sig(т <- т, n)= i SumInf(r,W)> V \0, else.

В этом выражении SumInf(т, N) = f Inf(T ^ T,N) — это суммарное сигнальное воздействие.

SumInf (т, N)

■Sig(r ^ tuN) ■Sig(T ^ Tn,N)

Inf (т ^ т\, N) Inf (т ^ Tn,N)

Рис. 1. Распределение суммарного сигнального воздействия.

Для аннулирования эффекта от клеток, чья численность меньше единицы6 введем приведенные численности по формуле: Ж(т) = а (](т) — 1 + с) N(т). В качестве функции а можно использовать как функцию Хевисайда, так и сигмоиду а(х) = 1/(1 + е-х/с), где 0 < с ^ 1. При этом отдельное элементарное сигнальное воздействие:

Inf (т ^ T, N) =

'Рг(т,т) ■ Ре(т, т) ■ Рг(т, t)N(t) > 0;

0,

Рг(т, T)N(T) = 0.

Коэффициент Dispf (Ж, ¿) позволяет учесть распределение сигнального раздражения между всеми подверженными ему клетками:

Disp(N )= ^ ]У(т V).

Т г*:рт(г*,Т)> 0

Учитывая потери и прирост клеток получаем7 выражение для N(т, ¿): ё] (т V)

dt

■ v^t) ■ ]Ч(т, t) ■ Sig^ ^ T,N^Л(т*,t) ■ N(т*,t). (4)

t f: d(t,f)=t *

6По умолчанию такие клетки считаются полностью вымершими в рамках модели.

7За элементарный промежуток времени ровно • Й(т, ¿) • dí клеток типа т делятся или

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

Аналогично строится выражение для пересчета интенсивностей деления/программируемой смерти Vd(r, t):

d Vd(T,t) dt

Y^ K (r,f) ■ Sig(r ^ f, N)

f: K(r,f)>0

+ a (vd(T,t) - c) ■ Y, K (T,f) ■ Sig(r ^ f, N).

f: K(r,f)<0

(5)

3. Описание базовой имитационной модели

Имитационная модель вводит клетки и сигналы как самостоятельные объекты, локализованные в пространстве модели, а также использует дискретное время. Для этого требуются дополнительные параметры:

• а — ребро элементарного куба;

• Н, Ш — размеры пространства в элементарных кубах;

• — элементарный такт времени.

H = 3

W = 1

Рис. 2. Пример пространства.

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

Cell(t) : {1,..., H} х {1,..., L} х {1,..., W} —► C U {0},

где C — это множество всех клеток. Отображение Cell(t) ставит в однозначное соответствие тройке дискретных координат (z,x,y) : z G 1, H, x G 1, L, y G 1, W содержимое8 ячейки на момент времени t.

Для упрощения работы с дискретным пространством кубов введем специальные обозначения (см. сводную таблицу 1):

• d G Directions = {|, j , 0, 0} — направления;

• d G Directions — противоположное направление для d;

• d(z, x, y) — соседний куб с (z, x, y) по направлению d;

8Для пустой ячейки будет выполняться Cell(Z;X>y)(t) = 0.

Dirsig(z,x,y)(t) — допустимые направления для передачи сигнала:

Dusig(z,x,y)(t) = Id : d(z, х, у) = (z,x,y) Л ( ¿еТл ) Л (Се%й>й) (i) ф 0)1 ;

I \y€l,W/ )

Dir¿iv (z,x, y)(t) — допустимые направления для деления клетки:

Dirdiv(z,x,y)(t) = id : d(z,x,y) = (z,x,y) Л ( ¿efj ) Л (Cell^^^)(í) = I \yei,w/

Таблица 1

Специальные обозначения для дискретного пространства

d т 1 —}■ <— © ©

d 1 т <— —}■ © ©

d(z, х, у) (z + 1, X, у) (z - 1, X, у) (z, X + 1, у) {z,x- 1, у) (z, X, у + 1) (z,x,y - 1)

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

• Is(At,a)(r) = Is(т) ■ At ■ (a3/V(т));

• K(At,a)(T,T*) = K(т,т*) ■ At ■ (V(т)/a3).

Переменными имитационной модели будут идеализированные клетки

Cell(z,x,y) (t) = (т,Уа,фа,Фз, Sin, Sig, Sout, Plan),

где t — тип клетки; Vd — интенсивность наступления деления/смерти; фd — фаза деления/смерти; ф8 — фаза синтеза сигналов; Sin — семейство сигналов из ближайшего окружения клетки; Sig — семейство сигналов, которые были приняты клеткой; Sout — семейство сигналов, ожидающих ретрансляции; Plan = (Tnew, Prnew) — план9 деления.

Динамика имитационной модели сводится к последовательной реализации для каждой Cell(z x y) (t) следующего набора правил:

Этап 1: Прием сигналов. Для каждого сигнала из входного семейства S Е Sin(t) разыгрывается случайное событие приема сигнала. С вероятностью Ре(т, тs) сигнал S = (ts) будет перемещен в семейство Sig(t), и с вероятностью 1 — Pe(T,Ts) сигнал считается прошедшим мимо клетки и попадает в семейство Sout(t).

Этап 2: Корректировка Vd. Для начала рассчитывается суммарный приоритет всех сигналов, поступивших внутрь клетки:

SumInf(t)= Y^ Pr (т, тs).

S 6 Sig(t)

Приращение скорости Vd можно посчитать по формуле:

Рг(т, TS )

Avd(t) =

S6 Sig(t)

SumInf (t)

K (At,a)(T,Ts).

9План деления задает какой тип (тпеш) получат дочерние клетки после деления и насколько приоритетным (Ртпеш) будет такое деление.

Корректировка интенсивности Vd будет осуществляться с учетом того, что Vd не должно становиться отрицательным:

ivd(t) + Avd(t), Vd(t) + Avd(t) ^ 0; Vd(t + 1)^0, else.

Этап 3: Корректировка плана деления. Для всех возможных новых типов Td (в том числе и для 0, моделирующего программируемую смерть) вычисляется суммарный приоритет

Sumpr^t) = ^ Рт(т,тв). sesig(i),

D(t,ts )=rd

Осуществляется поиск типа клеток тЦ, для которого суммарный приоритет максимальный:

Maxpr(i) = max (Sump,.^, i)), fd(t) : Sumpr(7^, t) = Maxpr(t).

td

В случае если максимальный приоритет Maxpr (t) превосходит Prnew(t), то происходит обновление плана деления:

Tnew((i + 1)) - Td(t), если Prnew(t).

Pr new ((t + 1)) = Maxpr (t), РГК' neWK '

Этап 4: Синтез и передача сигналов. Добавляем в семейство Бои^) сигналы Б = (т), количество которых задается параметром . Пусть для определенности (т) = т + к/п, где т,к,п € N. Тогда количество сигналов Б = (т) может быть вычислено по формуле т(£) = т + о{ф8($) — 1 + к/п), где <т — это функция Хевисай-да. Далее для каждого сигнала из семейства Бои^) разыгрывается случайное событие передачи: с равными вероятностями сигнал передается любой из клеток ближайшего10 окружения и помещается в семейство Бт(£ + 1) этой клетки. После завершения передачи пересчитываем фазу синтеза сигналов по формуле:

Этап 5: Деление/программируемая смерть. В зависимости от значения V^ рассчитывается количество элементарных актов деления. В случае, если тnew =0 и запланировано было деление, то клетка создает необходимое количество новых клеток в пустых соседних ячейках и меняет свой тип на тпе^ Пусть для определенности Vd(т) = т + к/п, где т,к,п € N. Тогда количество актов деления может быть вычислено по формуле т(£) = т + а(ф(1^) — 1 + к/п), где а — это функция Хевисайда. Для новых клеток используются параметры (тп^,Vd, 0, 0, 0, (тп^, 0)), где Vd берется от материнской клетки. После завершения этапа осуществляем пересчет фазы деления как:

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

10Ближайшим окружением клетки Се11(2>Х;У) будем считать все клетки из множества Б1г818(,г, х, у)(^), а также саму Се11(2>Х;У).

Заключение

В первую очередь нужно отметить, что обе базовых модели создавались в демонстрационных целях, и как следствие, их применимость на практике для описания динамики реальных организмов будет существенно ограничена. Так, например, базовая имитационная модель может дать приемлемо точное описание динамики формы лишь для организмов, чьи клетки не участвуют в миграции в пределах ткани. Проанализируем базовую количественную модель на предмет существования и единственности решения. Непрерывность по переменной t выполняется всегда для системы уравнений, так как прямой зависимости от времени в уравнениях (4) и (5) нет. Условие Липщица по переменным Vd и N будет11 выполняться, если в качестве a используется сигмоида a(x) = 1/(1 + e-x/c).

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

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

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

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

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

Рассмотренный в рамках данной работы набор универсальных параметров может быть обобщен как минимум двумя путями. Первое и наиболее очевидное обобщение набора универсальных параметров можно получить, если отказаться от усредненных сигналов (1) и ввести реальные химические сигналы в формализм моделей. При этом отпадает необходимость в использовании усредненных вероятностей (2) и приоритетов (3), поскольку вместо них можно будет работать напрямую с вероятностями и приоритетами настоящих химических сигналов. Второе обобщение сводится к включению альтернативных путей передачи информации между клетками в дополнение к химическим сигналам. Примером такой альтернативы является включение механических напряжений в модель (см. [9, 10]).

Литература

1. Лазарева Г. Г., Миронова В. В., Омельянчук Н. А., Шваб И. В., Вшивков В. А., Горпинченко Д. Н., Николаев С. В., Колчанов Н. А. Математическое моделирование морфогенеза растений // Сиб. журн. вычисл. мат.—2008.—Т. 11, вып. 2.—С. 151-166.

2. Diaz J., Alvarez-Buyalla E. R. A model of the ethylene signaling pathway and its gene response in Arabidopsis thaliana: pathway cross-talk and noise-filtering properties // Chaos.—2006.—Vol. 16, № 2.— P. 100-115.

11 Это не сложно показать, если учесть, что коэффициент DispT(N) либо равен нулю, либо строго больше единицы.

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

3. Mendoza L., Thieffry D., Alvarez-Buylla E. R. Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis // Chaos.-1999.-Vol. 15, № 7.—P. 593-606.

4. Lexa M., Genkov T., Malbeck J., Machakova I., Brzobohaty B. Dynamics of Endogenous Cytokinin Pools in Tobacco Seedlings: a Modelling Approach // Annals of Botany.—2003.—Vol. 91, № 5.—P. 585597.

5. Dumais J., Shaw S. L., Steele C. R., Long S. R., Ray P. M. An anisotropic-viscoplastic model of plant cell morphogenesis by tip growth // Int. J. Dev. Biol.-2006.-Vol. 50, № 2-3.—P. 209-222.

6. Geitmann A., Ortega J. K. E. Mechanics and modeling of plant cell growth // Trends in Plant Science.— 2009.—Vol. 14, № 9.—P. 467-478.

7. Mjolsness E., Sharp D. H. and Reinitz J. A connectionist model of development // J. Theor. Biol.— 1991.—Vol. 152, № 4.—P. 429-453.

8. Urdy S. Principles of morphogenesis: the contribution of cellular automata models // Acta Zoologica.— 2009.—Vol. 90, № 2.—P. 205-208.

9. Belousov L. V. The Dynamic Architecture of a Developing Organism: Dordrecht.—Kluwer Acad. Publ., 1998.—248 p.

10. Belousov L. V. Morphomechanics: goals, basic experiments and models // Int. J. Dev. Biol.—2006.— Vol. 50, № 2.—P. 81-92.

Статья поступила 5 июля 2012 г. Назаров Максим Николаевич

Национальный исследовательский университет МИЭТ, ассистент кафедры высшей математики РОССИЯ, 124498, Москва, Зеленоград, проезд 4806, д. 5 E-mail: [email protected]

THE SEARCH FOR UNIVERSAL PARAMETERS FOR THE MODELS OF MORPHOGENESIS

Nazarov M. N.

The parameters of models of morphogenesis can be considered universal in case they have biological meaning and their values can be freely exchanged between models of different classes. In the scope of this work exchange of chemical intercellular signals was chosen as the main mechanism controlling the activity of cells. For this case a set of universal parameters was obtained and it two models of morphogenesis, the discrete imitational and the average quantitative one, were built on this base.

Key words: morphogenesis modelling, intercellular signal exchange, coordination of different classes of models.