Популярное

Мифы о звукоизоляции



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



Как построить лестницы на садовом участке



Подбираем краску для ремонта



Каркасные дома из дерева


Главная » Астанин С.А.

Астанин С.А. (s.astanin@gmail.com) (1), Прециоси Л. (2), Тозин А. (2)

(1) МФТИ,

(2) Dipartimento di Matematica, Politecnico di Torino, Италия

1. Введение

Ожидается, что рак скоро выйдет на первое место среди причин смертности во многих странах мира, оттеснив на второе место заболевания сердечно-сосудистой системы. Масштаб проблемы ярко характеризует статистика заболеваемости по Великобритании, которая позволяет предположить, что в течение жизни заболевание разовьётся у каждого третьего жителя [1].

Математическому моделированию опухолевого роста посвящено большое количество работ. Обзоры современных математических моделей приведены в [2,3]. В последние годы появился целый ряд работ на русском языке, описывающих качественные особенности опухолевого роста в среде с гетерогенным распределением питательных веществ [4-6].

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

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

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

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

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

Модель роста опухолевого тяжа вдоль источника питательных веществ



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

2. Описание модели

Вводится система координат {x, y} в R2 такая, что прямая x = 0 совпадает с осью сосуда. Рассматриваются следующие функции: nT = nT(x, y, t) - объёмная доля клеточной фазы, vT = vT(x, y, t) - поле скоростей клеточной фазы, ni = n/(x, y, t) - объёмная доля жидкой фазы, vi = vi(x, y, t) - поле скоростей жидкой фазы, nECM -постоянная объёмная доля жёсткого внеклеточного матрикса. Условие насыщенности смеси nT +nl +nECM = 1 может быть также записано как nT + nl = nmax = const.

Уравнения баланса массы и напряжений

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

- div (nT vT ) = -, (1)

l - divinv,), (2)

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

div(nTvT + n,vl) = 0 . (3)

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

-div(nTTT ) = maT , (4)

-div(niTi) = Щ . (5)

Здесь TT и T, - тензоры напряжений Коши для клеточной и жидкой фазы соответственно, а m° и m t - силы воздействующие во время движения на каждый из компонентов со стороны среды. Предполагается следующий вид тензоров напряжений:

Tt =-( p + St ))E, (6)

T =-P E , (7)

где p = p(x, y, t) - давление в ткани, а Ът (nT) - функция, описывающая межклеточные взаимодействия, принимающая положительные значения при сжатии клеток. В зависимости от типа описываемых клеток, напряжения могут возникать не только при сжатии клеточной массы, но и при её растяжении (вплоть до разрыва межклеточных связей). Поэтому функция Ът (nT) может менять знак на (0; nmax), при этом точка перемены знака n0 - плотность состояния равновесия1.

1 Можно также рассматривать такие непрерывные функции 2T (nT), что 2T (nT) = 0 , при nT < n0. Эти функции описывают рыхлую ткань, без межклеточных связей.



Соответственно воздействия на компоненты смеси со стороны среды определяются следующими соотношениями:

тт = pVnT - ATvT + ATl(v; - vT), (8)

ml = pVn, -A (v, - vt ). (9)

Модель учитывает прикрепление клеток к матриксу вязким трением при их относительном движении, оно описывается слагаемым -ATvT. По сравнению с ним вязкое трение между жидкой фазой и матриксом считается пренебрежимо малым. Взаимодействие клеток и жидкости представлено слагаемыми +A (v, - vT).

Предполагая, что ATl = n72/k , где к - проницаемость среды [13], из (5), (7) и (9) получаем связь между скоростями жидкой и клеточной фазы:

V = Vt - - Vp , (10)

что совпадает по форме с законом фильтрации Дарси [14]. Аналогично, из (4), (6), (8) и (10), с учётом того, что n = nmax - nT , следует

VT =-A- (ax Vp + V(nT ЪТ (nT )) ) .

Подставляя полученные выражения для скоростей компонентов смеси в (3), получаем уравнение, выражающее связь между переменными nT и p:

(f nmax +1$ \]p + -Al V(nT Ът (nT) )) = 0. (11)

Подставляя выражение для vT в (1) получаем второе уравнение с теми же переменными:

Sb-. - div (nmax Vp +V(nT К )) )) .(12)

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

среде [15].

Чтобы определить скорость роста/гибели опухолевых клеток Гт , требуется учитывать не только структуру ткани, но и распределение в ней питательных веществ.

Упрощение модели

