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

ФИЗИКА

УДК 517.958:57

Брацун Дмитрий Анатольевич

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

Захаров Андрей Павлович

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

ФГБОУ ВПО «Пермский государственный гуманитарно-педагогический

университет», Пермь, Россия 614990, Пермь, Сибирская, 24, (342) 238-63-64, e-mail: dmitribratsun@rambler.ru

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РАКОВЫХ ОБРАЗОВАНИЙ ПРИ КОЛЛЕКТИВНОМ ВЗАИМОДЕЙСТВИИ КЛЕТОК ЭПИТЕЛИЯ

Dmitry A. Bratsun

DS, Head of the Theoretical Physics Department Andrej P. Zakharov PhD, Lecture assistant of the Theoretical Physics Department

Federal State Budget Educational Institution of Higher Professional Education «Perm State Humanitarian Pedagogical University»

24, Sibirskaja, 614990, Perm, Russia, e-mail: dmitribratsun@rambler.ru

MATHEMATICAL MODELLING CANCERS IN THE COLLECTIVE INTERACTION OF EPITHELIUM CELLS

Аннотация: предлагается математическая модель возникновения

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

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

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

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

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

Abstract: In this paper, we propose a mathematical model of cancer tumour occurrence in a quasi epithelial tissue. Basic model of the epithelium growth describes the appearance of intensive movement and growth of tissue when it is damaged. The model includes the effects of division of cells and intercalation. It is assumed that the movement of cells is caused by the wave of mitogen-activated protein kinase (MAPK), which in turn activated by the chemo-mechanical signal propagating along tissue due to its local damage. In this paper it is assumed that cancer cells arise from local failure of spatial synchronization of circadian rhythms. It is assumed that the study of the dynamic properties of the model could determine the spatial relationship between the occurrence of cancer cells and development of the entire tissue parameters coordinating its evolution through the exchange of chemical and mechanical signals.

Key words: mathematical modeling, cancer growth, circadian rhythms, gene regulation, synchronization, complex systems.

Моделирование процессов возникновения и роста злокачественных опухолей, несомненно, является одним из магистральных направлений в общем ряду математического моделирования в биологии. На данный момент существует обширная литература по этому вопросу, систематизация которой время от времени появляется в обзорах и монографиях [3]. Как отмечается автором, одной из главных проблем математического моделирования развития рака (как, впрочем, и вообще биологических систем) является тот факт, что разворачивающиеся при этом процессы являются разномасштабными. С одной стороны, они включают процессы генной регуляции, протекающие в масштабах клеточного ядра. Именно с этого уровня может прийти сигнал, который заставляет перерождаться здоровую клетку в раковую. На этом масштабном уровне может произойти сбой в «программируемой» смерти клетки, и она может перейти к неупорядоченному делению. С другой стороны, процессы межклеточного взаимодействия (масштаб клетки или группы клеток) также являются важными для понимания момента зарождения опухоли. Если иммунная система организма действует правильно, раковые клетки могут распознаваться и уничтожаться. Если происходит сбой, то повышается риск

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

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

Другой тип моделирования - когда в дополнение к популяционной динамике в рассмотрение вводятся определенные переменные, определяющие усредненную структуру популяции раковых клеток (например, возраст клеток),

- можно отнести к полуфеноменологическим моделям [6,7]. Хотя такие модели относят к более развитому классу, они имеют те же недостатки, что и модели выше. Существует ряд работ, посвященных изучению вопроса связи межклеточных взаимодействий и макроскопического уровня самой опухоли. Самый распространенный подход здесь - модели, в которых опухоль рассматривается как сплошная среда. Как правило, система уравнений в частных производных включает в себя уравнение сохранения массы для «клеточной среды» и уравнение «реакция-диффузия», описывающее поле химических сигналов, которыми обмениваются клетки [8,9]. Такие модели можно отнести к феноменологическим моделям сплошных сред со всеми недостатками, которые им присущи. Например, в работе [8] опухоль рассматривается в виде твердой матрицы пористой среды, которая взаимодействует с насыщающей ее «клеточной жидкостью» из здоровых клеток. Таким образом, авторам удается описать некоторые пространственные черты реальных опухолей (фрактальность структуры).

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

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

