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

Вычислительные технологии

Том 11, № 2, 2006

ИССЛЕДОВАНИЕ ОДНОМЕРНОЙ МОДЕЛИ РЕГУЛЯЦИИ РАЗМЕРОВ ВОЗОБНОВИТЕЛЬНОЙ ЗОНЫ В БИОЛОГИЧЕСКОЙ ТКАНИ*

С.В. Николаев, Н.А. Колчлнов Институт цитологии и генетики СО РАН, Новосибирск, Россия

С.И. Фадеев, В.В. Когай Институт математики СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]

Э. МйолснЕсс School of Information and Computer Science, and Institute for Genomics and Bioinformatics, University of California, USA

A cellular-oriented model for control of the size of renewing zone in cellular ensemble is studied numerically. By the example of soot apical meristem, it was shown that the model has stable solutions that correctly describe the homestasis of the meristem structure.

Введение

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

По современным представлениям, в процессе развития в ткани возникают "островки" возобновительной ткани. Такие возобновительные зоны имеют следующую структуру: в центре зоны находится одна или несколько стволовых клеток, которые делятся и создают пул "переходных делящихся клеток" (transit amplifying cells). Переходные клетки делятся

* Работа выполнена при частичной финансовой поддержке NSF: FIBR (грант № EF-0330786), Государственным контрактом № 02.467.11.1005 с Федеральным агентством по науке и инновациям.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

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

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

Как упоминалось выше, на протяжении всей жизни растение поддерживает определенные группы клеток в недифференцированном состоянии. Главнейшей среди них является группа клеток на верхушке вегетативного побега растения, которая называется апикальной меристемой побега (АМП). Она содержит стволовые клетки, которые постоянно делятся и в конечном счете дают начало всем клеткам растения. Хотя клетки АМП не дифференцированы, они детерминированы в отношении экспрессии определенных генов, и на этом основании АМП подразделяют на несколько компартменов, находящихся в определенном пространственном расположении друг относительно друга на протяжении всей жизни растения [2, 3, 4]. Клетки, расположенные вокруг вертикальной оси меристемы (рис. 1) в радиусе 2-4 клетки в самых верхних 3-4 слоях, экпрессируют (т. е. "включают" соответствующие гены и в результате синтезируют) некоторый белок, называемый СЬУ3, и принадлежат центральной зоне (ЦЗ). В нижнем слое клеток ЦЗ располагаются клетки, экспрессирующие WUS. Эти клетки относят к организационному центру (ОЦ), толщина которого в вертикальном направлении может составлять 2-3 клетки. Такое постоянство структуры АМП, как полагают, необходимо для поддержания пула стволовых клеток [2, 3, 4]. Механизмы, обеспечивающие постоянство структуры АМП, являются предметом интенсивных исследований, как экспериментальных, так и теоретических.

Рис. 1. Поперечный срез апикальной меристемы побега АгаЫйорв1в Т.

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

На рис. 1 изображен поперечный срез апикальной меристемы побега АгаЫАорзгз Т [2]. Здесь Ь1 — наружный слой клеток, Ь2 — второй слой клеток, Ь3 — условно называется третьим слоем. Фактически в результате того, что клетки, находящиеся ниже второго слоя, делятся во всех плоскостях, это уже не слой, а скопление клеток, CZ — центральная зона, PZ — периферическая зона, — риб-зона, здесь клетки начинают дифференцироваться в клетки сосудистой системы растения. Ось х направлена вниз от верхушки побега. Клетки вдоль оси принимаются за одномерный массив рассматриваемой модели.

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

1. Одномерная модель структуры меристемы

В принципе возможны два способа поддержания вертикальной компартментализации меристемы: первый способ — симметричное деление клеток на границах компартментов с их детерминацией в поле морфогенов, второй способ — асимметричное деление клеток на границах компартментов [2, 3, 4]. Тот факт, что в слое Ь3 деление происходит во всех плоскостях, является аргументом в пользу того, что скорее всего регуляция структуры меристемы осуществляется первым способом. Кроме того, мутанты, у которых нарушен паттерн ориентации делений на ранних стадиях, формируют проростки с нормальным основным планом строения [5]. Поэтому в качестве возможного механизма мы будем рассматривать детерминацию клеток, управляемую позиционной информацией.