Модель упрощается, исходя из предположения, что взаимодействия при относительном движении клеток и жидкости малы (членом +ATl (v - vT) в (8), (9) можно пренебречь). Следствием этого является Vp = 0 . Соответственно, (11) и (12) сводятся к одному уравнению для nT:

dt \A.

div [Al- V(nT Ът (nT ) )\=Гт- . (13)

Перенос и потребление питательных веществ

Центральный кровеносный сосуд, вокруг которого развивается тяж, снабжает клетки питательными веществами, необходимыми для их выживания и размножения. Вещества,



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

Распределение питательного вещества определяется следующим диффузионным процессом:

- = DAc-anTc ,

dt T

где c - концентрация питательного вещества, D - коэффициент диффузии питательного вещества, a - скорость потребления питательного вещества клетками опухоли. Концентрация питательного вещества внутри сосуда предполагается постоянной.

Лимитирующим питательным веществом в модели считается кислород. В этом случае производной по времени dc/dt можно пренебречь ввиду высокой подвижности растворённого в крови кислорода по сравнению с подвижностью опухолевых клеток. Таким образом, распределение питательного вещества с высокой точностью описывается решением уравнения Пуассона

DAc =a nTc . (14)

Пролиферация и гибель клеток

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

В модели учитывается

гибель клеток при концентрации питательного вещества ниже некоторого порогового значения,

пролиферация клеток при концентрации питательного вещества выше некоторого порогового значения,

существование предельной плотности клеток опухоли nmax, при которой дальнейшая локальная пролиферация невозможна.

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

Гт (С, nT ) =У nT (c - cx )(nmax - nT )

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

3. Граничные и начальные условия

Рассматривалась выпуклая связная область G в {x, y x > 0 и y > 0}, примыкающая к осям x = 0 и y = 0. Соответствующие границы области обозначались dGx=0 , dGy=0 и dGQE .

Начальное состояние определялось так:

I n0, если x е G, ., ч

пТ(x,0) =\ 0 (15)

[ 0, если x <£ G.

Свободная граница dGсв совпадала со границей проникновения клеточной фазы и двигалась со скоростью vT . Накладывались следующие граничные условия на nT :



nT (x, t) = n0, x 6 5Gсв, (16)

(x,t) = 0, x6dGx==0 u5Gy=0, (17)

где n - внешняя нормаль к области G.

Последнее условие означает непрерывность напряжений на внешней границе dG . В упрощённой модели условие непосредственно следует из (4), (6) и Vp = 0, если отсутствуют внешние силы, воздействующие на поверхности опухоли.

Заметим, что (16) означает разрывность потока

I =-- V(nT Ът (nT ) )

т

на границе dG, если граница неподвижна. Именно поэтому условия (16), (17) должны быть дополнены уравнением движения свободной границы:

=-A- V(nT Ът (nT )). (18)

т

Для питательного вещества полагалось

c(x,0,t) = и - (x,t) = 0 при x 6dGx=0 uдGсв. (19)

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

Уравнение (13) является по форме квазилинейным уравнением теплопроводности с источником (стоком) dnT/dt - div (k(nT )VnT ) = Q(nT, с). Соответствующий нелинейный

коэффициент диффузии будет иметь вид:

k(nT) = nT Ът (nT) - n-.-- (T (nT)).

Очевидно, что при выборе Ът (nT) = nT - n0, допустимыми значениями nT, при которых k(nT) > 0, будут nT > n0/2. При этом, если обеспечить nT > n0/2 + s/n0, то для k(nT) верна оценка k(nT) > е >0.

Предположим,что имеется избыток питательных веществ, тогда функция источника (стока) может быть представлена как Q(nT) = у nT (nmax - nT).

Достроим область G в второй, третьей и четвёртых координатных четвертях плоскости {x, y} симметрично относительно осей координат:

G = ((x,y): (Ix\,\y)6 G}.

Доопределим задачу (13), (14) для новой области G . Зададим симметричное начальное условие nT(x, 0) = n0, для x 6 G . Будем использовать те же условия (16), (18) на новой свободной границе dG .

Поскольку и Q(nT), и k(nT) являются гладкими при nT > 0, и Q(0) = 0, то при выполнении условия k(nT) > е = const > 0, существует локальное классическое решение рассматриваемой краевой задачи, а поскольку nT(x, 0) > 0 и n0 Ф 0, то nT(x, t) > 0 для любых t > 0 [16].

Однако поскольку в области пониженной концентрации питательных веществ возможны отрицательные значения Q(nT, с), то полная задача с учётом распределения питательных веществ (13), (14) требует дополнительного исследования. Рассуждения о



существовании классического решения верны, лишь если выполнено условие nT > n0/2 + s/n0. Данному условию удовлетворяют найденные численные решения.

5. Численный метод

Проводилось численное

моделирование решения для случая

£т (nT) = nT - n0. (20)

Масштабы входящих в задачу величин выбирались таким образом,

чтобы D = 1, nmax = 1, Р = 1, cb = 1.

В соответствие области G ставилась регулярная неравномерная сетка Gij = {gy} = {x, уу), i = 0,...,N Рис. 1. Трёхточечный шаблон {A,B,C} на j = 0,...,Nj так, что граничные узлы сетки регулярной неравномерной сетке Gtj.

лежали на границе области G. Для каждого внутреннего узла gi,j- соседними узлами считались узлы g+hj , g-hj, gtJ+1 и gyW.

Для произвольной сеточной функции fi j, заданной на Gi j , производные первого порядка по пространственным переменным аппроксимировались на трёхточечном шаблоне {A,B,C} (Рис. 1), с использованием следующих соотношений:

f (B) - f (A) = dx(A)(xB - xA) +df-(A)(yB -yA) + O((xB - xA)2 + (yB -yA)2),

f (C) - f (A) = £(A)(xc - xA) +dy (A)(yc - yA) + O((xc - xA)2 + (yc - yA)2).

Решение этой системы относительно df/dx(A) и df/dy(A) даёт:

~yA)fc, (21)

xB) fc, (22)

