Известия РАН. Механика жидкости и газа, 2020, № 6, стр. 3-16

КОНТИНУАЛЬНОЕ МОДЕЛИРОВАНИЕ БИОЛОГИЧЕСКОЙ СРЕДЫ, СОСТАВЛЕННОЙ АКТИВНО ВЗАИМОДЕЙСТВУЮЩИМИ КЛЕТКАМИ ДВУХ РАЗНЫХ ТИПОВ

С. А. Логвенков a*, А. А. Штейн b**

a Национальный исследовательский университет “Высшая школа экономики”
Москва, Россия

b МГУ им. М.В. Ломоносова, Научно-исследовательский институт механики
Москва, Россия

* E-mail: logv@bk.ru
** E-mail: stein.msu@bk.ru

Поступила в редакцию 22.05.2020
После доработки 18.06.2020
Принята к публикации 18.06.2020

Полный текст (PDF)

Аннотация

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

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

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

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

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

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

Бо́льшая часть работ, посвященных континуальному моделированию активных клеточных систем, использует уравнения диффузионного типа для объемных плотностей клеток и матрикса (внеклеточных структур ткани) [4, 611]. Потоки плотностей определяются соотношениями, учитывающими, помимо эффективной диффузии, соответствующей случайным клеточным блужданиям, также активную нелокальную составляющую, которая интерпретируется как проявление межклеточной адгезии. Модели этой группы применялись, в том числе, и при моделировании сортировки клеток в гетерогенных клеточных смесях [4, 6, 10]. В качестве граничных условий при решении модельных задач использовались условия периодичности или отсутствия потока через неподвижную границу. Недостатком такого рода моделей оказывается игнорирование связи потоков с параметрами, характеризующими развиваемые клетками усилия и напряженное состояние среды.

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

Нами разработан общий подход к моделированию процессов в механически активных клеточных средах, наиболее полно изложенный в [16] и основанный на последовательном применении методов механики многофазных сред с учетом как локальной, так и нелокальной зависимости активных напряжений от характеризующих среду параметров. Этот подход был применен к анализу деформационных процессов в тканях различной структуры. В [17, 18] изучались активные переупаковки в двумерном слое, составленном плотно прилегающими одна к другой клетками (эмбриональный эпителий), а в [19, 20] – активные деформации в трехмерной среде, состоящей из клеток одного сорта и межклеточной жидкости. В настоящей работе в рамках такого подхода разработана модель механически активной среды, образованной клетками двух различных типов в присутствии жидкой фазы. Представление о “поверхностном натяжении” в клеточной среде, на котором базируются многие теоретические исследования, при построении модели непосредственно не используется. Объясняемые с его помощью эффекты удается описать с помощью физически более обоснованного механизма активных межклеточных взаимодействий. Полученная модель применена к решению задачи о сортировке клеток в бесконечном плоском слое.

1. ОБЩАЯ МНОГОФАЗНАЯ МОДЕЛЬ

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

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

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

Далее среда предполагается состоящей из двух клеточных (обозначаются индексами 1 и 2) и жидкой (индекс w) фаз, а также трех дополнительных, посредством которых учитываются активные взаимодействия как между клетками внутри одного типа (индексы 1 и 2), так и между клетками, относящимися к разным клеточным фазам (индекс 12). Далее будем говорить о клеточных фазах 1 и 2, жидкой фазе w и дополнительных фазах 1, 2 и 12. Истинные плотности фаз, т.е. плотности жидкости и клеток, будем считать постоянными и равными между собой. Объемные концентрации трех основных фаз, которые равны, таким образом, их массовым концентрациям и пропорциональны распределенным плотностям, обозначим ${{\phi }_{1}}$, ${{\phi }_{2}}$ и ${{\phi }_{w}}$ соответственно. Далее, говоря о плотностях фаз, всегда будем иметь в виду распределенную плотность. Объемом дополнительных фаз будем пренебрегать по сравнению с объемами клеточных и жидкой фаз, что соответствует выполнению соотношения ${{\phi }_{1}} + {{\phi }_{2}} + {{\phi }_{w}} = 1$.

Тензор полных напряжений в среде T представляется суммой напряжения во внеклеточной жидкости, характеризуемого гидростатическим давлением p, напряжений в двух клеточных фазах ${{\sigma }^{{(1)}}}$ и ${{\sigma }^{{(2)}}}$, а также напряжений в дополнительных фазах ${{\tau }^{{(1)}}}$, ${{\tau }^{{(2)}}}$ и ${{\tau }^{{(12)}}}$. Тензоры ${{\tau }^{{(1)}}}$ и ${{\tau }^{{(2)}}}$ характеризуют напряжения в дополнительных фазах, обеспечивающих активные межклеточные взаимодействия между клетками только первой и только второй клеточных фаз соответственно, а тензор ${{\tau }^{{(12)}}}$ описывает взаимодействие клеток, принадлежащих к разным фазам:

(1.1)
${{T}^{{ij}}} = - {{\phi }_{w}}p \times {{g}^{{ij}}} + {{\phi }_{1}}{{\sigma }^{{(1)ij}}} + {{\phi }_{2}}{{\sigma }^{{(2)ij}}} + {{\tau }^{{(1)ij}}} + {{\tau }^{{(2)ij}}} + {{\tau }^{{(12)ij}}}$

Определяющие соотношения для активных напряжений в дополнительных фазах учитывают хаотическую активность клеток, характеризуемую соответствующими давлениями, которые проявляются в отталкивании клеток. Будем также предполагать, что в активных напряжениях присутствует стягивающая составляющая, определяемая неоднородностью распределения клеток по сферической поверхности $S$ некоторого фиксированного для данной среды радиуса Rs (радиуса дальнодействия) с центром в рассматриваемой точке. Этот радиус характеризует эффективное расстояние, на котором происходит силовое взаимодействие клеток и имеет порядок одного или нескольких клеточных размеров [20].

В локальной декартовой системе координат с центром в рассматриваемой точке примем следующие зависимости для тензоров активных напряжений от объемных плотностей клеток

(1.2)
${{\tau }^{{(1)}}}(x) = - {{\Pi }_{1}}{\rm I} + {{m}_{1}}\frac{3}{{4\pi R_{s}^{2}}} \cdot {{\phi }_{1}}(x)\int\limits_S {{{\phi }_{1}}} (x + r)N \otimes Nds$
(1.3)
${{\tau }^{{(2)}}}(x) = - {{\Pi }_{2}}{\rm I} + {{m}_{2}}\frac{3}{{4\pi R_{s}^{2}}} \cdot {{\phi }_{2}}(x)\int\limits_S {{{\phi }_{2}}} (x + r)N \otimes Nds$
(1.4)
${{\tau }^{{(12)}}}(x) = - {{\Pi }_{{12}}}{\rm I} + {{\tau }^{{[12]}}}(x) + {{\tau }^{{[21]}}}(x)$
(1.5)
$\begin{gathered} {{\tau }^{{[12]}}}(x) = \frac{3}{{4\pi R_{s}^{2}}} \cdot {{m}_{{12}}}{{\phi }_{1}}(x)\int\limits_S {{{\phi }_{2}}(x + r)} N \otimes Nds \\ {{\tau }^{{[21]}}}(x) = \frac{3}{{4\pi R_{s}^{2}}} \cdot {{m}_{{21}}}{{\phi }_{2}}(x)\int\limits_S {{{\phi }_{1}}(x + r)} N \otimes Nds \\ \end{gathered} $

Здесь $r$ – радиус-вектор точки на поверхности сферы, I – единичный тензор, N – единичный вектор нормали к поверхности сферы. Для точек, находящихся на расстоянии, меньшем Rs от поверхности, ограничивающей рассматриваемый объем среды, интегрирование ведется только той части поверхности сферы, которая лежит внутри этого объема.

Первые слагаемые в (1.2)–(1.4) описывают нелинейное активное сопротивление клеток их избыточному сближению, обусловленное хаотическим отталкиванием клеток. Будем полагать положительно определенные скаляры ${{\Pi }_{1}}$, ${{\Pi }_{2}}$ и ${{\Pi }_{{12}}}$ функциями объемных плотностей клеток ${{\phi }_{1}}$ и ${{\phi }_{2}}$, определяющих степень такого сближения:

(1.6)
${{\Pi }_{1}} = {{\Pi }_{1}}({{\phi }_{1}},{{\phi }_{2}}),\quad {{\Pi }_{2}} = {{\Pi }_{2}}({{\phi }_{1}},{{\phi }_{2}}),\quad {{\Pi }_{{12}}} = {{\Pi }_{{12}}}({{\phi }_{1}},{{\phi }_{2}})$
причем эти функции должны возрастать с увеличением объемных плотностей клеток и стремиться к бесконечности при ${{\phi }_{1}} + {{\phi }_{2}} \to 1$, т.е. при приближении к плотной упаковке. Кроме того, функция ${{\Pi }_{1}}({{\phi }_{1}},{{\phi }_{2}})$ стремится к нулю при ${{\phi }_{1}} \to 0$, функция ${{\Pi }_{2}}({{\phi }_{1}},{{\phi }_{2}})$ – при ${{\phi }_{2}} \to 0$, а функция ${{\Pi }_{{12}}}({{\phi }_{1}},{{\phi }_{2}})$ – при стремлении к нулю любого из ее аргументов.

Последующие слагаемые в (1.2)–(1.4) имеют смысл напряжений, создаваемых в результате направленных активных межклеточных взаимодействий. Вид этих слагаемых в (1.2) и (1.3) основан на предположении, что активное напряжение в дополнительной фазе 1 или 2 пропорционально концентрации фазы в данной точке и ее концентрации на расстоянии ${{R}_{s}}$ от этой точки с учетом неоднородности распределения этой концентрации по сфере $S$. В уравнении (1.4) соответствующее слагаемое построено, как видно из (1.5), аналогично, но с учетом присутствия двух различных механизмов, описываемых тензорами ${{\tau }^{{[12]}}}$ и ${{\tau }^{{[21]}}}$. Первый из них соответствует активной роли клеток фазы 1 (первое соотношение (1.5)), а второй – клеток фазы 2 (второе соотношение (1.5)). Если взаимодействие клеток осуществляется посредством псевдоподий, активную роль исполняют клетки, которые выбрасывают псевдоподии и с их помощью распластываются и ползут по другим клеткам, а пассивную – клетки, играющие роль субстрата для распластывания активных [20]. Коэффициенты ${{m}_{1}}$, ${{m}_{2}}$, ${{m}_{{12}}}$ и ${{m}_{{21}}}$ положительны в соответствии с тем, что направленные активные напряжения (в отличие от хаотических) стягивают, а не расталкивают клетки.

Для среды, образованной одной клеточной фазой, в [20] аналогичное (1.2) уравнение для активных напряжений сводилось к дифференциальному посредством разложения подынтегральной функции в ряд Тейлора до членов второго порядка. В этом случае анизотропия активных напряжений определяется вторыми производными объемной плотности клеток по координатам. В настоящей работе в расчетах используется непосредственно исходная интегральная форма.

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

(1.7)
$\begin{gathered} {\text{div}}({{\phi }_{1}}{{\sigma }^{{(1)}}}) + {{F}_{{12}}} + {{F}_{{1w}}} + {{F}_{{1,\tau 1}}} + {{F}_{{1,\tau 12}}} = 0 \\ {\text{div}}({{\phi }_{2}}{{\sigma }^{{(2)}}}) - {{F}_{{12}}} + {{F}_{{2w}}} + {{F}_{{2,\tau 2}}} + {{F}_{{2,\tau 12}}} = 0 \\ - \nabla ({{\phi }_{w}}p) - {{F}_{{1w}}} - {{F}_{{2w}}} = 0 \\ {\text{div}}\,{{\tau }^{{(1)}}} - {{F}_{{1,\tau 1}}} = 0 \\ {\text{div}}\,{{\tau }^{{(2)}}} - {{F}_{{2,\tau 2}}} = 0 \\ {\text{div}}\,{{\tau }^{{(12)}}} - {{F}_{{1,\tau 12}}} - {{F}_{{2,\tau 12}}} = 0 \\ \end{gathered} $

Сумма всех уравнений дает с учетом (1.1) уравнение равновесия для среды в целом ${\text{div}}\,T = 0$. Здесь F12 – сила взаимодействия между первой и второй клеточными фазами, ${{F}_{{kw}}}$($k = 1,2$) – сила взаимодействия между k-й клеточной фазой и жидкостью, ${{F}_{{k,\tau k}}}$ – сила взаимодействия между k-й клеточной фазой и связанной только с ней дополнительной фазой, ${{F}_{{k,\tau 12}}}$ – сила взаимодействия между k-й клеточной фазой и “перекрестной” дополнительной фазой.

В общем случае для биологических тканей, образованных перемещающимися одна относительно другой клетками, тензоры скоростей деформации клеточных фаз можно представить в виде двух составляющих: связанной с деформацией самих клеток и определяемой их переупаковкой [17]. В настоящей работе примем допущение о малости слагаемого, описывающего скорость деформации этих фаз за счет деформации клеток в сравнении со скоростью неупругого деформирования, связанного с переупаковками. Будем полагать, что течение клеточных фаз определяется активными напряжениями, причем течением каждой фазы непосредственно управляют активные напряжения как в связанной с ней дополнительной фазе, так и в перекрестной дополнительной фазе 12. Кроме того, на это течение оказывают влияние все напряжения, присутствующие в основных фазах w, 1 и 2: давление жидкости p и клеточные напряжения ${{\sigma }^{{(1)}}}$ и ${{\sigma }^{{(2)}}}$.

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

(1.8)
$\begin{gathered} {{{\mathbf{e}}}^{{(1)}}} = - \frac{1}{{2{{\mu }_{{11}}}}}{{\tau }^{{(1)}}} - \frac{1}{{2{{\mu }_{{112}}}}}{{\tau }^{{[12]}}} - \frac{1}{{2{{\mu }_{{121}}}}}{{\tau }^{{[21]}}} + \frac{1}{{2{{\mu }_{{1w}}}}}p{\mathbf{I}} + \frac{1}{{2{{\mu }_{{1c1}}}}}{{\sigma }^{{(1)}}} - \frac{1}{{2{{\mu }_{{1c2}}}}}{{\sigma }^{{(2)}}} \\ {{{\mathbf{e}}}^{{(2)}}} = - \frac{1}{{2{{\mu }_{{22}}}}}{{\tau }^{{(2)}}} - \frac{1}{{2{{\mu }_{{212}}}}}{{\tau }^{{[12]}}} - \frac{1}{{2{{\mu }_{{221}}}}}{{\tau }^{{[21]}}} + \frac{1}{{2{{\mu }_{{2w}}}}}p{\mathbf{I}} + \frac{1}{{2{{\mu }_{{2c2}}}}}{{\sigma }^{{(2)}}} - \frac{1}{{2{{\mu }_{{2c1}}}}}{{\sigma }^{{(1)}}} \\ \end{gathered} $

Коэффициенты, входящие в зависимости (1.8), могут, вообще говоря, зависеть от фазовых концентраций ${{\phi }_{1}}$ и ${{\phi }_{2}}$ и положительны. Развиваемые в дополнительных фазах активные растягивающие напряжения стимулируют сближение клеток, т.е. уплотнение (сжатие) клеточной фазы, с которой они связаны. Давление в жидкости, растягивающие клеточные напряжения в рассматриваемой фазе и сжимающие клеточные напряжения в соседней фазе этому уплотнению препятствуют. Раздельное в общем случае присутствие в (1.8) двух составляющих “перекрестной” активной силы ${{\tau }^{{(12)}}}$ связано с тем, что эти составляющие соответствуют разному характеру взаимодействия между клетками, принадлежащими к разным фазам.

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

Поскольку в рассматриваемой модели другие деформации, кроме описываемых зависимостями (1.8), как указано выше, не учитываются, тензоры e(1) и e(2) играют роль полных скоростей деформаций фаз и подчиняются кинематическим соотношениям