Физическим носителем позиционной информации являются поля концентраций веществ, распространяющихся из различных источников (например, путем диффузии). В самом простом случае размер некоторой "выделенной" зоны от "начала системы отсчета" может определяться пороговым значением концентрации некоторого диффундирующего из начала отсчета вещества У (например, зона активации Ш). Но если скорость синтеза вещества, "определяющего систему координат", изменится, то изменится и размер зоны.

В случае, когда при некоторой концентрации (довольно узкой) индуцируется "стабильный" источник Ш, его сигнал в начале координат может быть использован для "измерения" его расстояния от начала отсчета. Узкую зону индукции синтеза сигнала можно создать при наличии его ингибитора при больших концентрациях У. В такой системе уровень W в начале координат может интерпретироваться как расстояние от начала координат до его источника: чем слабее сигнал, тем дальше зона его синтеза. Если скорость синтеза У пропорциональна сигналу, то удаление сигнала будет приводить к снижению скорости синтеза У, а следовательно, к "приближению" зоны синтеза сигнала. Увеличение уров-

Зона активации С - ингибитора для \Л/

0.2

0.6

0.4

0.8

Y

0

Зона активации И/

Стационарное распределение У как функция расстояния от верхушки АМП

5 х

Зона экспрессии И/

Рис. 2. Схема одномерной модели. Распределение У в зависимости от расстояния от верхушки меристемы и его пороговые значения, при которых активируется экспрессия С и Ш.

ня сигнала в начале координат вызовет увеличение скорости синтеза У, что приведет к "отодвиганию" источника сигнала от начала координат.

На рис. 2 "начало отсчета" связано с зоной активации У. Здесь по оси х отложено расстояние в условных единицах вниз по вертикальной оси, проходящей через центр верхушки побега, по оси у — концентрация некоторого морфогена У, распространяющегося из верхушки побега (из точки 0). В результате диффузии У и его постоянного распада устанавливается некоторое стационарное распределение (убывающая от х функция). Вещество У при концентрациях выше пороговых может активировать экспрессию генов веществ С и Ш. Причем порог активации для С выше порога активации для Ш. Поскольку предполагается, что С является репрессором экспрессии гена Ш, там, где происходит экспрессия С, экспрессия Ш подавляется, а экспрессия гена Ш реально происходит в зоне, удаленной от верхушки побега (от начала координат).

2. Уравнения модели

Рассмотрим одномерный массив из п клеток, в котором между клетками может происходить перенос веществ У, Ш, Еу и — сопротивления, обратные переносу. В клетках могут происходить реакции синтеза веществ со скоростью, зависящей от присутствия других веществ:

Здесь Ец — коэффициенты чувствительности регуляции, которые больше нуля, если вещество ] стимулирует синтез вещества к, и меньше нуля — если угнетает; Тк — коэффициенты, обратные максимальной скорости экспрессии. Параметры Нк, как и Ец, определяют пороговые значения функции д(х):

Для описания скоростей реакций в предлагаемой модели используется сигмоидная функция [6]:

д(х) = 2 (1 +

X

V/ГT

X2

Далее, предполагается, что диффундируют вещества У и Ш. Вещество У синтезируется в клетке 1 и диффундирует в другие клетки массива. Скорость его синтеза зависит от концентрации вещества Ш в клетке 1. В зависимости от концентрации У в других клетках может синтезироваться вещество С, которое не диффундирует, а только распадается. В тех же клетках в зависимости от концентрации веществ У и C может синтезироваться вещество Ш, которое диффундирует по клеточному массиву и, достигая клетки 1, регулирует синтез вещества У.

В результате предлагаемая модель процесса формулируется как задача Коши для автономной системы уравнений следующего вида:

— = -аууг + Бу(у2 - уг) +--д(хг), хг = Ну + ЕуШщ

¿Ух

(г

(и аг

Г

Ту

= -ау уг + Ву (у— - 2уг + уг+г), г = 2, 3,... ,п - 1,

(Уп п ( ,

= -аууп + Бу(уп-1 - упУ,

¿Сг 1

= -агСг +--д(и), иг = Не + Есууг, г = 1, 2,...,п;

аг те

= -ашщ + Бы(щ2 - Щг) + — д(Уг), Уг = Нш + Е^уг + Е^Сг, аг тш

(щ- 1

= -атчщ + Бы(щг-г - 2щг + щг+г) +--д(У), У = К + ЕШууг + Е^Сг

аг ты

г = 2, 3,...,п - 1,

(щ 1

—- = -ашщп + Бш(щп-г - Ып) +--д(Уп), Уп = Нт + ЕШууп + ЕШеСп.

аг ты

Начальные условия при г = 0:

уг = у0, Сг = С0, щг = щ0, г = 1, 2,...,п.

(1)

(2)

(3)

(4)

Здесь уг, Сг, щг — концентрации веществ в г-й клетке; ау, ае, аы — коэффициенты распада; Б у, Бы — коэффициенты диффузии; определение параметров Ту, Те, тш , Не, Нт, Ееу, Еыу, Еые дано выше.

Исследование стационарных решений автономной системы (1)-(3) показало, что при определенном подборе параметров предложенная модель действительно описывает простой механизм экспрессии, при котором обеспечивается устойчивое положение максимума концентрации вещества Ш в той области пространства, где другое вещество У находится в определенной концентрации. При этом, если заданный уровень концентрации вещества У изменяет положение в пространстве, то соответственно перемещается максимум концентрации вещества Ш. Оказалось, что положение максимума устойчиво к возмущениям стационарных концентраций вещества Ш в достаточно широком диапазоне, т. е. модель отражает высокое "качество удержания зоны Ш".

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

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

dY „ „ 1 ^ dC „ 1 ^

-¡7 + Qy Y =— Fy, — + acC =- Fc, dt ту dt тс

dW 1

— + Qw W = — Fw, dt tw

(5)

Y = Y0, C = C0, W = W0 при t = 0.

Здесь У, С, Ш — векторы с компонентами у^, е^, г = 1, 2,... , п, соответственно; У0, С0, Ш0 — векторы начальных условий; Qy и — трехдиагональные матрицы:

Qy

ay + Dy J^y

-Dy

-Dy

2Dy -Dy

-Dy ay + 2Dy

D

Dy

y

Dy

Qw

aw + Dw Dw

Dw aw + 2Dw Dw

Dw aw + 2Dw Dw

Dw

обладающие диагональным преобладанием и, следовательно, хорошо обусловленные. Компонентами векторов , ^ являются сигмоидные функции с указанными в (1)-(3) аргументами:

Fy

" g(Xi) " g(Ui) g(Vi)

0 , Fc = g(U2) F = w g(V2)

0 g(Un-1) g(Vn-1)

0 g(Un ) g(Vn)

3. Стационарные решения

Рассмотрим численное определение стационарных решений автономной системы (5), обращаясь для этого непосредственно к системе нелинейных уравнений:

Qy Y = - Fy, acC = - Fc, Qw W = — Fw.

Ty

(6)

Напомним, что в (5) первая компонента вектора зависит от первой компоненты вектора Ш.

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

a

y

a

y

c

w

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

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

Рассмотрим систему (6) с параметром р > 0, задание которого приводит к уравнению

—д(Х 1) = р, Х-1 = Ну + ЕуШ иг.

(7)

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

1

^ (р)= рЯ-1е 1, С (р)

тгаг

-Ре(р), Ш(р)

- Я-^ (Р),

тт

где 61

первый столбец единичной матрицы;

д(и1(р))

Рс(р)

д(ВД)

д(ип-1(р))

Ръ (р)

дУ(р)) д(ВД)

д(Уп-1(р))

д(ип(р)) \ I д(Уп(р))

Щрр) = Не + ЕСуУг(р), Уг('р) = К + ЕъуУг(р) + ЕшсОг(р), г =1, 2,..., п.

В результате становится известным вектор Ш(р), а следовательно, его первая компонента и1(р). Отсюда значение параметра ту, соответствующее заданному значению параметра р, определяется по формуле

ту(р) = д(Х1(р)), Х1(р) = Ну + Еуади1(р).

р

(8)

Тем самым найдено решение нелинейной системы уравнений (6) при ту, вычисляемом по формуле (8).

График функции и1 = и1(ту), заданной параметрически в виде и1 = и1(р), ту = ту(р), р0 <р<р*, будем называть диаграммой стационарных решений системы (6). После построения диаграммы число стационарных решений определяется числом пересечений графика функции и1 = и1(ту) с прямой ту = ту, где ту — заданное значение параметра.

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

Устойчивость стационарных решений определяется численно в результате интегрирования задачи Коши (5) с начальными данными в виде "возмущенного" стационарного решения.

у

4. Полунеявный метод интегрирования

Для интегрирования автономной системы с постоянным шагом А воспользуемся простейшей полунеявной схемой точности порядка А. Введем обозначения:

г е Мт], ¿¿+1 = и + А, г = 1, 2,..., = 0,

У * « У (¿¿), С* « С (¿¿), Ш * « Ш (¿¿).

Приближенно автономную систему представим в виде разностных уравнений:

У*+1 _ У*

+ Qy Yi+1 = F

A

C i+1 C i 1

+ acCi+1 = - Fi, A tc

W^1 _ Wi . , ! 1 .

- + Qw Wi+1 = - FW

A

Tw

где

" g(Xi) " g(Uî)

0 , Fci = c g(U) , Fwi = w g(Vi)

0 -1) g(Vi-1)

0 g(Un) . g(Vi)

Fi

XI = йс + Ц, и] = йс + ЕСуу}, V/ = + у* + , = 1, 2,..., п. В результате получаем задачу Коши для следующей системы разностных уравнений:

А

[I + AQ„ ]Yi+1 = Yi +--AFy, Y1 = Y0,

Ty

(1 + Aac)Ci+1 = Ci + A Fc\ C1 = C0,

rc

[I + AQw ]Wi+1 = Wi + — FW, W1 = W'

(9)

где I — единичная матрица. Как следует из вида системы разностных уравнений (9), значения векторов У¿+1, С*+1, Ш¿+1 при г ^ 1 определяются из решения систем линейных алгебраических уравнений с матрицами [I + ], (1 + Аас)1 и [I + ] соответственно, т. е.

Yi+1 = [I + AQy ]-1(Yi + AFy), Ci+1

1

y

Wi+1 = [I + AQw]-1 (Wi + AFW^ .

1 + Aac

Ci + A F!

Tc

(10)

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

w

5. Примеры расчетов стационарных решений

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

ау = 0.1, Бу = 6, ту = 1, Ну = -5, Еуы = 40, ас = 1, тс = 1, Нс = -20.5, Есу = 20,

аш = 0.75, Б,ш = 1.5, = 1, = -30, ЕТОу = 60, Е.шс = -80.

(11)

На рис. 3 представлена диаграмма стационарных решений, из которой следует, что при ту = 1 система (6) определяет три стационарных решения.

На рис. 4-6 представлены распределения концентраций веществ в клетках с номерами 1-30: рис. 4 — первое решение, рис. 5 — второе решение, рис. 6 — третье решение (п — число клеток).

Исследование устойчивости показало, что первое и третье решения — асимптотически устойчивы, а второе решение — неустойчиво. На рис. 7 представлен переходной процесс, в котором неустойчивое второе решение стремится к асимптотически устойчивому первому решению.

Рис. 3. Диаграмма стационарных решений. Рис. 4. Первое стационарное решение при ту = 1.

Рис. 5. Второе стационарное решение при ту = 1. Рис. 6. Третье стационарное решение при ту = 1.

Рис. 7. Потеря устойчивости второго стацио- Рис. 8. Пятое стационарное решение при нарного решения. Ту = 0.6.

Рис. 9. Третье стационарное решение при Рис. 10. Потеря устойчивости третьего стаци-Ту = 0.6. онарного решения.

При ту = 0.6, как следует из диаграммы стационарных решений, система (6) определяет пять стационарных решений. Из них решения первое и пятое — асимптотически устойчивы (рис. 8), а второе, третье и четвертое решения — неустойчивы. При этом второе решение, потеряв устойчивость, стремится к первому решению, а четвертое решение — к пятому решению. Третье решение (рис. 9) после потери устойчивости переходит в устойчивые автоколебания (рис. 10).

6. Модель с непрерывным распределением веществ

Формально систему (1)-(3) можно рассматривать как результат дискретизации следующей системы уравнений, описывающей непрерывное распределение веществ у(£, г), с(Ь, г) и г):

= ПуЙ - ауу, 1> 0, г е [0,Я], Я> 0; (12)

де £(и) тт ^ ^

— = — асе + —^, и = йс + ЕСу у; дт д2т )

"д£ = "дТ^ — т + ——, V = + ЕадуУ + Еадсе.

(13)

(14)

Параметры Ду, , ау, ас, имеют тот же смысл, что и в модели (1)-(3). Значение К достаточно велико. Экспрессия у, определяемая значением 0), задается граничным условием на левом конце отрезка по г:

ду

1

дг = — -2(Х), X = йу + Еуад0).