где dA = (xc - xA)(yB-yA) - (xb - xA)(yC - yA), а xp, yp, p = A, B, C - координаты узлов и значения сеточной функции для точек A, B и C соответственно. Такая аппроксимация имеет первый порядок точности и допустима для произвольного расположения узлов A, B, C, не лежащих на одной прямой. Оператор численного дифференцирования на шаблоне {gy, gi-1,j, gij-A обозначим D-, а на шаблоне {gy, gt+hj, gij+г} - D+.

Использованный при построении аппроксимации (21-22) метод неопределённых коэффициентов широко известен и часто используется на практике [17].

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

Для численного решения квазилинейного уравнения параболического типа диффузионное слагаемое приводилось к следующему консервативному виду [18]:

div (nT V(nT ST (nT ) )) = div (V K(nT ) ),


df (A) = ~г(Ус -Ув) fA + -tCVa - Ус) fB +г(Ув

xdAdAdA

ydAdAdA



где

K(nT) = j uLT (u) + u 2Y!T (u) du. 0

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

div(nTV(nTЪт (nT))) = D- (D+ K(nT)) + O(hi j),

аналогично записывается аппроксимация оператора Лапласа в уравнении для концентрации питательного вещества с.

o.ai

0.79 (°7л ш747

0.75


0 о. ?д.2оStadtfaeМол6.9

0.9 В.85

0.8 0.75

0.7 0.65

0.6 0.55

0.5 0.45

0 0.5 1

max(x) -

<(V) -

----

..............S

....................>

...............>

1.5 t

2 2.5 3

Рис. 2. Рост опухоли вдоль сосуда. Рис. 3. Рост опухоли вдоль сосуда. Положение границы опухоли в Изменение линейных размеров

последовательные моменты времени опухоли по оси x (сплошная линия)

и по оси y (пунктирная линия).

Рис. 4. Профиль опухолевого фронта вдоль оси y = 0 в момент времени t = 2,9. На внешней границе dnT/dx < 0 , граница смещается вправо.

0.81 0.8 0.79 0.78 0.77 0.76 0.75 0.74

0 0.1 0.2 0.3 0.4 0.5 0.6 У

Рис. 5. Профиль опухолевого фронта вдоль оси x = 0 в момент времени t = 2,9. На внешней границе dnT/dx 0 , граница почти неподвижна.


В ходе численных расчётов было установлено, что удлинение опухоли опухоли вдоль оси x (сосуда), происходит значительно быстрее, чем изменение её радиуса - максимального линейного размера вдоль оси y (Рис. 2, 3).



0.50

-0.25

0.00

и г / / / / / S /

Рис. 6. Поле скоростей клеточной фазы растущего тяжа (t = 2,9). В области постоянного радиуса (слева) преобладает радиальная компонента.

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

Заключение

Как и модели [4, 5], рассмотренная модель описывает анизотропный опухолевый рост в среде с неоднородным распределением питательных веществ.

Ключевой особенностью модели (13)-(20) является отсутствие предположений о собственной подвижности клеток (хаотической или направленной), вместо это делаются пред положения о характере взаимодействия между клетками опухоли и механических свойствах ткани. Следствием непрерывности напряжений являются конечная скорость движения клеток

vT =-- V(nT 2L, (nT))

и сохранение компактной формы опухоли.

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

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

Работа выполнена при финансовой поддержке фонда Марии Кюри в рамках проекта Modelling Mathematical Methods and Computer Simulation of Tumour Growth and Therapy ,

Профили движущихся фронтов опухоли (в разрезе вдоль осей y = 0 и x = 0) представлены на рис. 4, 5.



MCRTN-CT-2004-503661. Авторы выражают благодарность д.ф.-м.н., проф. А. И. Лобанову за помощь в подготовке статьи.

Литература

1. Scientific yearbook 2004/2005: Teoh. rep.: Cancer Research UK, 2005.

2. Cancer Modelling and Simulation / Ed. by L. Preziosi.- CRC Press, 2003.- 456 p.

3. Araujo R. P., McElwain D. L. S. A history of the study of solid tumour growth: The contribution of mathematical modelling Bulletin of Mathematical Biology.- 2004.- no. 66.- Pp. 10391091.

4. Астанин С. А., Лобанов А. И. Трехмерная модель роста неваскуляризованной опухоли в ткани Математика. Компьютер. Образование . XII Международная конференция: Сб. научных трудов.- М. - Ижевск: Регулярная и хаотическая динамика, 2005.- Т. 2.- С.

759-769.

5. Колобов А. А., Полежаев А. В. Направленный рост и инвазия опухоли в отсутствии хемотактической подвижности ее клеток Математика. Компьютер. Образование . XII Международная конференция: Сб. научных трудов.- М. - Ижевск: Регулярная и хаотическая динамика, 2005.- Т. 3.- С. 979-989.

6. Влияние пространственной гетерогенности среды на рост и инвазию опухоли. Анализ методами математического моделирования / С. А. Астанин, А. В. Колобов, А. И. Лобанов и др. Медицина в зеркале информатики.- 2006.- в печати.

7. Moore J. V., Hopkins H. A., Looney W. B. Tumour-cord parameters in two rat hepatomas that differ in their radiobiological oxygenation status Radiat. Envir. Biophys.- 1984.- no. 23.- Pp.

213-222.

8. Migration and internalization of cells and polystyrene microspheres in tumour cell spheroids / M. J. Dorie, R. F. Kallman, D. Antwerp, Y. R. Huang Exp. Cell Res.- 1982.- no. 141.- Pp. 201-209.

9. Bertuzzi A., Gandolfi A. Cell kinetics in a tumour cord J. Theor. Biol.- 2000.- no. 204.-

Pp. 587-599.

10. Cell kinetics in tumour cords studied by a model with variable cell cycle length / A. Bertuzzi, A. Fasano, A. Gandolfi, D. Marangi Math. Biosci.- 2002.- no. 177-178.- Pp. 103-125.

11. Bertuzzi A., Fasano A., Gandolfi A. A free boundary problem with unilateral constraints describing the evolution of a tumor cord under the influence of cell killing agents SIAM J. Math. Anal.- 2004.- Vol. 36, no. 3.- Pp. 882-915.

12. Byrne H., Preziosi L. Modelling solid tumour growth using the theory of mixtures Mathematical Medicine and Biology.- 2003.- no. 20.- Pp. 341-366.

13. Ambrosi D., Preziosi L. On the closure of mass balance models for tumor growth Mathematical Models and Methods in Applied Science.- 2002.- Vol. 12, no. 5.- Pp. 737-754.

14. Полубаринова-Кочина П. Я. Теория движения грунтовых вод.- М.: Наука, 1977. - 664 с.

15. Механика насыщенных пористых сред / В. Н. Николаевский, К. С. Басниев, А. Т. Горбунов, Г. А. Зотов.- М.: Недра, 1970.- 339 с.

16. Режимы с обострением в задачах для квазилинейных параболических уравнений /

А. А. Самарский, В. А. Галактионов, С. П. Курдюмов, А. П. Михайлов.- М.: Наука. Гл. ред. физ.-мат. лит., 1987.- 480 с.

17. Магомедов К. М., Холодов А. С. Сеточно-характеристические численные методы.- М.: Наука, 1988.- 287 с.



18. Evje S., Karlsen K. H. Monotone difference approximations of BV solutions to degenerate convection-diffusion equations SIAM J. Numer. Anal.- 2000.- Vol. 37, no. 6.- Pp. 18381860.



© 2017 РубинГудс.
Копирование запрещено.