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

УДК 519.63

DOI: 10.14529/ mmp170208

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

Л. С. Исаченко, А.И. Лобанов

Московский физико-технический институт (государственный университет), г. Долгопрудный

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

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

Постановка задачи

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

В данной работе предлагается модель с учетом диффузионной и хемотаксической подвижности опухолевых клеток, а также гетерогенного окружения, основанная на подходе 111. Супцественнъхми недостатками [1] являются возможность роста суммарной концентрации клеток до любого уровня, а также бесконечная скорость распро-отрешения опухолевых подвижных клеток. Одним из путей решения этих проблем является постановка вариационной задачи, аналогично [2], с учетом ограничений типа неравенств на полный поток и суммарную численность клеток. Другой путь — использование феноменологического ограничения потока подвижных клеток. Такой подход использовался в [3] для моделирования релаксационного теплопереноса в плазме. Здесь он применяется для биологической задачи. С учетом ограничения потока

система уравнений в одномерном случае имеет вид:

+ Р (Б) т - Р2 (Б) п - ¿пщ (1)

дп д дЬ дх

с \ ( дп дБ 1--) I -О— + ап—

стах ) \ дх дх

дт „ / с

дтБ (1 - - Рг (Б) т + Р2 (Б) п, (2)

дЬ \ к

с

к

— = Ds- 11Б1 - 1тБт - 7иБп, (4)

д] / с \

— = в]Б (1 - к) - Н] (п + т) , (3)

пт

]

Б — концентрация питательного субстрата, с = п + т + ] — суммарная плотность клеток, О — коэффициент диффузии подвижных опухолевых клеток, а — коэффициент хемотаксиса по отношению к субстрату, q — скорость роста неподвижных опухолевых клеток, Рг (Б) = К 1-гк(€(8-£!сг)) ^ р2 (б) = К — функции итенсивности перехода из неподвижного состояния в подвижное и обратно.

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

ЗЯПИСЫВавТСЯ В ВИД6!

™ = 1-- - О——+ ап-—

V стах) V дх дх) '

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

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

стах

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

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

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

этой границы. Ее положение можно описать функцией £ (Ь). Найдем скорость движения гр&ницы. Предположим без ограничения общности, что она движется слева направо. Пусть в момент времени Ь граница находится в £ (Ь), за некоторое малое время ДЬ она переместится в точку £ (Ь + ДЬ) и пройдет расстояние VДЬ, где V — скорость подвижной границы. Так как время ДЬ мало, скорость V можно считать постоянной.

V

_L

£(1+А1) х

Рис. 1. К оценке фазовой скорости границы Просуммировав (1) - (3), получим уравнение баланса для суммарной плотности

клеток

dc dt

d_ dx