Имеется большая группа моделей, в которых особое внимание уделяется взаимодействию на микроскопическом и клеточном уровне. Например, в работе [13] рассматривается вопрос, как слишком сильная экспрессия определенных генов может привести к сбою в работе клетки, т.е. сделать ее раковой. При этом может быть использовано стохастическое описание процессов генной регуляции [14]. Слабым местом этого моделирования является то, что за скобками рассмотрения остается процесс формирования самой опухоли. Учитывая, что такое заболевание, как рак, - явление разномасштабное, наиболее реалистичный подход к моделированию требует учета в модели процессов на всех уровнях описания. В силу трудности этого подхода в литературе имеется не так много попыток подобного анализа. Одной из первых в осуществлении этого подхода обычно называют работу [15]. В которой, модель представляет собой гибридный подход, включающий расчет клеточных автоматов, состояние которых определяется непрерывным распределением кислорода вокруг кровяного сосуда вблизи возникновения опухоли. Еще одну примечательную попытку представляет работа [16]. Авторы рассмотрели сферически растущую опухоль, которая включала в себя латтис-модель дискретных клеток. Для каждой клетки отдельно производился расчет системы однородных дифференциальных уравнений, описывающих процессы генной регуляции. Клетки могли механически взаимодействовать между собой в рамках решеточного газа. Из-за большого числа клеток (более миллиона) попытка оказалась не очень удачной.

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

Разработка таких моделей

сложных систем (individual based models) является нетривиальной задачей. Развитие этого

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

исследователей к возможности полноценного реалистичного

моделирования функционирования живой ткани.

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

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

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

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

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

Рис. 1. Элементы хемомеханической модели эпителия, предложенной в работе [15]

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

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

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

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

dm aH(o)o _ /1Ч

^ — = С— Pmm , (1)

dt 1+bo

где m - локальная концентрация MAPK в клетке; o = S/S0-1 - нормированное

отклонение площади клетки от среднего размера клетки в сообществе; H -функция Хэвисайда; pm - скорость деградации MAPK; т0, a и b - параметры задачи. Важную роль в активации MAPK играет концентрация сигнального белка С, который может распространяться от клетки к клетке:

= X Jik - РсСк , (2)

dt ieadj(k)

где Рс - скорость деградации сигнального белка; сумма вычисляется по соседним клеткам, а межклеточный поток Jik белка С из клетки i в клетку k вычисляется как

Jik = aLik(Ci -Ck) . (3)

Здесь a - коэффициент диффузии; Lik - длина общей границы между клетками i и k (см. рис.1). Химический сигнал в клеточной ткани может распространяться вследствие её повреждения (порез, рана). Уравнения (1-3)

показывают, каким образом сигнальный белок C активирует внутри клетки МАРК. В свою очередь протеинкиназа ответственна за формирование так называемой активной силы

¥ас‘ = (т. Рк) , (4)

1 \ к к / кеа^О) ^ ’

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

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