(15)

На правом конце отрезка по г ставится условие

ду

— = 0 при г = К.

дг

Краевые условия для уравнения (14) имеют вид

дт л

-7— = 0 при г = 0 и г = К. дг

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

Рассмотрим краевую задачу на отрезке [0, К], описывающую стационарное распределение веществ в виде графиков функций у(г), е(г) и т(г):

—2у

£у — ау у = 0, г € [0, К], Ду -¡у =--^(Х^) при г = 0,

—г

Ту

(16)

ду

— = 0 при г

дг

К;

—2т 1

^ — т + — ) = 0, г € [0,К],

—г2 Гад

—— = 0 при г = 0, —г

—— = 0 при г = К, —г

(17)

где

1

Гса,

-£(и), и = йс + ЕСу у,

X! = йу + Еуад т(0), V = + Еаду у + Е^С.

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

Как и ранее, следуя идее параметризации, будем считать, что в уравнении

—#(Х1) = р, Х1 = йу + Еуад т(0),

(18)

е

у

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

р 1 + е-2шу (Я-т) Гл~

у(г,р) = 1+ ^ е~Шуг, ^у = \ . (19)

^/аущ 1 - е-2^ " ' ^ ув

Отсюда следует, что

с(т,р) = —д(и (г,р)), и (г,р) = Не + Ееу у(г,р). (20)