(

/-i \ I т^ди dS (1 - c) -D— + an— 1 dx dx

+ g,

(5)

где g = qmS (1 — c/k) + [US (1 — c/k) — Hl (и + m) — dnn. Интегрируем уравнение баланса (5):

i+Ai i+Ai

[ —dt = ( dW

I dt I dx

i+Ai

i+Ai

?(i+Ai)

1 f dW , .

dt + g dt =—— — dx + g dt.

dx

i m i

ДЬ

раздела фаз, можно записать

с (Ь + ДЬ) - с (Ь) = - V (V (£ (Ь + ДЬ)) - Ж(£ (Ь)}) + (д (Ь + ДЬ) - д (Ь))ДЬ.

Считая, что при малых значениях ДЬ значение потока (Ь)) соответствует потоку слева от границы, аЖ(£ (Ь + ДЬ)) — справа, получим выражение для скорости гряницы

V--

[W ]

[g] At — [c]

(6)

где [Ж] — скачок полного потока, [с] — скачок суммарной концентрации клеток на границе, [д] — скачок функции д.

С

Также, исходя из балансовых соотношении, получим условие на разрыве для потока

Г Ш 1

0.

1с

(7)

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

пх (0,Ь) = 0, пх (Ь,Ь) = 0, тх (0,Ь) = 0, тх (Ь,Ь) = 0, ]х (0,Ь) = 0, ]х (Ь,Ь) = 0,

Бх (0,Ь) = 0, Б (Ь,Ь) = 1.

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

1. Разностная схема

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

Система уравнений приводится к безразмерному виду. В области решения вводится равномерная по пространству расчетная сетка, выбор шага по времени при послойном переходе осуществляется автоматически. Используется двухслойная неявная разностная схема. Для значений на верхнем слое используются обозначения с Для построения консервативной разностной схемы необходимо пользоваться дивергентной формой записи уравнений. С этой целью преобразуем выражение -дх [(1 - с) {-О^дп + апдхУ\ , воспользовавшись аппроксимацией Шарфеттера [8]. Получим

д_ дх

(

дп дп (1 - с)\ -о— + ап— дх дх

д_

дх

/-, ч ^ аЗ д ( (1 - с) ^ -

Разностное уравнение для подвижных опухолевых клеток имеет вид

п — п

О

т

2к

к

+

^ _ а с _а с

п+ге о- пе о° л

1 -;--К_

2

_ а с л. _ а с

пе о° - п_1е о°-1

к

(1 - с+1 - с+1)

к

+1 (п,т,С

п, т, Б) - йп, (8)

е - ^ е о'

f (п,т, Б^

Kqm

г-т^-всг)

2

где к+1 = _ , _ ..

^ * ' а с а с

2 е + е о Д+1

Kqn-(-2г-——.Значения велич ин на п - 1и п +1 шагах по пространству обозначаются индексами —1 и +1.

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

= (Б+г - 2 Б + ¡3—1^ - 7эБп - 72Бт - 71 Бке,

2

/ -l

= в IS (1 - - ) - HI (n + m),

m — m

(1 - i)

(1 - i)-'('

ФЯ (1 - - ) - / (п,т,8) .

Модифицируем разностные уравнения вблизи подвижной границы. Рассмотрим случай, когда граница движется слева направо. Для аппроксимации функции £ (¿) используется кусочно-линейная экстраполяция

I = £ + Ут,

где У — скорость границы на нижнем слое.

Так как в £ разрыв, запишем производные по пространству с использованием дополнительной точки £ . Обозначим значения величин справа от разрыва индексом £+, слева — индексом £_. Тогда слева от границы разностное уравнение для подвижных опухолевых клеток имеет вид

nn

т

2D

(1 - с + 1 - с—)1

- neD (SZ--S)'

xn + h £ - xn ( 1 + eDD(s«--s) 2D (1 - с +1 - с -i)ineD(S-1-S) - n-t'

£ - xn + h

Справа от границы

h

1 + e d

(S-!-S)

+ f (n, m, sj - dnn.

nn

2D

(1 - с + 1 - с+i) ( n+1

- neD (S+1-S)'

xn - £ + h

2D

h

1 + e D -S)

(1 - с +1 - с^+) / ne

DD &+-S) _

щ+

xn - £ + h xn - £

1 + e DD (S«+-S)

+ f (n, m, sj - dnn.

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

nn

т

2D (1 - с +1 - с+i) ^ e-D^^

£ - xn + ^ 1 + eDD(s+I-S)

2D (1 - с +1 - с+1) Адт

h + D (с - с+i)

£ - Xn + h

(

1 + e Dl

(S+-S) %^ (h + V (с - с+i))

+

2D (1 - с +1 - с -i)ineDD^ - n-1

£ - Xn + h (1 + e DD^

h

+ f (n,m, s) - dnn,

где Ag = в IS (1 - §) + qmS (1 - §) - dnn - вf+i^+i (1 - - qm+iS+i (1 - + dnn+i.

Справа:

* ОП ,Л - , 1 - N (п+1-пеЪе-Ъ^^

с -п 2Б (1 - с+1 - с+г) ^ п+1 Се )е

Т С - Хп + Ь ^ + еъ (5+1-5) Ь + V (С - С+1)

2Б (1 - С+1 - С+1) Адт

с - Хп + Ь (1 + еЪ (5+1-5)+ V (С - С+1))

^1 + е ъ (¿+1-5)^^

2Б ( - С +1 - С-)) . (СеЪ(5-1-5) -С-Л + . (й7Н§]

С - Хп + Ь (1 + еъ (5-1-5А ^ ь )+ ПС'т

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

2. Результаты расчетов

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

Рис. 2. Плотность клеток и концентрация субстрата при отсуствии подвижных границ. Значения параметров Б = 0, 026, а = 0, 2, q = 100

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

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

Рис. 3. Плотность опухолевых клеток и распределение субстрата при наличии подвижных границ (отсуствует вырождение в гиперболический тип уравнения). Значения параметров Б = 0,026, а = 2, д = 100

Рис. 4. Плотность опухолевых клеток и распределение при наличии подвижных границ (вырождение в гиперболический тип уравнения). Значения параметров Б = 0,026, а = 2, д =150

Заключение

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

Литература

1. Астанин, С.А. Влияние пространственной гетерогенности среды на рост и инвазию опухоли. Анализ методами математического моделирования / С.А. Астанин, А.В. Колобов,

A.И. Лобанов, Т.П. Пименова, А.А. Полежаев, Г.И. Соляник // Медицина в зеркале информатики. - М.: Наука, 2008. - С. 188-223.

2. Албу, А.Ф. О процессе плавления с ограничением на скорость остывания / А.Ф. Албу,

B.И. Зубов // Математическое моделирование. - 2002. - Т. 14, № 5. - С. 119-123.

3. Волосевич, П.П. Динамика и нагрев плазмы с учетом релаксации теплового потока / П.П. Волосевич, Н.В. Змитренко, Е.И. Леванов, Е.В. Северина // Математическое моделирование. - 2008. - Т. 20, № 4. - С. 57-68.

4. Freeman, А.Е. In Vivo-Like Growth of Human Tumors in Vitro / A.E. Freeman, R.M. Hoffman // Proceedings of the National Academy of Sciences of the United States of America. - 1986. - V. 83, № 8. - P. 2694-2698.

5. Hsiao, A.Y. Microfluidic System for Formation of PC-3 Prostate Cancer Co-Culture Spheroids / A.Y. Hsiao, Y. Torisawaa, Y.C. Tunga, S. Sudb, R.S. Taichmanc, K.J. Pientab, S. Takayamaet // Biomaterials. - 2009. - V. 30. - P. 3020-3027.

6. Федоренко, Р.П. Разностная схема для задачи Стефана / Р.П. Федоренко // Журнал вычислительной математики и математической физики. - 1975. - Т. 15, № 5. - С. 1339-1344.

7. Низьев, ВТ. Численное моделирование плавления двухкомпонентных порошков при лазерном спекании / В. Г. Низьев, А. В. Колдоба, Ф. X. Мирзаде, В. Я. Панченко, Ю. А. Повещенко, М. В. Попов// Математическое моделирование. - 2011. - Т. 23, № 4. -

C. 90-102.

8. Scharfetter, D.L. Large-Signal Analysis of Silicon Read Diode Oscillator / D.L. Scharfetter, H.K. Gummel // IEEE Transactions on Electron Devices. - 1969. - V. 16, № 1. - P. 64-77.

Лариса Сергеевна Исаченко, аспирант, кафедра «Информатика и вычислительная математика:», Московский физико-технический институт (государственный университет) (Московская область, г. Долгопрудный, Российская Федерация), isachenko.larisa@gmail.com.

Алексей Иванович Лобанов, доктор физико-математических наук, профессор, кафедра «Информатика и вычислительная математика», Московский физико-технический институт (государственный университет) (Московская область, г. Долгопрудный, Российская Федерация), alexey@crec.mipt.ru.

Поступила в редакцию 30 июня 2016 г.

MSC 65N06 DOI: 10.14529/mmp170208

A DIFFERENCE SCHEME FOR SOLVING THE EQUATIONS OF TUMOR GROWTH SUBJECT TO THE RESTRICTED FLOW OF MOTILE CELLS L.S. Isachenko, L.I. Lobanov

Moscow Institute of Physics and Technology (State University), Dolgoprudny, Moscow Region, Russian Federation

E-mail: isachenko.larisa@gmail.com, alexey@crec.mipt.ru

The article investigates one-dimensional mathematical model of tumor growth represented by a system of quasi-linear parabolic equations. We assume certain restrictions

on the full flow of the motile tumor cells, leading to the possible degeneration of the system into a hyperbolic type and emergence of discontinuous (weak) solution. To find weak solution we consider tumor growth as the emergence of a new phase. Thus we have a generalized (nonlinear) Stefan problem. The authors propose and implement a difference scheme with the explicit statement of the phase-change moving boundary to solve the problem. It is shown that this approach allows to describe different regimes of tumor growth.

Keywords: difference scheme; substrate taxis; the problem with movable boundary; break allocation.

References

1. Astanin S.A., Kolobov A.V., Lobanov A.I., Pimenova T.P., Polezhaev A.A., Solynik G.I. [An

Effect of Spatial Heterogeneity of the Environment on the Tumor Growth and Invasion. Analysis Using Methods of Mathematical Modelling]. Meditsina v zerkale informatiki [Medicine in the Mirror of Informatics]. Moscow, Nauka, 2008, pp. 188-223.

2. Albu A.F., Zubov V.I. [On a Melting Process with Restriction on a Cooling Velocity]. Mathematical Models and Computer Simulations, 2002, vol. 14, no. 8, pp. 119-123. (in Russian)

3. Volosevich P.P., Zmitrenko N.V., Levanov E.I., Severina E.V. [Dynamic and Heating of Plasma Subject to Heat Flux Relaxation]. Mathematical Models and Computer Simulations, 2008, vol. 20, no. 4, pp. 57-68. (in Russian)

4. Hsiao A.Y., Torisawaa Y., Tunga Y.-C., Sudb S., Taichmanc R.S., Pientab K.J., Takayamaet S. Microfluidic System for Formation of PC-3 Prostate Cancer Co-Culture Spheroids. Biomaterials, 2009, vol. 30, pp. 3020-3027. DOI: 10.1016/j.biomaterials.2009.02.047

5. Freeman A.E., Hoffman R.M. In Vivo-Like Growth of Human Tumors in Vitro. Proceedings of the National Academy of Sciences of the United States of America, 1986, vol. 83, pp. 2694-2698. DOI: 10.1073/pnas.83.8.2694

6. Fedorenko R.P. A Difference Scheme for Stefan's Problem. USSR Computational Mathematics and Mathematical Physics, 1975, vol. 15, no. 5, pp. 246-251. DOI: 10.1016/0041-5553(75)90122-6

7. Niziev V.G., Koldoba A.V., Mirzade F.Kh., Panchenko V.Ya., Poveschenko Yu.A., Popov M.V. [Numerical Modelling of Laser Sintering of Two-Component Powder Mixtures]. Mathematical Models and Computer Simulations, 2011, vol. 23, no. 4, pp. 90-102. (in Russian)

8. Scharfetter D.L., Gummel H.K. Large-Signal Analysis of Silicon Read Diode Oscillator. IEEE Transactions on Electron Devices, 1969, vol. 16, no. 1, pp. 64-77. DOI: 10.1109/T-ED.1969.16566

Received 30 June, 2016