(1.9)
$e_{{ij}}^{{(1)}} = \frac{1}{2}\left( {{{\nabla }_{i}}{{{v}}_{{1j}}} + {{\nabla }_{j}}{{{v}}_{{1i}}}} \right),\quad e_{{ij}}^{{(2)}} = \frac{1}{2}\left( {{{\nabla }_{i}}{{{v}}_{{2j}}} + {{\nabla }_{j}}{{{v}}_{{2i}}}} \right)$

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

(1.10)
$\frac{{\partial {{\phi }_{1}}}}{{\partial t}} + \nabla \cdot ({{\phi }_{1}}{{{\mathbf{v}}}_{1}}) = 0,\quad \frac{{\partial {{\phi }_{2}}}}{{\partial t}} + \nabla \cdot ({{\phi }_{2}}{{{\mathbf{v}}}_{2}}) = 0,\quad \frac{{\partial {{\phi }_{w}}}}{{\partial t}} + (\nabla \cdot {{\phi }_{w}}{\mathbf{w}}) = 0$

2. О ВЫРАЖЕНИЯХ ДЛЯ МЕЖФАЗНЫХ СИЛ

В соотношениях (1.7) присутствуют 7 векторов, представляющих силы, действующие между фазами. Для дополнительных фаз 1 и 2 нет необходимости постулировать какие-либо определяющие соотношения. Это связано с тем, что они “привязаны” к соответствующим клеточным фазам, т.е. для них не рассматриваются самостоятельные законы движения (и соответственно кинематические характеристики). Дополнительная фаза 12 также не осуществляет самостоятельного движения; тем не менее, поскольку она привязана к обеим клеточным фазам, возникает необходимость указать, каким образом действующая на нее сила распределена между этими фазами. Приняв во внимание, что на микроуровне структурные элементы дополнительной фазы 12 находятся в квазиравновесии и их размер мал в сравнении с характерными размерами задачи, будем полагать, что

(2.1)
${{F}_{{1,\tau 12}}} = {{F}_{{2,\tau 12}}} = {{F}_{{\tau 12}}}$

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

Будем полагать, что силы ${{F}_{{kw}}}$ (k = 1, 2), действующие между жидкой и твердыми (клеточными) фазами, физически сходны с аналогичными силами в любой пористой среде. Для случая одной твердой фазы вопрос подробно исследован и обосновано давно ставшее традиционным представление межфазной силы как суммы вязкой силы, пропорциональной разности скоростей фаз, и силы давления, равной произведению давления в жидкости на градиент ее концентрации [2123]. Обобщим соответствующую формулу на случай двух твердых фаз, предположив пропорциональное концентрациям распределение силы давления между твердыми фазами и отсутствие перекрестных межфазовых эффектов для вязкой силы:

(2.2)
${{F}_{{1w}}} = - p\frac{{{{\phi }_{1}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}\nabla {{\phi }_{w}} + {{{\mathbf{k}}}_{1}} \cdot (w - {{v}_{1}}),\quad {{F}_{{2w}}} = - p\frac{{{{\phi }_{2}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}\nabla {{\phi }_{w}} + {{{\mathbf{k}}}_{2}} \cdot (w - {{v}_{2}})$

Здесь ${{v}_{1}}$ и ${{v}_{2}}$ – скорости первой и второй клеточных фаз соответственно, w – скорость жидкости, а ${{{\mathbf{k}}}_{1}}$ и ${{{\mathbf{k}}}_{2}}$ – тензоры вязких коэффициентов. В частном случае постоянной плотности одной из двух твердых фаз выражение такой структуры использовалось в [15].

Наибольшую сложность представляет получение выражения для силы F12 с компонентами F12i, действующей между двумя клеточными (твердыми) фазами. Для рассматриваемой среды существенна сильная неоднородность напряжений на микроуровне, возникающая вследствие различной механической активности клеток разного типа и проявляющаяся, в частности, в существенном несовпадении средних напряжений в клеточных фазах. Вопрос о выражении для межфазной силы в такой среде в литературе не рассматривался. Можно говорить лишь о некоторых частичных аналогиях [24, 25].

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

${{p}_{1}} = - \frac{1}{3}{{I}_{1}}({{\sigma }^{{(1)}}})$ и ${{p}_{2}} = - \frac{1}{3}{{I}_{1}}({{\sigma }^{{(2)}}})$
где ${{I}_{1}}({{\sigma }^{{(1)}}})$ и ${{I}_{1}}({{\sigma }^{{(2)}}})$ – первые инварианты соответствующих тензоров.

Вследствие существенной неоднородности напряжений среднее объемное давление в фазе не совпадает со средним давлением в ней на межфазной поверхности. Микроскопические давления в фазах на межфазной границе равны между собой, вследствие чего совпадают и осредненные давления на границе раздела фаз. Давление на межфазной поверхности будем обозначать ${{p}_{*}}$.

Рассмотрим малый представительный объем среды V с границей Σ. Объем состоит из трех объемов: $V = {{V}_{1}} + {{V}_{2}} + {{V}_{w}}$ (слагаемые соответствуют двум клеточным фазам и жидкости соответственно). Пусть сначала жидкая фаза отсутствует: $V = {{V}_{1}} + {{V}_{2}}$. Таким образом, ${{\varphi }_{1}} + {{\varphi }_{2}} = 1$, где ${{\varphi }_{1}}$ и ${{\varphi }_{2}}$ – концентрации клеточных фаз. Объемы предполагаются связными. Рассмотрим осреднение микроскопического распределения давлений в фазах в два этапа. Сначала будем осреднять микроскопические переменные в масштабе характерного размера дисперсности, т.е. в таком масштабе, на котором фазы присутствуют раздельно. Второй этап осреднения – по объему $V$, размер которого много больше. Величины, полученные в результате “мелкомасштабного” осреднения, для основной (“полномасштабной”) процедуры осреднения являются микроскопическими. Будем помечать такие величины штрихом.

Связь между давлениями, средними в масштабе дисперсности по объему фазы $p_{k}^{'}$ (k = 1, 2) и по ее границе с другой фазой $p_{*}^{'}$, определяется распределением давлений внутри этой фазы. Распределение давлений внутри включений неизвестно, но естественно предположить, что чем больше объем частиц k-й фазы (т.е. больше ее локальная объемная доля $\varphi _{k}^{'}$), тем бо́льшая часть межфазового перепада средних давлений ${\text{|}}p_{1}^{'} - p_{2}^{'}{\text{|}}$ приходится на эту фазу. Предположим, что доля перепада средних давлений, приходящаяся на каждую из фаз 1 и 2 (${\text{|}}p_{k}^{'} - p_{*}^{'}{\text{|}}$, k = 1, 2), пропорциональна одинаковой для обеих фаз функции их концентраций. Тогда, положив для определенности, что $p_{1}^{'} > p_{2}^{'}$ и учтя, что общее для обеих фаз давление на границе между ними имеет некоторое промежуточное значение ($p_{1}^{'} > p_{*}^{'} > p_{2}^{'}$), имеем зависимость

$\frac{{p_{1}^{'} - p_{*}^{'}}}{{p_{*}^{'} - p_{2}^{'}}} = \frac{{\psi (\varphi _{1}^{'})}}{{\psi (\varphi _{2}^{'})}}$
где ψ – возрастающая функция, стремящаяся к нулю при стремлении к нулю ее аргумента. Отсюда получаем для $p_{*}^{'}$ выражение

(2.3)
$p_{*}^{'} = \frac{{p_{1}^{'}\psi (\varphi _{2}^{'}) + p_{2}^{'}\psi (\varphi _{1}^{'})}}{{\psi (\varphi _{1}^{'}) + \psi (\varphi _{2}^{'})}}$

Далее рассматриваем полномасштабное осреднение по всему объему V, размер которого много больше характерного размера мелкомасштабного осреднения. При этом осредняются величины, полученные в результате этого первоначального осреднения. Сглаженные давления $p_{1}^{'}$ и $p_{2}^{'}$ определены в каждой точке объема. В соответствии с этим в каждой точке по формуле (2.3) формально определено и межфазное давление $p_{*}^{'}$.

Сила давления $P_{{12i}}^{'}$, действующая в объеме V со стороны второй фазы на первую равна

(2.4)
$P_{{12i}}^{'} = - \int\limits_{{{\Sigma }_{{12}}}} {p_{*}^{'}{{n}_{i}}} ds = - \int\limits_{{{\Sigma }_{{12}}} + {{\Sigma }_{1}}} {p_{*}^{'}{{n}_{i}}} ds + \int\limits_{{{\Sigma }_{1}}} {p_{*}^{'}{{n}_{i}}} ds$
где Σ1 – часть внешней границы Σ, приходящаяся на первую фазу, а ni – компоненты вектора нормали к соответствующим поверхностям.

Предполагая, что доля фаз в площади внешней границы асимптотически равна их объемным концентрациям и преобразуя в (2.4) поверхностные интегралы к объемным, а затем асимптотически преобразуя интеграл по объему V1 к интегралу по объему V, получаем

(2.5)
$P_{{12i}}^{'} = - \int\limits_{{{\Sigma }_{{12}}} + {{\Sigma }_{1}}} {p_{*}^{'}{{n}_{i}}ds} + \int\limits_\Sigma {p_{*}^{'}\varphi _{1}^{'}{{n}_{i}}ds} = - \int\limits_{{{V}_{1}}} {{{\nabla }_{i}}p_{*}^{'}d{v}} + \int\limits_V {{{\nabla }_{i}}(p_{*}^{'}\varphi _{1}^{'})d{v}} = \int\limits_V {p_{*}^{'}{{\nabla }_{i}}\varphi _{1}^{'}d{v}} $

Относя окончательное выражение (2.5) к объему V и принимая, что среднее от произведения равно произведению средних, получаем следующее выражение для осредненной по объему $V$ межфазной силы $P_{{12i}}^{{}}$:

(2.6)
${{P}_{{12i}}} = {{p}_{*}}{{\nabla }_{i}}{{\varphi }_{1}}$
связывающее уже средние по объему межфазное давление и фазовую концентрацию.

Примем, что связь между тремя учитываемыми давлениями (2.3) сохраняется для полномасштабных средних:

(2.7)
${{p}_{*}} = \frac{{{{p}_{1}}\psi ({{\varphi }_{2}}) + {{p}_{2}}\psi ({{\varphi }_{1}})}}{{\psi ({{\varphi }_{1}}) + \psi ({{\varphi }_{2}})}}$

Подставляя соотношение (2.7) в (2.6), получаем для межфазной силы зависимость от давлений ${{p}_{1}}$ и ${{p}_{2}}$, которая при ${{p}_{1}} = {{p}_{2}} = p$ для любого вида функции $\psi $ тождественна традиционной зависимости этой силы в заполненной жидкостью пористой среде от давления в жидкости:

${{P}_{{12i}}} = p{{\nabla }_{i}}{{\varphi }_{1}} = - p{{\nabla }_{i}}{{\varphi }_{2}}.$

Ограничимся максимально простым предположением о линейности функции ψ. Тогда, учитывая, что $\psi ({{\varphi }_{k}}) \to 0$ при ${{\varphi }_{k}} \to 0$, имеем зависимость $\psi ({{\varphi }_{k}}) = m{{\varphi }_{k}}$ (m – положительная константа), а выражение (2.7) с учетом равенства ${{\varphi }_{1}} + {{\varphi }_{2}} = 1$ преобразуется к виду

(2.8)
${{p}_{*}} = {{\varphi }_{2}}{{p}_{1}} + {{\varphi }_{1}}{{p}_{2}}$

Подставляя (2.8) в (2.6) и приводя полученное выражение к форме, в которой градиенты концентраций фаз присутствуют равноправно, получаем

(2.9)
${{P}_{{12i}}} = {{\varphi }_{2}}{{p}_{1}}{{\nabla }_{i}}{{\varphi }_{1}} - {{\varphi }_{1}}{{p}_{2}}{{\nabla }_{i}}{{\varphi }_{2}}$

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

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

${{\varphi }_{1}} = \frac{{{{\phi }_{1}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}},\quad {{\varphi }_{2}} = \frac{{{{\phi }_{2}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}$
межфазная сила давления будет падать при росте концентрации жидкости из-за уменьшения площадей контакта между клетками. Приняв зависимость этой силы от концентрации жидкой фазы ϕw в простейшем случае линейной и учитывая, что при ${{\phi }_{w}} = 0$ имеет место формула (2.9), получаем, что сила пропорциональна $1 - {{\phi }_{w}} = {{\phi }_{1}} + {{\phi }_{2}}$:

(2.10)
${{P}_{{12i}}} = ({{\phi }_{1}} + {{\phi }_{2}})({{\varphi }_{2}}{{p}_{1}}{{\nabla }_{i}}{{\varphi }_{1}} - {{\varphi }_{1}}{{p}_{2}}{{\nabla }_{i}}{{\varphi }_{2}}) = {{\phi }_{2}}{{p}_{1}}{{\nabla }_{i}}\left( {\frac{{{{\phi }_{1}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}} \right) - {{\phi }_{1}}{{p}_{2}}{{\nabla }_{i}}\left( {\frac{{{{\phi }_{2}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}} \right)$

Для случая равных давлений в твердых фазах (${{p}_{1}} = {{p}_{2}}$) формулы, сходные по структуре (но не вполне совпадающие) с (2.10), ранее предлагались без детального обоснования и использовались при моделировании биологических сред в [24, 25]).

Для полной силы взаимодействия двух клеточных фаз ${{F}_{{12i}}}$, состоящей из силы давления и вязкой компоненты, имеем

(2.11)
${{{\mathbf{F}}}_{{12}}} = {{{\mathbf{P}}}_{{12i}}} + {{{\mathbf{k}}}_{{12}}} \cdot ({{{\mathbf{v}}}_{2}} - {{v}_{1}}) = {{\phi }_{2}}{{p}_{1}}\nabla \left( {\frac{{{{\phi }_{1}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}} \right) - {{\phi }_{1}}{{p}_{2}}\nabla \left( {\frac{{{{\phi }_{2}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}} \right) + {{{\mathbf{k}}}_{{12}}} \cdot ({{{\mathbf{v}}}_{2}} - {{v}_{1}})$
где k12 – тензор вязких коэффициентов при проскальзывании клеточных фаз одна по другой.

В формулах (2.2) и (2.11) вязкие коэффициенты k1, k2 и k12 далее будут предполагаться шаровыми (в соответствии с общим допущением об изотропии среды) и пропорциональными плотностям взаимодействующих фаз:

(2.12)
${{{\mathbf{k}}}_{1}} = {{k}_{1}}{{\phi }_{1}}{{\phi }_{w}}{\mathbf{I}},\quad {{{\mathbf{k}}}_{2}} = {{k}_{2}}{{\phi }_{2}}{{\phi }_{w}}{\mathbf{I}},\quad {{{\mathbf{k}}}_{{12}}} = {{k}_{{12}}}{{\phi }_{1}}{{\phi }_{2}}{\mathbf{I}}$

С учетом соотношений (2.1), (2.2), (2.11) и (2.12) система уравнений (1.1)–(1.10) становится замкнутой, если заданы входящие в нее функции и коэффициенты.

3. ПОСТАНОВКА ЗАДАЧИ О РАССЛОЕНИИ КЛЕТОК В БЕСКОНЕЧНОМ ПЛОСКОМ СЛОЕ

В качестве нелинейных функций ${{\Pi }_{1}}$, ${{\Pi }_{2}}$ и ${{\Pi }_{{12}}}$, описывающих хаотическую активность клеток, препятствующую их избыточному сближению, могут быть предложены различные соотношения, которые удовлетворяют сформулированным выше требованиям. Возможной формой таких зависимостей являются дробно-рациональные функции концентраций ${{\phi }_{1}}$ и ${{\phi }_{2}}$, имеющие в знаменателе $1 - {{\phi }_{1}} - {{\phi }_{2}}$, а в числителе однородные функции одинаковой степени. В дальнейшем будем использовать следующие выражения:

(3.1)
${{\Pi }_{1}} = {{E}_{1}}\frac{{\phi _{1}^{2}}}{{1 - {{\phi }_{1}} - {{\phi }_{2}}}},\quad {{\Pi }_{2}} = {{E}_{2}}\frac{{\phi _{2}^{2}}}{{1 - {{\phi }_{1}} - {{\phi }_{2}}}},\quad {{\Pi }_{{12}}} = {{E}_{{12}}}\frac{{{{\phi }_{1}}{{\phi }_{2}}}}{{1 - {{\phi }_{1}} - {{\phi }_{2}}}}$

Сходные по форме с (3.1) зависимости использовались в [15, 24], а также в нашей работе [20].

В уравнениях (1.8) сократим число независимых коэффициентов, предполагая, что влияние на неупругое деформирование клеточной фазы силового воздействия других фаз (жидкой и другой клеточной) пропорционально их концентрациям. Кроме того, примем для простоты, что компоненты ${{\tau }^{{[12]}}}$ и ${{\tau }^{{[21]}}}$ “перекрестного” активного напряжения ${{\tau }^{{(12)}}}$ не влияют на этот процесс раздельно, а лишь через суммарное перекрестное напряжение ${{\tau }^{{(12)}}}$:

$\begin{gathered} {{{\mathbf{e}}}^{{(1)}}} = \frac{1}{{2{{\mu }_{1}}}}[{{\sigma }^{{(1)}}} + {{\nu }_{1}}({{\phi }_{w}}p{\mathbf{I}} - {{\phi }_{2}}{{\sigma }^{{(2)}}})] - \frac{1}{{2{{\mu }_{{11}}}}}{{\tau }^{{(1)}}} - \frac{1}{{2{{\mu }_{{12}}}}}{{\tau }^{{(12)}}} \\ {{{\mathbf{e}}}^{{(2)}}} = \frac{1}{{2{{\mu }_{2}}}}[{{\sigma }^{{(2)}}} + {{\nu }_{2}}({{\phi }_{w}}p{\mathbf{I}} - {{\phi }_{1}}{{\sigma }^{{(1)}}})] - \frac{1}{{2{{\mu }_{{22}}}}}{{\tau }^{{(2)}}} - \frac{1}{{2{{\mu }_{{21}}}}}{{\tau }^{{(12)}}} \\ \end{gathered} $
где все коэффициенты положительны. В последующих расчетах примем ${{\nu }_{1}} = {{\nu }_{2}} = 1$.

Рассмотрим применение сформулированных выше уравнений к решению наиболее элементарной в механическом отношении задачи о сортировке двух типов клеток в бесконечном плоском слое.

Пусть в начальный момент времени клетки двух разных типов имеют пространственно-однородные распределения и заполняют плоский слой толщиной $2{{h}_{0}}$, окруженный жидкостью, в которой давление равно нулю. Толщина слоя h в дальнейшем отождествляется с областью, занятой клетками, и, в зависимости от конкретной постановки задачи, остается постоянной ($h = {{h}_{0}}$) или меняется со временем ($h = h(t)$). Введем декартову систему координат с осями Oy и Oz, лежащими в средней плоскости слоя и осью Ox, направленной перпендикулярно этой плоскости.

Решение задачи о перераспределении объемных плотностей клеток и жидкости в слое будем искать в следующем виде: ${{\phi }_{1}} = {{\phi }_{1}}(t,\;x)$, ${{\phi }_{2}} = {{\phi }_{2}}(t,\;x)$, $p = p(t,\;x)$. Будем полагать, что скорости фаз имеют единственную, отличную от нуля компоненту, направленную по оси Ox: $w = w(t,\;x)$, ${{{v}}_{1}} = {{{v}}_{1}}(t,x)$, ${{{v}}_{2}} = {{{v}}_{2}}(t,x)$ для межклеточной жидкости и клеточных фаз 1 и 2 соответственно. У всех тензоров напряжений смешанные компоненты считаются равными нулю, а диагональные компоненты – функциями t и x. Будем отмечать компоненты векторов и тензоров в соответствующих направлениях индексами x, y, z, двукратный индекс у компонент тензоров заменять однократным и опускать индекс x. После исключения ${{F}_{{1,\tau 1}}}$, ${{F}_{{1,\tau 12}}}$, ${{F}_{{2,\tau 2}}}$ и ${{F}_{{2,\tau 12}}}$ с помощью последних трех уравнений (1.7) и условия (2.1) уравнения (1.1)(1.10), (2.2), (2.11), (2.12) принимают следующий вид:

$\frac{{\partial {{\phi }_{1}}}}{{\partial t}} + \frac{{\partial {{\phi }_{1}}{{{v}}_{1}}}}{{\partial x}} = 0,\quad \frac{{\partial {{\phi }_{2}}}}{{\partial t}} + \frac{{\partial {{\phi }_{2}}{{{v}}_{2}}}}{{\partial x}} = 0,\quad \frac{{\partial {{\phi }_{w}}}}{{\partial t}} + \frac{{\partial {{\phi }_{w}}w}}{{\partial x}} = 0,\quad {{\phi }_{w}} + {{\phi }_{1}} + {{\phi }_{2}} = 1$
$\begin{gathered} \frac{\partial }{{\partial x}}\left( {{{\phi }_{1}}{{\sigma }^{{(1)}}} + {{\tau }^{{(1)}}} + \frac{1}{2}{{\tau }^{{(12)}}}} \right) + {{\phi }_{2}}{{p}_{1}}\frac{\partial }{{\partial x}}\left( {\frac{{{{\phi }_{1}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}} \right) - {{\phi }_{1}}{{p}_{2}}\frac{\partial }{{\partial x}}\left( {\frac{{{{\phi }_{2}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}} \right) + {{k}_{{12}}}{{\phi }_{1}}{{\phi }_{2}}({{{v}}_{2}} - {{{v}}_{1}}) - \\ - \;p\frac{{{{\phi }_{1}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}\frac{{\partial {{\phi }_{w}}}}{{\partial x}} + {{k}_{1}}{{\phi }_{1}}{{\phi }_{w}}(w - {{{v}}_{1}}) = 0 \\ \end{gathered} $
$\begin{gathered} \frac{\partial }{{\partial x}}\left( {{{\phi }_{2}}{{\sigma }^{{(2)}}} + {{\tau }^{{(2)}}} + \frac{1}{2}{{\tau }^{{(12)}}}} \right) + {{\phi }_{1}}{{p}_{2}}\frac{\partial }{{\partial x}}\left( {\frac{{{{\phi }_{2}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}} \right) - {{\phi }_{2}}{{p}_{1}}\frac{\partial }{{\partial x}}\left( {\frac{{{{\phi }_{1}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}} \right) + {{k}_{{12}}}{{\phi }_{1}}{{\phi }_{2}}({{{v}}_{1}} - {{{v}}_{2}}) - \\ - \;p\frac{{{{\phi }_{2}}}}{{{{\phi }_{1}} + {{\phi }_{2}}}}\frac{{\partial {{\phi }_{w}}}}{{\partial x}} + {{k}_{2}}{{\phi }_{2}}{{\phi }_{w}}(w - {{{v}}_{2}}) = 0 \\ \end{gathered} $
$ - \frac{{\partial p}}{{\partial x}} - {{k}_{1}}{{\phi }_{1}}(w - {{{v}}_{1}}) - {{k}_{2}}{{\phi }_{2}}(w - {{{v}}_{2}}) = 0$
${{e}^{{(1)}}} = \frac{1}{{2{{\mu }_{1}}}}({{\sigma }^{{(1)}}} - {{\phi }_{2}}{{\sigma }^{{(2)}}} + {{\phi }_{w}}p) - \frac{1}{{2{{\mu }_{{11}}}}}{{\tau }^{{(1)}}} - \frac{1}{{2{{\mu }_{{12}}}}}{{\tau }^{{(12)}}}$
${{e}^{{(2)}}} = \frac{1}{{2{{\mu }_{2}}}}({{\sigma }^{{(2)}}} - {{\phi }_{1}}{{\sigma }^{{(1)}}} + {{\phi }_{w}}p) - \frac{1}{{2{{\mu }_{{22}}}}}{{\tau }^{{(2)}}} - \frac{1}{{2{{\mu }_{{21}}}}}{{\tau }^{{(12)}}}$
(3.1)
$\frac{1}{{2{{\mu }_{1}}}}(\sigma _{\alpha }^{{(1)}} - {{\phi }_{2}}\sigma _{\alpha }^{{(2)}} + {{\phi }_{w}}p) - \frac{1}{{2{{\mu }_{{11}}}}}\tau _{\alpha }^{{(1)}} - \frac{1}{{2{{\mu }_{{12}}}}}\tau _{\alpha }^{{(12)}} = 0$
$\frac{1}{{2{{\mu }_{2}}}}(\sigma _{\alpha }^{{(2)}} - {{\phi }_{1}}\sigma _{\alpha }^{{(1)}} + {{\phi }_{w}}p) - \frac{1}{{2{{\mu }_{{22}}}}}\tau _{\alpha }^{{(2)}} - \frac{1}{{2{{\mu }_{{21}}}}}\tau _{\alpha }^{{(12)}} = 0$
${{\tau }^{{(1)}}} = - {{\Pi }_{1}} + {{m}_{1}}\frac{3}{{4\pi R_{s}^{4}}} \cdot {{\phi }_{1}}\int\limits_S {{{\phi }_{1}}} \cdot {{x}^{2}}ds,\quad {{\tau }^{{(2)}}} = - {{\Pi }_{2}} + {{m}_{2}}\frac{3}{{4\pi R_{s}^{4}}} \cdot {{\phi }_{2}}\int\limits_S {{{\phi }_{2}}} \cdot {{x}^{2}}ds$
$\tau _{{}}^{{(12)}} = - {{\Pi }_{{12}}} + \frac{3}{{4\pi R_{s}^{4}}}\left( {{{m}_{{12}}} \cdot {{\phi }_{1}}\int\limits_S {{{\phi }_{2}}} \cdot {{x}^{2}}ds + {{m}_{{21}}} \cdot {{\phi }_{2}}\int\limits_S {{{\phi }_{1}}} \cdot {{x}^{2}}ds} \right)$
$\tau _{\alpha }^{{(1)}} = - {{\Pi }_{1}} + {{m}_{1}}\frac{3}{{4\pi R_{s}^{4}}} \cdot {{\phi }_{1}}\int\limits_S {{{\phi }_{1}}} \cdot r_{\alpha }^{2}ds,\quad \tau _{\alpha }^{{(2)}} = - {{\Pi }_{2}} + {{m}_{2}}\frac{3}{{4\pi R_{s}^{4}}} \cdot {{\phi }_{2}}\int\limits_S {{{\phi }_{2}}} \cdot r_{\alpha }^{2}ds$
$\tau _{\alpha }^{{(12)}} = - {{\Pi }_{{12}}} + \frac{3}{{4\pi R_{s}^{4}}}\left( {{{m}_{{12}}} \cdot {{\phi }_{1}}\int\limits_S {{{\phi }_{2}}} \cdot r_{\alpha }^{2}ds + {{m}_{{21}}} \cdot {{\phi }_{2}}\int\limits_S {{{\phi }_{1}}} \cdot r_{\alpha }^{2}ds} \right)$
${{p}_{1}} = - \frac{1}{3}{{I}_{1}}({{\sigma }^{{(1)}}}),\quad {{p}_{2}} = - \frac{1}{3}{{I}_{1}}({{\sigma }^{{(2)}}})$
$e_{{}}^{{(1)}} = \frac{{\partial {{{v}}_{1}}}}{{\partial x}},\quad e_{{}}^{{(2)}} = \frac{{\partial {{{v}}_{2}}}}{{\partial x}}$

Здесь индекс α принимает значения y и z, ${{r}_{y}} \equiv y$, ${{r}_{z}} \equiv z$. Уравнения баланса импульса для каждой фазы по осям Oy и Oz выполнены тождественно.

Распределение плотностей клеток и жидкости на отрезке $[ - h,h]$ считается симметричным, и система уравнений (3.1) решается на отрезке $[0,h]$. Будут исследованы решения задач с граничными условиями, учитывающими два типа нагружения клеточных фаз. В первом случае толщина слоя считается постоянной. Тогда в качестве граничных условий для скоростей примем, что ${{{v}}_{1}} = {{{v}}_{2}} = 0$ при x = 0 и x = h, и $w = 0$ при x = 0. Давление во внеклеточной жидкости удовлетворяет условию p = 0 при x = h. Граничные условия для ${{\phi }_{1}}$ и ${{\phi }_{2}}$ не ставятся, так как x = 0 и $x = h$ являются характеристиками для первого и второго уравнений системы (3.1). Соответствующие граничные условия получаются из решения этих уравнений, записанных в точках x = 0 и x = h с учетом условий ${{{v}}_{1}} = {{{v}}_{2}} = 0$. Однако более удобным при использовании численных методов является условие симметрии объемных плотностей клеток $\partial {{\phi }_{1}}{\text{/}}\partial x = \partial {{\phi }_{2}}{\text{/}}\partial x = 0$ при x = 0. Это условие приводит к соотношению на характеристике в точке x = 0. Новая клеточная структура формируется в результате эволюции некоторого начального распределения клеток. В качестве начальных условий для решения системы (3.1) примем распределения ${{\phi }_{1}}(x) = {{\phi }_{{10}}} + \Delta {{\phi }_{1}}(x)$ и ${{\phi }_{2}}(x) = {{\phi }_{{20}}} + \Delta {{\phi }_{2}}(x)$, где $\Delta {{\phi }_{1}}(x)$ и $\Delta {{\phi }_{2}}(x)$ – малые возмущения некоторых пространственно-однородных распределений ${{\phi }_{1}}(x) = {{\phi }_{{10}}}$ и ${{\phi }_{2}}(x) = {{\phi }_{{20}}}$.

Во втором случае будет решаться задача с подвижными внешними границами $x = \pm h(t)$, перемещение которых определяется условиями нагружения. Будет рассмотрен случай равенства нулю внешней силы. Изменение толщины слоя связывается с перемещениями клеток, для которых предполагается совпадение скоростей двух клеточных фаз на внешних границах $x = \pm h$. В этом случае в качестве граничных условий примем ${{{v}}_{1}} = {{{v}}_{2}} = 0$, $\partial {{\phi }_{1}}{\text{/}}\partial x = \partial {{\phi }_{2}}{\text{/}}\partial x = 0$ при x = 0 и $ - (1 - {{\phi }_{1}} - {{\phi }_{2}})p + {{\phi }_{1}}{{\sigma }^{{(1)}}} + {{\phi }_{2}}{{\sigma }^{{(2)}}} + {{\tau }^{{(1)}}} + {{\tau }^{{(2)}}} + {{\tau }^{{(12)}}} = 0$, ${{{v}}_{1}} = {{{v}}_{2}} = dh{\text{/}}dt$ при x = h. В качестве начальных условий будем использовать соотношение $h = {{h}_{0}}$ и, как и в первой постановке, распределения ${{\phi }_{1}}(x) = {{\phi }_{{10}}} + \Delta {{\phi }_{1}}(x)$, ${{\phi }_{2}}(x) = {{\phi }_{{20}}} + \Delta {{\phi }_{2}}(x)$.

Рассматриваемая постановка соответствует предположению об отсутствии боковых смещений в направлениях осей Oy и Oz на бесконечно удаленных границах. В слое присутствуют клеточные и активные напряжения, параллельные слою, которые могут быть рассчитаны по тем уравнениям системы (3.1), которые содержат компоненты с индексом $\alpha $. В процессе перераспределения клеток в задаче с подвижными границами жидкость может перетекать через границу слоя (т.е. области, занятой клетками).

Преобразование координаты $x* = x{\text{/}}h(t)$ позволяет перейти от решения задачи на отрезке [0, h] (в общем случае с подвижной верхней границей) к решению задачи на отрезке [0, 1]. При введении безразмерных величин в качестве масштаба длины используется начальная толщина слоя ${{h}_{0}}$. Безразмерное время вводится как $t* = tE{\text{/}}\mu $, где присутствуют некоторые характерные значения коэффициентов $E$ и μ (с разными индексами). Те же характерные значения E и μ используются при определении безразмерных напряжений и скоростей.

Коэффициент межфазного трения между клеточной фазой и жидкостью оценивался как 5 × 107 Па ∙ с ∙ м–2, а коэффициента вязкости для эмбриональных тканей как 1 × 104 Па ∙ с [27]. Тогда для многоклеточного биологического объекта с характерным размером порядка 100 мкм [28] значения безразмерных коэффициентов межфазного трения между клеточными фазами и жидкостью $k_{1}^{*} = {{k}_{1}}h_{0}^{2}{\text{/}}\mu $ и $k_{2}^{*} = {{k}_{2}}h_{0}^{2}{\text{/}}\mu $ составят величину порядка 5 × 10–5. В соответствии с этими оценками в уравнении системы (3.1), определяющем $\partial p{\text{/}}\partial x$, можно пренебречь изменением гидростатического давления и считать его постоянным и равным нулю.

Исключая напряжения в клеточных фазах из уравнений импульсов с помощью уравнений для ${{e}^{{(1)}}}$, ${{e}^{{(2)}}}$ и кинематических соотношений и пренебрегая силами вязкого взаимодействия между клеточными фазами и внеклеточной жидкостью, получим систему уравнений для ${{\phi }_{1}}$, ${{\phi }_{2}}$ (типа уравнений переноса в дивергентной форме) и ${{{v}}_{1}}$, ${{{v}}_{2}}$ (обыкновенные дифференциальные уравнения второго порядка).

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Исследуем сначала эволюцию плотности клеток в бесконечном плоском слое, образованном клетками, для которых активные взаимодействия внутри фаз доминируют над межфазными активными взаимодействиями. При этом активные взаимодействия внутри клеточной фазы 1 превышают таковые внутри клеточной фазы 2. В последующих расчетах, если не оговорено иное, приняты следующие значения безразмерных параметров, характеризующих напряжения, развивающиеся в результате активных межклеточных взаимодействий: $m_{1}^{*} = {{m}_{1}}{\text{/}}E = 12$, $m_{2}^{*} = {{m}_{2}}{\text{/}}E$ = 6, $m_{{12}}^{*} = {{m}_{{12}}}{\text{/}}E = m_{{21}}^{*} = {{m}_{{21}}}{\text{/}}E$ = 0. Для остальных безразмерных параметров приняты значения:

$\mu _{k}^{*} = \frac{{{{\mu }_{k}}}}{\mu } = 1,\quad E_{k}^{*} = \frac{{{{E}_{k}}}}{E} = 1,\quad \mu _{{km}}^{*} = \frac{{{{\mu }_{{km}}}}}{\mu } = 1,\quad E_{{km}}^{*} = \frac{{{{E}_{{km}}}}}{E} = 1\quad (k = 1,2;m = 1,2)$
$k_{{12}}^{*} = \frac{{{{k}_{{12}}}h_{0}^{2}}}{\mu } = 0.3,\quad {{\phi }_{{10}}} = {{\phi }_{{20}}} = {{\phi }_{0}} = 0.3$

В дальнейшем звездочку при безразмерных параметрах будем опускать, но сохраним ее при безразмерных координате и времени. Существенным параметром, влияющим на эволюцию начальных распределений плотностей клеток, является радиус дальнодействия ${{R}_{s}}$, определяющий характерное расстояние, на котором клетки взаимодействуют одна с другой. Считая ${{R}_{S}}$ величиной порядка характерного радиуса клетки, будем рассматривать плоский слой с начальной толщиной порядка100 клеточных диаметров, т.е. положим ${{R}_{S}} = 0.01$. Значение толщины слоя соответствует характерному размеру кластера клеток, получаемого в экспериментах смешиванием клеток двух типов [3, 28].

При расчетах начальные распределения плотностей клеток принимаются либо в виде ${{\phi }_{1}}(0,x) = {{\phi }_{0}} + 0.001 \times (2 \times {\text{random(}}x{\text{)}} - 1)$ и ${{\phi }_{2}}(0,x) = {{\phi }_{0}} + 0.001 \times (2 \times {\text{random(}}x{\text{)}} - 1)$, где функция random(x) генерирует значения случайной величины, равномерно распределенной на отрезке [0; 1], либо в форме регулярных синусоидальных распределений с различным периодом.

Решение системы уравнений для ${{\phi }_{1}}$, ${{\phi }_{2}}$ и ${{{v}}_{1}}$, ${{{v}}_{2}}$ на каждом временно́м шаге начинается с решения уравнений для ${{{v}}_{1}}$, ${{{v}}_{2}}$; при этом значение ${{\phi }_{1}}$ и ${{\phi }_{2}}$ берутся с предыдущего временного шага. Разностная схема, полученная интегро-интерполяционным методом [29], монотонна и имеет первый порядок аппроксимации по x. Монотонность достигается использованием направленных разностей при аппроксимации $\partial {{{v}}_{1}}{\text{/}}\partial x$ и $\partial {{{v}}_{2}}{\text{/}}\partial x$ в зависимости от знаков стоящих при них коэффициентов. Полученная система разностных уравнений, позволяющая найти ${{{v}}_{1}}$ и ${{{v}}_{2}}$, решалась методом матричной прогонки.

Затем при полученных значениях ${{{v}}_{1}}$ и ${{{v}}_{2}}$ по отдельности решались уравнения для ${{\phi }_{1}}$ и ${{\phi }_{2}}$. При аппроксимации дивергентных слагаемых использовалась двухслойная разностная схема [30] c весами 0.5. Эта неявная схема обладает свойством консервативности, однако не является монотонной. Поэтому корректность вычислений проверялась путем сравнения решений при уменьшении шага по пространственной переменной $x$. Для граничных узлов мы ограничивались аппроксимацией направленными разностями; поэтому разностные уравнения имеют порядок аппроксимации $O(\Delta x) + O(\Delta t)$. Уравнения решались методом прогонки. Вычисление интегралов в выражениях для активных напряжений проводилось приближенно с помощью разложения подынтегральных функций в ряд по формуле Тейлора с точностью до слагаемых второго порядка. При этом учитывалось различие пределов интегрирования для внутренних и приграничных точек, т.е. точек, удаленных от границы слоя на расстояние, не превышающее радиус дальнодействия.

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

Эволюция распределений плотностей клеток, полученных в результате решения задачи о сортировке клеток в слое с фиксированной границей (постоянной толщины), показана на рис. 1. Видно, что клеточная фаза 1, в которой активные взаимодействия между клетками сильнее, чем во второй (${{m}_{1}} > {{m}_{2}}$), стремится занять центральную область слоя. Сформировалась отделенная резкой границей центральная область, в которой плотность клеток фазы 1 значительно превосходит начальную, тогда как вне этой области она быстро падает практически до нуля. Периферическая область занята практически только клетками фазы 2; плотность клеток этой фазы возрастает при приближении к границе. Резкое изменение плотности фазы 2 около внешней границы x* = 1 связано с выбором начального возмущения в результате реализации значений случайной величины.

Рис. 1.

Распределения плотностей клеточных фаз (${{m}_{1}} > {{m}_{2}}$) в задаче с фиксированной границей при $t* = 1.89$.

Распределения плотности клеток в момент времени $t* = 1.1$ и изменение толщины слоя при подвижной границе и отсутствии внешней нагрузки представлены на рис. 2. Полученные результаты, как и в задаче с фиксированной границей, описывают тенденцию клеточной фазы, характеризуемой более сильными активными взаимодействиями между клеток, занять центральную область и быть окруженной клетками другой фазы (рис. 2а). При этом толщина слоя уменьшается (рис. 2б) и суммарная плотность клеток обеих фаз возрастает. В последующие моменты времени распределения плотности клеток стремятся принять вид, аналогичный рис. 1.

Рис. 2.

Распределения плотностей клеточных фаз по координате $x* = x{\text{/}}h(t)$ (${{m}_{1}} > {{m}_{2}}$) в задаче с подвижной границей при $t* = 1.1$ (а) и эволюция толщины слоя по отношению к ее начальному значению (б).

Важным параметром, оказывающим значительное влияние на сортировку клеток и формирование новых клеточных паттернов, является радиус дальнодействия ${{R}_{S}}$, отвечающий за степень нелокальности взаимодействия между клетками. На рис. 3 представлены предельные при ${{R}_{S}} \to 0$ распределения плотностей клеток для задач с фиксированной и подвижной границами. В этом случае интеграл заменяется значением подынтегральной функции в рассматриваемой точке, т.е. взаимодействие локально. Начальные распределения плотностей клеток взяты в следующем виде: ${{\phi }_{1}}(0,x) = {{\phi }_{0}} + 0.001 \times \cos \left( {10\pi \times x} \right)$ и ${{\phi }_{2}}(0,x) = {{\phi }_{0}} - 0.001 \times \cos \left( {10\pi \times x} \right)$.

Рис. 3.

Распределения плотностей клеточных фаз (${{m}_{1}} > {{m}_{2}}$) в результате эволюции начальных распределений ${{\phi }_{1}}(0,x) = {{\phi }_{0}} + 0.001 \times \cos \left( {10\pi \times x} \right)$ и ${{\phi }_{2}}(0,x) = {{\phi }_{0}} - 0.001 \times \cos \left( {10\pi \times x} \right)$ в задаче с неподвижной границей при $t* = 1.89$ (а), и подвижной границей при $t* = 1.4$ (б).

Представленные на рис. 3 графики позволяют заключить, что нелокальность активных взаимодействий между клетками является необходимым условием, позволяющим получать распределения плотностей клеток, качественно описывающие их сортировку. Локальные законы активных межклеточных взаимодействий (при ${{R}_{S}} \to 0$) не приводят к сортировке. В случае регулярных синусоидальных начальных возмущений возникает, как видно на рис. 3, чередование областей с различным соотношением между плотностями клеток разных типов. Расчеты, проведенные для случайных начальных возмущений, приводили к формированию неупорядоченных паттернов.

ЗАКЛЮЧЕНИЕ

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

Работа поддержана РФФИ (проект № 20-01-00329).

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

  1. Steinberg M.S., Wiseman L.L. Do morphogenetic tissue rearrangements require active cell movements? // J. Cell Biol. 1972. V. 55. P. 606–615.

  2. Mehes E., Viscek T. Segregation mechanisms of tissue cells: from experimental data to models // Complex Adapt. Syst. Model. 2013. V. 1. P. 4.

  3. Mehes E., Viscek T. Collective motion of cells: from experiments to models // Integr. Biol. 2014. V. 6. № 9. P. 831–854.

  4. Painter K.J., Bloomfield J.M., Sherratt J.A.,·Gerisch A. A nonlocal model for contact attraction and repulsion in heterogeneous cell populations // Bull. Math. Biol. 2015. V. 77. P. 1132–1165.

  5. Pawlizak S., Fritsch A.W., Grosser S., Ahrens D., Thalheim T., Riedel S., Kiebling T.R., Oswald L., Zink M., Manning M.E., Kas J.A. Testing the differential adhesion hypothesis across the epithelial-mesenchymal transition // New Journal of Physics. 2015. V. 17. № 8. 083049.

  6. Armstrong N.J., Painter K.J., Sherratt J.A. A continuum approach to modelling cell–cell adhesion // J. Theor. Biol. 2006. V. 243. № 1. P. 98–113.

  7. Gerisch A., Chaplain M.A.J. Mathematical modelling of cancer cell invasion of tissue: Local and non-local models and the effect of adhesion // J. Theor. Biol. 2008. V. 250. № 4. P. 684–704.

  8. Painter K.J., Armstrong N.J., Sherratt J.A. The impact of adhesion on cellular invasion processes in cancer and development // J. Theor. Biol. 2010. V. 264. № 3. P. 1057–1067.

  9. Domschke P., Trucu D., Gerisch A., Chaplain M. Mathematical modelling of cancer invasion: Implications of cell adhesion variability for tumour infiltrative growth patterns // J. Theor. Biol. 2014. V. 361. P. 41–60.

  10. Murakawa H., Togashi H. Continuous models for cell–cell adhesion // J. Theor. Biol. 2015. V. 374. P. 1–12.

  11. Bitsouni V., Trucu D., Chaplain M., Eftimie R. Aggregation and travelling wave dynamics in a two-population model of cancer cell growth and invasion // Math. Med. Biol. 2018. V. 35. № 4. P. 541–547.

  12. Jackson T.L., Byrne H.M. A mechanical model of tumor encapsulation and transcapsular spread // Math. Biosci. 2002. V. 180. P. 307–328.

  13. Green J.E., Waters S.L., Shakesheff K.M., Byrne H.M. A mathematical model of liver cell aggregation in vitro // Bull. Math. Biol. 2009. V. 71. P. 906–930.

  14. Lemon G., King J.R., Byrne H.M., Jensen O.E., Shakesheff K.M. Mathematical modelling of engineered tissue growth using a multiphase porous flow mixture theory // J. Math. Biol. 2006. V. 52. P. 571–594.

  15. O’Dea R.D., Waters S.L., Byrne H.M. A multiphase model for tissue construct growth in a perfusion bioreactor // Math. Med. Biol. 2010. V. 27. № 2. P. 95–127.

  16. Stein A.A., Logvenkov S.A., Volodyaev I.V. Continuum modeling of mechano-dependent reactions in tissues composed of mechanically active cells // BioSystems. 2018. V. 173. P. 225–234.

  17. Белоусов Л.В., Логвенков С.А., Штейн А.А. Математическая модель активной биологической сплошной среды с учетом деформаций и переупаковки клеток // Известия РАН. МЖГ. 2015. № 1. С. 3–14.

  18. Логвенков С.А., Штейн А.А. Математическое моделирование индуцированного растяжением удлинения слоя эмбрионального эпителия при отсутствии внешней нагрузки // Биофизика. 2015. Т. 60. № 6. С. 1174–1179.

  19. Логвенков С.А., Штейн А.А. Математическая модель пространственной самоорганизации в механически активной клеточной среде // Биофизика. 2017. Т. 62. № 6. С. 1123–1133.

  20. Логвенков С.А. Математическая модель биологической среды с учетом активных взаимодействий и взаимных перемещений составляющих ее клеток // Изв. РАН. МЖГ. 2018. № 5. С. 3–16.

  21. Нигматулин Р.И. Основы механики гетерогенных сред. М: Наука, 1978. 336 с.

  22. Whitaker S. Flow in Porous Media I: A Theoretical Derivation of Darcy’s Law // Transport in Porous Media. 1986. V. 1. P. 3–25.

  23. Drew D.A., Segel L.A. Averaged equations for two-phase flows // Stud. Appl. Math. 1971. V. 50. № 3. P. 205–231.

  24. Lemon G., King J.R., Byrne H.M., Jensen O.E., Shakesheff K.M. Mathematical modelling of engineered tissue growth using a multiphase porous flow mixture theory // J. Math. Biol. 2006 V. 52. P. 571–594.

  25. Preziosi L., Tosin A. Multiphase modeling of tumor growth and extracellular matrix interaction: Mathematical tools and applications // J. Math. Biol. 2009. V. 58. P. 625–656.

  26. Green A.E., Naghdi P.M. A dynamical theory of interacting continua // Int. J. Engng. Sci. 1965. V. 3. P. 231–241.

  27. Lemon G., King J.R. Multiphase modelling of cell behaviour on artificial scaffolds: effects of nutrient depletion and spatially nonuniform porosity // Math. Med. Biol. 2007. V. 24. P. 57−83.

  28. Foty R.A., Steinberg M.S. Cadherin-mediated cell-cell adhesion and tissue segregation in relation to malignancy // Int. J. Dev. Biol. 2004. V. 48. P. 397–409.

  29. Самарский А.А. Теория разностных схем. М.: Наука, 1977. 656 с.

  30. Самарский А.А., Вабищевич П.Н. Разностные схемы для уравнения переноса // Диф. уравн. 1998. Т. 34. № 12. С. 1675–1685.

Дополнительные материалы отсутствуют.