8еТе

Решение линейной краевой задачи относительно щ(т) можно представить в интегральном виде, используя выражения функций Грина:

-ш-(т-з) (1 + е-2ш™(д-г))(1 + е-2ш™в)

= , Ki(r,s) =---—2—R-, r>s

Dw ' ' 2\Jaw Dw 1 - e

K2(r, s) = —. ------, r < s.

A ' J 2VÖD 1 - e-2owR '

Для удобства введем обозначения

С(т,'р) = — д(У (г,р)), У (г,р) = К + ЕШу у (г, р) + ЕтеС(т,р). При этом решение принимает вид

т Я

щ(т,р) = ^ Кг(т,8)С(8}р)(18 ^У К2(т,в)С(в,р)ав. (21)

0т

Отсюда, в частности, имеем

Я

/е-ш-з 1 + е-2ш™ (я-з)

К2(0, 8)0(8,р)ав, К2(0, 8) = -== 1+—.

V а-ш Б'Ш 1 е -

0

Таким образом, заданному значению параметра р соответствует значение параметра ту (р), вычисляемое по формуле

Ту(р) = д(ХЛр)), Хг(р) = Ну + Еуыщ(0,р). (22)

р

Тем самым найдено решение нелинейной краевой задачи (16), (17) при ту = Ту(р).

Отметим, что при К ^ то решение краевой задачи (16), (17) определяется по формулам (18)-(22), где

у(т,р) = ,р =е-ШуТ;

ауБу

e-Ш-w (r-s) e-LOw (s-r)

K\(r, s) = —= (1 + e~2üJws), K2(r, s) = —= (1 + e~2üJwr).

2\j aw Dw 2\j aw Dw

Рассмотрим дискретный аналог краевой задачи (16), (17). С этой целью введем равномерное разбиение отрезка [0, Л] на п — 1 частей с узлами гг:

гг = к(г — 1), к

Л

п1

г = 1, 2,..., п.

Обозначим через уг, сг, иг приближенные сеточные значения у(гг), с(гг), и(гг) соответственно. Используя вместо производных по г их разностные аналоги

дУ , . _ Уг — Уг-1

дГ(Гг)--^