Подмодель циркадианных колебаний. Как известно, злокачественные опухоли имеют множество типов и могут возникать благодаря очень большому числу причин. Среди них называются обычно такие мутагенные процессы, как радиация, воздействие токсических веществ и т.д. В последние годы, однако, стала проявлять себя тенденция более глубокого взгляда на причины возникновения рака. Например, в ряде недавних статей [18,19] утверждается, что до 30 % заболеваний могут быть связаны с нарушением суточного ритма в организме. По крайней мере этот факт твердо установлен для людей, чья профессия требует нарушения естественного суточного ритма. Используя этот факт, наша модель возникновения раковых клеток основывается на сбое процессов синхронизации циркадианных колебаний в живой ткани. Для этого установим для каждой клетки генератор циркадианных ритмов и исследуем вопрос о пространственной синхронизации этих биоритмов на уровне всего ансамбля клеток растущего эпителия. В качестве механизма колебаний будем использовать динамическую модель циркадианных колебаний, предложенную авторами ранее [1]. Несмотря на то что модель предлагалась для описания суточных ритмов организма Neurospora crassa, она имеет достаточно общий характер и может быть использована для описания биоритмов других организмов. Главным элементом механизма колебаний является эффект запаздывания реакций синтеза белков в процессах транскрипции и трансляции генов (рис. 2). Эти процессы не просто медленные, но еще и состоят из

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

На рис.2 представлена схема взаимодействия двух белков, формально обозначенных как F и W, кодируемых двумя генами f и w.

У ряда организмов были выделены такие пары принципиальных генов, отвечающих за работу механизма циркадианных колебаний. Например, в случае нейроспоры это гены frq и wcc, у дрозофилы - per и dclock. Это не значит, что механизм колебаний поддерживается работой

исключительно этой пары генов, так как циркадианные ритмы

Рис. 2. Схема взаимодействий белков в модели циркадианных колебаний. Здесь:/ м> _ гены;

^ соответствующие им белки

обладают широким набором свойств (автономность, способность компенсировать температурные изменения среды, а также изменять фазы колебаний под действием внешнего освещения или температуры и т.д.), задаваемых десятками генов. Однако указанные пары являются принципиальными для поддержания ритмов. Как видно из Рис. 2, обобщенная модель включает в себя как положительную, так и отрицательную петлю обратной связи и является симметричной по отношению к обеим переменным. Кроме того, она учитывает процессы димеризации и деградации белков. Полный список реакций, происходящих при транскрипции генов, приведен в табл. 1. Концентрации белков-мономеров в таблице обозначены ^ и Ж, димеров _ ^2 и Ж2. Динамика оператор-сайта промотора описывается бинарной функцией В е {В0, В1}, которая принимает значение В0 в случае его открытия и В1 в случае закрытия. Наиболее важными для поддержания колебаний являются реакции синтеза белков, протекающие с характерным временем запаздывания т. Именно этот механизм определяет динамику системы и её устойчивость к внешним и внутренним возмущениям.

Таблица 1

Список реакций транскрипции генов

1 Процессы димеризации: f + f k1 > f2 W + W kW > W2

2 Процессы редимеризации: f2 kF > f + f W2 kW > W + W

3 Динамика оператор-сайтов: bFF + w2 kF > b1F bF—bFF + w2 bW + f2 kW > bW bW kW > bWW + f2

4 Процессы синтеза белков: bF (t) kr > bF + Ft+Tr bW (t) —w> bW + wt+Tw

Окончание табл. 1

5 Процессы деградации белков: р — > 0 Ж Уж > 0

6 Процесс образования гетеродимера: р + Ж—> 0

Примечание: здесь, к, к], к-и к2, к-2, ур-, уЖ - скорости соответствующих реакций.

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

(1 + 4КрР) ^ = кр

ш

1 -

1 + К^КрЖ2(1-г)

- урр - крЖ

(1 + 4КЖЖ)

ёЖ

= к

ж

1-

1 + КрКЖ р 2(t -г)

-уЖЖ-крЖ,

(5)

(6)

где Кр = кр / кр, КЖ = кЖ / кЖ. В численных расчетах были использованы значения параметров, приведенные в табл. 2, которые позволяют получить колебания с периодом 22,65 часа.

Таблица 2

Параметры модели

г к кр кЖ кр кр Кшх КЖ ур Уж

6 h 30 пМ1К1 8 пМ^ 4 пМ^ 5 пМ-1 5 пМ-1 5 пМ-1 5 пМ-1 0.3 И1 0.4 И1

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

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