д2У { \ ~ Уг-1 — 2Уг + У+1

(Гг) ~ -

дг2

к2

ди иг — иг-1

"зГ(Гг) * к ,

д2и иг-1 — 2иг + иг+1

-(гг) « -

дг2

к2

получим уравнения дискретной модели типа (6):

ПуУ = , асС = - П« W = — ^,

где У, С, W, , — векторы с теми же компонентами, что и в (6); Пу и П«

трехдиагональные матрицы, зависящие от параметров ау, а«, 9у и :

(23)

9у

к2

9«

к2

ау + 5у —5у

—5у ау + 25у —5у

—5у ау + 25у

-Зу а

9у ■у Ь ^у

П«

а« +

—

——

и,« +

— /

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

В качестве примера рассмотрим решение краевой задачи (16), (17) с параметрами (11), где вместо п задано Л = 30. На рис. 11 представлена диаграмма стационарных решений, из которой следует, что при ту = 1 краевая задача имеет три решения. Различие в диаграммах на рис. 3 и рис. 11 связано с тем, что формально система (6), представляющая дискретную модель краевой задачи, недостаточно хорошо ее аппроксимирует. Однако это не имеет "практического" значения, поскольку сам процесс моделируется решениями системы (6).

Решение 3 приведено на рис. 12. Сопоставление с рис. 6 показывает, что концентрации веществ в г-й клетке определяются сеточными значениями краевой задачи в г-м узле сетки при равномерном разбиении отрезка [0, 30] по г на 30 частей.

W,

1.4-1

0.00

0.50

0.25

0.75

1.25

1.00

0.0 0.2 0.4 0.6 0.8 l!0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 т.

•у

0

5

10

15

20

25

30 Г

Графики, представленные на рис. 11 и 12, практически останутся теми же, если решение краевой задачи определить на большем отрезке, К > 30. В этом случае концентрации веществ в г-й клетке будут определяться сеточными значениями решения краевой задачи в г-м узле г г, г = 1, 2,30. При К > 30 концентрации веществ монотонно стремятся к нулю.

Заключение

Рассмотренная в статье модель показывает, что положенный в ее основу механизм взаимной регуляции детерминации пространственно распределенных клеток может обеспечить стабильную пространственную локализацию всех зон моделируемого клеточного ансамбля. Это проявляется в том, что возмущения, вносимые в пространственное распределение веществ У, С, W, в свою очередь, возмущают "правильное" распределение зон в клеточном ансамбле. Предложенный механизм приводит к "исправлению" таких возмущений и тем самым к стабилизации пространственной локализации зон. Это хорошо согласуется с тем фактом, что размеры и расположение компартментов апикальной меристемы побега остаются стабильными на протяжении всей жизни растения, несмотря на постоянные возмущения, которые апикальная меристема испытывает со стороны как среды, так и остального растения.

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

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

Авторы выражают благодарность Н.А. Омельянчук, Н.Л. Подколодному, В.В. Мироновой за весьма полезные дискуссии.

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

[1] Merks R.M.H., Glazier J.A. A cell-centered approach to developmental biology // Physica A. 2005. Vol. 352. P. 113-130.

[2] Gross-Hardt R., Laux T. Stem cell regulation in the shoot meristem // J. of Cell Sci. 2003. Vol. 116. P. 1659-1666.

[3] Bowman J.L., Eshed Y. Formation and maintenance of the shoot apical meristem // Trends in Plant Sci. March 2000. Vol. 5, N 3.

[4] Sharma V.K., Carles C., Fletcher J.C. Maintenance of stem cell populations in plants // PNAS. Sept. 30. 2003. Vol. 100. P. 11823-11829.

[5] Berleth T., Chatfield S. Embryogenesis: Pattern Formation from a Single Cell. The Arabidopsis Book. ©2002 American Society of Plant Biologists. (http://www.bioone.org/pdfserv/i1543-8120-007-01-0001.pdf)

[6] Mjolsness E., Sharp D.H., Reinitz J. A connectionist model of development // J. of Theoretical Biology. 1991. Vol. 152. P. 429-454.

[7] Фадеев С.И. О решении системы трансцендентных уравнений с параметром методом Ньютона // Сплайн-аппроксимация и численный анализ. Вычислительные системы. Новосибирск, 1985. Вып. 108. С. 78-93.

[8] Фадеев С.И. Программа численного решения нелинейных краевых задач для систем обыкновенных дифференциальных уравнений с параметром // Вычислительные методы линейной алгебры / Под ред. С.К. Годунова. Новосибирск: Наука, Сиб. отд-ние. (Тр. Ин-та математики) 1990. T. 17. С. 104-198.

[9] Fadeev S.I. Organization of numerical experiment for investigation of nonlinear boundary value problems by the method of continuation of solution with respect to parameter // Sib. J. Diff. Equation. 1998. Vol. 1, N 4. P. 321-350. N.Y.: Nova Sci. Publ., Inc., 1998.

Поступила в редакцию 22 декабря 2005 г.