С , Л

1

(1+4КрР)

кр -

кр

1+К^КрЖ,2($-г)

-гРр, - крщ

+аЕЦк(рк -р1). (7)

Здесь механизм распространения записан по аналогии с формулами (2,3). Коэффициент переноса белка р через клеточную мембрану равен а = 0,1.

1

1

к

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

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

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

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

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

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

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

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

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

Ф, = (к -4ке«ф(о , (8)

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

Если клетка , находится в полностью синхронизированном поле, то величина (8), очевидно, будет равняться нулю, так как фаза колебаний клетки не будет отличаться от фазы колебаний соседней. С другой стороны, значительной величины расфазировка будет достигаться на границе смены кластеров (рис.4). Тогда как максимальное значение расфазировки будет возникать в случае появления обособленной клетки внутри кластера, фаза циркадианных ритмов которой отлична от фаз соседних клеток (рис.5). Именно такие клетки, согласно нашей гипотезе, подвергаются наибольшей опасности дифференцироваться в раковое состояние. Уравнение состояния каждой клетки будем описывать с помощью простой феноменологической модели:

= -лх,(1 -X,Ха -X,)+NФ,£(О, (9)

М

где X - функция состояния , -ой клетки; Л - скорость затухания возмущения; N - амплитуда шума; ) - генератор случайных чисел на отрезке от 0 до 1. Динамическое уравнение (9) устроено так, что в отсутствие шума у него имеется три положения равновесия: 0, А и 1 (число А подбирается между 0 и 1). Крайние положения равновесия устойчивы и притягивают к себе все траектории из соответствующей области притяжения, ограниченной неустойчивым равновесием А. Будем считать, что состояние клетки X = 0 соответствует здоровой клетке, а состояние X = 1 - раковой. Параметры модели (9) калибруются таким образом, чтобы переход клетки в раковое состояние

осуществлялся под действием шума, но только при высоком значении параметра расфазировки (8). В расчетах использовались значения Х = 10, A = 0,15, N = 2,5, которые устанавливают, что переход клетки в раковое состояние происходит достаточно редко, но если он случился, то обратный переход принципиально не возможен.

Отметим еще одно важное обстоятельство, которое учитывается при построении модели (9). Хотя большинство динамических процессов в природе демонстрируют стохастическое поведение, по-настоящему роль флуктуаций имеет смысл учитывать только в случаях, когда количество степеней свободы у рассматриваемой системы невелико. Взрывной рост интереса исследователей к малоразмерным стохастическим системам возник в начале 2000-х во многом благодаря тому, что количество молекул, вступающих в реакции во время генных процессов транскрипции/трансляции, весьма невелико - скорость работы РНК-полимеразы, например, всего около 50 нуклеотидов в секунду. Таким образом, даже небольшие флуктуации концентраций рибонуклеиновой кислоты и белка могут иметь значительное влияние на общую динамику системы. В этом состоит принципиальное отличие этих систем от, скажем, гидродинамических систем с шумом, в эволюции которых принимают участие настолько огромное количество молекул, что флуктуации могут проявить себя только вблизи точки бифуркации. В настоящий момент существует уже большой список литературы, посвященный экспериментальным доказательствам, что шум является важнейшим регулятором генных процессов. Отметим, что в генетике принято различать два вида шума, которые могут возникать в системе. Первый вид шума связан со случайной природой самих химических реакций между небольшим количеством молекул внутри одной клетки (в англоязычной литературе используют термин «внутренний шум» -intrinsic noise). Все, что в этих условиях работает на усиление стохастичности, может быть отнесено к этому типу шума. Например, флуктуации, возникающие при нагреве системы, определенно служат усилению внутренне присущего системе шума. Другой тип шума, получивший название «внешний шум» (extrinsic noise), генерируется за пределами клетки и связан с межклеточными различиями. Феноменологическая модель (9), очевидно, позволяет учесть влияние как внешнего шума (за счет межклеточных различий), так и внутреннего шума (для каждой клетки производится расчет своего собственного стохастического уравнения).

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

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

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

Калибровка параметров модели Письмена-Салма

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

отличаются от большинства клеток сообщества. Раковые

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

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

Рис. 4. Результаты расчета возникновения и эволюции раковых опухолей в ткани эпителия. Значения шкалы соответствуют изменениям размера клетки относительно нормального

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

1. Брацун Д.А., Захаров А.П. Моделирование пространственно-временной динамики циркадианных ритмов Neurospora crassa // КиМ. 2011. - Т. 3, № 2. - С. 191-213.

2. Захаров А.П., Брацун Д.А. Синхронизация циркадианных ритмов в масштабах гена, клетки и всего организма // КиМ. 2013. - Т.5, № 2. -С. 255-270.

3. Alarcon T. A cellular automaton model for tumour growth in inhomogenous environment // J. Theor. Biol. 2003. - Vol.225. - P.257-274.

4. Ambrosi D., Preziosi L. On the closure of mass balance models for tumour growth // Math. Mod. Meth. Appl. Sci. 2002. - Vol.12. - P.737-754.

5. Anderson A. A hybrid mathematical model of solid tumour invasion: The importance of cell adhesion // Math. Med. Biol. 2005. - Vol.22. - P.163-186.

6. Anderson A., Weaver A., Cummings P., Quaranta V.. Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. // Cell. 2006. - Vol.127. - P.905-915.

7. Byrne H.M., King J.R., McElwain D.L.S. A two-phase model of solid tumor

growth // Appl. Math. Lett. 2003. - Vol.16. - P. 567-573.

8. Dyson J., Villella-Bressan R., Webb G. The steady state of a maturity

structured tumor cord cell population // Discr. Cont. Dyn. Syst. B. 2004. -

Vol.4. - P.115-134.

9. Greene M.W. Circadian rhythms and tumor growth // Cancer Letters. 2012. -Vol. 318. -P. 115-123.

10. Gyllenberg M., Webb G. A nonlinear structured population model of tumour growth with quiescence // J. Math. Biol. 1990. - Vol.28. - P. 671-684.

11. Kim Y., Stolarska M.A., Othmer H.G. A hybrid model for tumour spheroid growth in vitro I: Theoretical development and early results // Math. Mod. Meth. Appl. Sci. 2007. - Vol.17. - P.1773-1798.

12. Kimmel M., Lachowicz M., Swierniak A. Cancer growth and progression, mathematical problems and computer simulations // Int. J. Appl. Math. Comput. Sci. - 2003. - Vol.13. - P.279-429.

13. Komarova N. Stochastic modeling of loss- and gain-of-function mutation in cancer // Math. Mod. Meth. Appl. Sci. 2007. - Vol.17. - P. 1647-1674.

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

15. Rossetti S., Esposito J., Corlazzoli F., Gregorski A., Sacchi N. Entrainment of breast (cancer) epithelial cells detects distinct circadian oscillation patterns for clock and hormone receptor genes. // Cell Cycle. 2012. - Vol. 11. - P.350-360.

16. Salm M., Pismen L.M. Chemical and mechanical signaling in epithelial spreading // Phys. Biol. 2012. - VOL. 9, N.2. - P.026009-026023.

17. Smolle J., Stettner H. Computer simulation of tumour cell invasion by a stochastic growth model // J. Math. Biol. 1993. - Vol.160. - P.63-72.

18. Vogelstein B., Kinzler K.W. Cancer genes and the pathways they control // Nature Med. 2004. - Vol.10. - P.789-799.

19. Webb G.F. Theory of Nonlinear Age-Dependent Population Dynamics. -Marcel Dekker, 1985 - 245 P.

20. Weber G.F. Molecular Mechanisms of Cancer. - Springer, 2007. - 645 P.