Derivation and application of effective interface conditions for continuum mechanical models of cell invasion through thin membranes

TL;DR

This paper develops a mathematical framework to simplify models of cell invasion through thin membranes by replacing them with effective interfaces, validated through simulations and applied to cancer invasion scenarios.

Contribution

It introduces a formal asymptotic method to derive biophysically consistent interface conditions for continuum models of cell invasion through thin membranes.

Findings

The derived interface conditions accurately approximate the original problem as membrane thickness approaches zero.

Numerical simulations confirm the validity of the limiting transmission problem.

Application examples demonstrate the model's relevance to cancer cell invasion and metastasis.

Abstract

We consider a continuum mechanical model of cell invasion through thin membranes. The model consists of a transmission problem for cell volume fraction complemented with continuity of stresses and mass flux across the surfaces of the membranes. We reduce the original problem to a limiting transmission problem whereby each thin membrane is replaced by an effective interface, and we develop a formal asymptotic method that enables the derivation of a set of biophysically consistent transmission conditions to close the limiting problem. The formal results obtained are validated via numerical simulations showing that the relative error between the solutions to the original transmission problem and the solutions to the limiting problem vanishes when the thickness of the membranes tends to zero. In order to show potential applications of our effective interface conditions, we employ the limiting transmission problem to model cancer cell invasion through the basement membrane and the metastatic spread of ovarian carcinoma.

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=\Gamma (\rho ),(t,\mathbf{x})\in {\mathbb{R}}^{+}\times \mathcal{D}$$

$$\mathbf{v}:=-\mu \nabla p,$$

$$[[\rho \mathbf{v}\cdot {\mathbf{n}}_{ij}]]=0\text{mbox}on{\Sigma }_{ij}\text{ with }i=1,2\text{ and }j=i+1,$$

$$[[p]]=0\text{mbox}on{\Sigma }_{ij}\text{ with }i=1,2\text{ and }j=i+1.$$

$$[[\mu \rho \nabla p\cdot {\mathbf{n}}_{ij}]]=0\text{mbox}on{\Sigma }_{ij}\text{ with }i=1,2\text{ and }j=i+1.$$

$$\left\{\begin{array}[]{ll}\displaystyle{\frac{\partial\rho_{i}}{\partial t}-\nabla\cdot(\mu_{i}\,\rho_{i}\,\nabla p)=\Gamma_{i}(\rho_{i})}&{\rm in}\ {\cal D}_{i},\;i=1,2,3,\\[12.0pt] \mu_{i}\,\rho_{i}\,\nabla p\cdot{\bf n}_{ij}=\mu_{j}\,\rho_{j}\,\nabla p\cdot{\bf n}_{ij}&{\rm on}\ \Sigma_{ij},\;i=1,2,\\[12.0pt] [\![p]\!]=0&{\rm on}\ \Sigma_{ij},\;i=1,2,\end{array}\right.$$

$$[[\rho ]]=0\text{mbox}on{\Sigma }_{ij}\text{ with }i=1,2\text{ and }j=i+1.$$

$$[[\mathbf{v}\cdot {\mathbf{n}}_{ij}]]=0\text{or}[[\mu \nabla p\cdot {\mathbf{n}}_{ij}]]=0\text{mbox}on{\Sigma }_{ij}\text{ with }i=1,2\text{ and }j=i+1.$$

$$\varepsilon:=\max_{\hat{{\bf x}}_{12}\in\Sigma_{12}}\big{\{}\min\{a>0:\hat{{\bf x}}_{12}+a\,{\bf n}_{12}\in\Sigma_{23}\}\big{\}}$$

$${\cal P}_{\varepsilon}\equiv\left\{\begin{array}[]{ll}\displaystyle{\frac{\partial\rho_{i\varepsilon}}{\partial t}-\nabla\cdot(\mu_{i\varepsilon}\rho_{i\varepsilon}f^{\prime}(\rho_{i\varepsilon})\nabla\rho_{i\varepsilon})=\Gamma_{i\varepsilon}(\rho_{i\varepsilon})}&{\rm in}\ {\cal D}_{i\varepsilon},\;i=1,2,3,\\[12.0pt] \mu_{i\varepsilon}\nabla\rho_{i\varepsilon}\cdot{\bf n}_{ij}=\mu_{j\varepsilon}\nabla\rho_{j\varepsilon}\cdot{\bf n}_{ij}&{\rm on}\ \Sigma_{ij\varepsilon},\;i=1,2,\\[12.0pt] \rho_{i\varepsilon}=\rho_{j\varepsilon}&{\rm on}\ \Sigma_{ij\varepsilon},\;i=1,2,\end{array}\right.$$

$${\cal P}_{0}\equiv\left\{\begin{array}[]{ll}\displaystyle{\frac{\partial\tilde{\rho}_{i}}{\partial t}-\nabla\cdot(\tilde{\mu}_{i}\,\tilde{\rho}_{i}\,f^{\prime}(\tilde{\rho}_{i})\nabla\tilde{\rho}_{i})=\tilde{\Gamma}_{i}(\tilde{\rho}_{i})}&{\rm in}\ \tilde{\cal D}_{i},\;i=1,3,\\[12.0pt] \text{transmission conditions}&{\rm on}\ \tilde{\Sigma}_{13},\end{array}\right.$$

$${\tilde{\mathcal{D}}}_1=\underset{\epsilon \to 0}{\lim }{\mathcal{D}}_{1\epsilon },{\tilde{\mathcal{D}}}_3=\underset{\epsilon \to 0}{\lim }{\mathcal{D}}_{3\epsilon },{\tilde{\Sigma }}_{13}=\underset{\epsilon \to 0}{\lim }{\Sigma }_{12\epsilon }=\underset{\epsilon \to 0}{\lim }{\Sigma }_{23\epsilon },$$

$${\tilde{\rho }}_i=\underset{\epsilon \to 0}{\lim }{\rho }_{i\epsilon },{\tilde{\mu }}_i=\underset{\epsilon \to 0}{\lim }{\mu }_{i\epsilon }\text{and}{\tilde{\Gamma }}_i({\tilde{\rho }}_i)=\underset{\epsilon \to 0}{\lim }{\Gamma }_{i\epsilon }({\rho }_{i\epsilon }).$$

$${\mu }_{2\epsilon }\underset{\epsilon \to 0}{\to }0\text{ in such a way that }\frac{{\mu }_{2\epsilon }}{\epsilon }\underset{\epsilon \to 0}{\to }{\tilde{\mu }}_{13},\text{mbox}with{\tilde{\mu }}_{13}:{\mathbb{R}}^{+}\times {\mathcal{D}}_{2\epsilon }\to {\mathbb{R}}^{+}$$

$$\underset{\epsilon \to 0}{\lim }\frac{\nabla {\mu }_{2\epsilon }}{\epsilon }\cdot {\mathbf{n}}_{12}=\underset{\epsilon \to 0}{\lim }\frac{\nabla {\mu }_{2\epsilon }}{\epsilon }\cdot {\mathbf{n}}_{23}=\nabla {\tilde{\mu }}_{13}\cdot {\tilde{\mathbf{n}}}_{13}=0,$$

$$Q=-\frac{\kappa }{\nu }\frac{\Delta P}{\Delta x},$$

$$\kappa \underset{\Delta x\to 0}{\to }0\text{ in such a way that }\frac{\kappa }{\Delta x}\underset{\Delta x\to 0}{\to }\tilde{\kappa },\text{mbox}with\tilde{\kappa }\in {\mathbb{R}}^{+},$$

$${\tilde{\mu }}_1{\tilde{\rho }}_1{f}^{\prime }({\tilde{\rho }}_1)\nabla {\tilde{\rho }}_1\cdot {\tilde{\mathbf{n}}}_{13}={\tilde{\mu }}_3{\tilde{\rho }}_3{f}^{\prime }({\tilde{\rho }}_3)\nabla {\tilde{\rho }}_3\cdot {\tilde{\mathbf{n}}}_{13}\text{on }{\tilde{\Sigma }}_{13}.$$

$${\tilde{\mu }}_{13}[[\Pi ]]={\tilde{\mu }}_1{\tilde{\rho }}_1{f}^{\prime }({\tilde{\rho }}_1)\nabla {\tilde{\rho }}_1\cdot {\tilde{\mathbf{n}}}_{13}={\tilde{\mu }}_3{\tilde{\rho }}_3{f}^{\prime }({\tilde{\rho }}_3)\nabla {\tilde{\rho }}_3\cdot {\tilde{\mathbf{n}}}_{13}\text{on }{\tilde{\Sigma }}_{13},$$

$${\Pi }^{\prime }(\rho ):=\rho f^{\prime }(\rho ).$$

$$\frac{\partial {\rho }_{2\epsilon }}{\partial t}-{\nabla }_{{\mathbf{x}}_{\Sigma }}\cdot ({\mu }_{2\epsilon }{\rho }_{2\epsilon }{f}^{\prime }({\rho }_{2\epsilon }){\nabla }_{{\mathbf{x}}_{\Sigma }}{\rho }_{2\epsilon })-\frac{1}{\epsilon }\frac{\partial }{\partial \eta }\left(\frac{{\mu }_{2\epsilon }}{\epsilon }{\rho }_{2\epsilon }{f}^{\prime }({\rho }_{2\epsilon })\frac{\partial {\rho }_{2\epsilon }}{\partial \eta }\right)={\Gamma }_{2\epsilon }({\rho }_{2\epsilon })$$

$$\displaystyle{\dfrac{\mu_{2\varepsilon}}{\varepsilon}\rho_{2\varepsilon}f^{\prime}(\rho_{2\varepsilon})\frac{\partial\rho_{2\varepsilon}}{\partial\eta}\Big{|}_{\eta=0}=\mu_{1\varepsilon}\rho_{1\varepsilon}f^{\prime}(\rho_{1\varepsilon})\nabla\rho_{1\varepsilon}\cdot\mathbf{n}_{12}\big{|}_{\Sigma_{12\varepsilon}}\,},$$

$$\displaystyle{\dfrac{\mu_{2\varepsilon}}{\varepsilon}\rho_{2\varepsilon}f^{\prime}(\rho_{2\varepsilon})\frac{\partial\rho_{2\varepsilon}}{\partial\eta}\Big{|}_{\eta=1}=\mu_{3\varepsilon}\rho_{3\varepsilon}f^{\prime}(\rho_{3\varepsilon})\nabla\rho_{3\varepsilon}\cdot\mathbf{n}_{23}\big{|}_{\Sigma_{23\varepsilon}}\,}.$$

$$\frac{\partial }{\partial \eta }\left(\frac{{\mu }_{2\epsilon }}{\epsilon }{\rho }_{2\epsilon }{f}^{\prime }({\rho }_{2\epsilon })\frac{\partial {\rho }_{2\epsilon }}{\partial \eta }\right)=\epsilon \left(\frac{\partial {\rho }_{2\epsilon }}{\partial t}-{\nabla }_{{\mathbf{x}}_{\Sigma }}\cdot ({\mu }_{2\epsilon }{\rho }_{2\epsilon }{f}^{\prime }({\rho }_{2\epsilon }){\nabla }_{{\mathbf{x}}_{\Sigma }}{\rho }_{2\epsilon })-{\Gamma }_{2\epsilon }({\rho }_{2\epsilon })\right).$$

$$\rho_{2\varepsilon}\left(\dfrac{x_{\perp}}{\varepsilon},\mathbf{x}_{\Sigma}\right)=\rho_{2}^{0}(\eta,\mathbf{x}_{\Sigma})+\varepsilon\rho_{2}^{1}(\eta,\mathbf{x}_{\Sigma})+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(\varepsilon)$$

$$f^{\prime}(\rho_{2\varepsilon})=f^{\prime}(\rho_{2}^{0})+\varepsilon f^{\prime\prime}(\rho_{2}^{0})\rho_{2}^{1}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(\varepsilon),\qquad\Gamma(\rho_{2\varepsilon})=\Gamma(\rho_{2}^{0})+\varepsilon\Gamma^{\prime}(\rho_{2}^{0})\rho_{2}^{1}+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}(\varepsilon).$$

$$\frac{\partial }{\partial \eta }\left({\tilde{\mu }}_{13}{\rho }_2^0{f}^{\prime }({\rho }_2^0)\frac{\partial {\rho }_2^0}{\partial \eta }\right)=0⟹{\tilde{\mu }}_{13}{\rho }_2^0{f}^{\prime }({\rho }_2^0)\frac{\partial {\rho }_2^0}{\partial \eta }=const.\forall \eta \in (0,1).$$

$$\displaystyle{\tilde{\mu}_{13}\,\rho^{0}_{2}\,f^{\prime}(\rho^{0}_{2})\,\frac{\partial\rho^{0}_{2}}{\partial\eta}\Big{|}_{\eta=0}=\tilde{\mu}_{1}\tilde{\rho}_{1}f^{\prime}(\tilde{\rho}_{1})\nabla\tilde{\rho}_{1}\cdot\tilde{\mathbf{n}}_{13}\big{|}_{\tilde{\Sigma}_{13}}\,},$$

$$\displaystyle{\tilde{\mu}_{13}\,\rho^{0}_{2}\,f^{\prime}(\rho^{0}_{2})\,\frac{\partial\rho^{0}_{2}}{\partial\eta}\Big{|}_{\eta=1}=\tilde{\mu}_{3}\tilde{\rho}_{3}f^{\prime}(\tilde{\rho}_{3})\nabla\tilde{\rho}_{3}\cdot\tilde{\mathbf{n}}_{13}\big{|}_{\tilde{\Sigma}_{13}}\,}.$$

$$\tilde{\mu}_{13}\,\rho^{0}_{2}\,f^{\prime}(\rho^{0}_{2})\,\frac{\partial\rho^{0}_{2}}{\partial\eta}=\tilde{\mu}_{1}\,\tilde{\rho}_{1}\,f^{\prime}(\tilde{\rho}_{1})\,\nabla\tilde{\rho}_{1}\cdot\tilde{{\bf n}}_{13}\big{|}_{\tilde{\Sigma}_{13}}=\tilde{\mu}_{3}\,\tilde{\rho}_{3}\,f^{\prime}(\tilde{\rho}_{3})\,\nabla\tilde{\rho}_{3}\cdot\tilde{{\bf n}}_{13}\big{|}_{\tilde{\Sigma}_{13}}.$$

Full text

1 Introduction

Biological background

Cell migration is crucial to maintain normal homeostasis [58] and sustain many physiological and pathological processes [32, 52, 61, 63]. During migration phenomena, cells encounter a variety of barriers encompassing other cells, cell-cell junctions, and extracellular matrices (ECMs) of different densities and compositions [52].

One of the most difficult barriers for the cells to cross is the basement membrane. This is a thin, dense and highly cross-linked sheet-like network of ECM macromolecules that underlies, among others, all epithelial and endothelial layers [50, 52]. With its pore size being in the order of 50 nm, only small molecules such as nutrients (e.g. oxygen and glucose) and other chemical factors are able to passively diffuse across the basement membrane [50, 76]. Nonetheless, such a structural barrier is crossed daily by billions of cells in healthy tissues in the course of normal immune cell trafficking [45], epithelial-to-mesenchymal transition [79], collective cell migration [32, 61, 63], and tissue development and morphogenesis [80]. Recent empirical studies [52, 83] suggest that during these physiological processes cells can invade the basement membrane and other thin ECM barriers in a variety of ways, including either active removal (e.g. through invadopodia breaching and barrier disruption mediated by the down-regulation of adhesion receptors) or structural remodelling leading to the creation of gaps in the barrier, or even physiological enlargement of preexisting openings that facilitate, for instance, leukocyte trafficking in the vasculature [67].

Similar mechanisms of cell invasion are likely to be activated in pathological conditions, including fibrotic diseases (most commonly affecting the lungs or kidneys), inflammatory diseases, arteriosclerosis and neoplastic processes [52]. In particular, many types of tumours originate and develop in body regions that are separated from the surrounding environment by the basement membrane. This is, for instance, the case of breast tumours (ductal carcinoma) [24], ovary tumours [3], and exocrine or endocrine pancreatic tumours [17]. During the first stages of cancer progression, non-invasive dysplastic cells proliferate locally and form a carcinoma in situ. At some later stage of tumour development, such a localised cancer lesion may acquire the capacity to invade the adjacent tissues by perforating the basement membrane, thus becoming an invasive carcinoma [23, 77]. The transition from carcinoma in situ to invasive carcinoma is sustained by the ability of cancer cells to produce matrix metalloproteinases (MMPs). These are enzymes capable of digesting the collagen fibers that constitute the extracellular environment and the basement membrane [47, 83]. The MMPs’ action widens the pores of the fibre networks and enable cancer cells to spread from the primary site to the surrounding tissues. Notably, experimental studies on cancer cell mobility in MMP-degradable collagen lattices and non-degradable substrates of various porosity have revealed the existence of an ECM critical pore size below which cancer cell migration is entirely hampered in the absence of MMP secretion. Such a critical pore size was termed “the physical limit of migration” [83].

Mathematical modelling background

Despite our growing knowledge about the underpinnings of cell invasion during physiological and pathological processes [32, 44, 52, 76, 83, 84], a number of key aspects still remain unclear. This is mainly due to the difficulty of examining in vivo the interactions occurring between cells and the basement membrane or other ECM barriers during cellular invasion, as well as to the wide range of diverse mechanisms that cells can use to cross different extracellular structures [52]. As a consequence of our partial understanding of this complex biological phenomenon, there has been little prior work on the mathematical modelling of cell invasion through thin membranes. In fact, classical mathematical models of tumour growth [9, 31, 66, 73] and cell migration on two-dimensional flat substrates [26] do not take into account the effect of cell invasion through ECM barriers nor the transition from carcinomas in situ to invasive tumours.

Only more recently physiological and pathological processes involving the migration of single cells in the presence of obstacles or barriers have been mathematically described by means of discrete models [39, 48, 62], and different aspects of tumour growth in confined environments have been investigated in silico using both discrete and hybrid models [41, 53, 54]. These models can be easily tailored to capture fine details of the changes in cell-cell and cell-ECM adhesion properties observed during cell migration. However, their computational cost can become prohibitive for large cell numbers. Therefore, to model cell migration through the basement membrane and other thin ECM barriers at the scale of larger portions of tissues, it is desirable to use continuum models, which offer the possibility to carry out efficient numerical simulations for large cell numbers that are biologically and clinically relevant.

In this regard, focussing on breast cancer, which originates in the epithelial lining of the milk ducts, Ribba et al. [75] have proposed a mathematical model of cancer cell invasion whereby the basement membrane of the ducts is explicitly represented as a weakly permeable thin region. Although it has provided some interesting biological insights, such a modelling approach could become computationally inefficient in the presence of multiple thin membranes, as they would still be modelled as finite regions. Moreover, Gallinato et al. [35] have proposed a mixture model of breast cancer cell invasion whereby the presence of the basement membrane of the milk ducts is taken into account by imposing nonlinear Kedem-Katchalsky interface conditions [19, 28, 29, 51, 55, 70] at the interface between the tumour and the host region. In the setting of Gallinato et al. [35], such transmission conditions lead the normal velocity of the cells and the cell volume fraction to be continuous across the basement membrane, which is not necessarily the case. Finally, Arduino & Preziosi [11] and Giverso et al. [37] have presented a number of multiphase models of cancer cell migration and invasion through the ECM. In agreement with the biological experiments of Wolf et al. [83], in these models the cellular mobility vanishes when the ECM pore size decreases below a certain critical value. These models effectively capture the fact that the ECM critical pore size is relative to the geometrical and mechanical characteristics of the cells (e.g. the size and elasticity of the nucleus, the stiffness of the nuclear membrane, cellular adhesion and traction), and they have been proven useful to study cancer cell invasion in cases where the morphological characteristics of the ECM are spatially heterogeneous, or even discontinuous. However, such models do not apply to biological scenarios where ECM regions with different mechanical and structural properties (i.e. different cell mobilities) are separated by thin membranes.

Contents of the paper

In this paper, we consider a continuum mechanical model of cell movement and proliferation in a spatial domain that is divided into subdomains by one or multiple thin membranes. The model is formulated in terms of a transmission problem defined by a system of nonlinear partial differential equations for the cell volume fraction complemented with mass-continuity and stress-continuity conditions on the interfaces between the membranes and the rest of the domain.

Nonlinear partial differential equations describing reaction-diffusion processes and transport phenomena in spatial domains that comprise different parts separated by thin layers (i.e. films or membranes) arise in the mathematical modelling of various chemical, physical and biological systems [1, 2, 4, 8, 13, 14, 15, 16, 20, 21, 27, 34, 36, 43, 49, 57, 60, 64, 65, 68, 71, 72]. Due to the analytical and numerical challenges posed by the presence of such layers [12], it is often convenient to approximate the original problem by an equivalent transmission problem whereby each thin layer is replaced by an effective interface. The equivalent problem is then closed by imposing appropriate transmission conditions on the effective interfaces.

In this spirit, we develop a formal procedure to derive a set of biophysically consistent interface conditions to close the limiting problem. Specifically, we find that the mass flux across the effective interfaces must be continuous, as one would expect, and proportional to the jump of a term linked to the cell pressure. The biophysical interest lies in the fact that this proportionality coefficient can be related to the size of the pores of the thin membrane, as well as to the geometrical and mechanical characteristics of the cells as in [11, 37, 39]. This makes the limiting transmission problem suitable for providing a possible macroscopic description of cell invasion through thin membranes that takes explicitly into account cell microscopic characteristics, such as the mechanical constraints imposed by the cell nuclear envelope and the solid material surrounded by it [83].

The transmission condition identified by the limiting procedure can be regarded as a nonlinear generalisation of the classical Kedem-Katchalsky interface condition, as it reduces to it for a peculiar (logarithmic) choice of the constitutive relation between the cell pressure and the cell volume fraction. In contrast to other nonlinear Kedem-Katchalsky interface conditions that have been previously employed to model cell invasion through the basement membrane [35], our transmission condition allows the cell volume fraction to be discontinuous across the equivalent interface, while ensuring mass conservation.

The remainder of the paper is organised as follows. In Section 2, we present the original transmission problem and we introduce the related limiting problem. In Section 3, we formally derive a set of effective interface conditions to close the limiting problem. In Section 4, we present sample numerical solutions that illustrate the formal results established in Section 3 and show their potential applications. In particular, we use the limiting transmission problem to describe cancer cell invasion through the basement membrane and to model the metastatic spread of ovarian carcinoma. Section 5 concludes the paper and provides a brief overview of possible research perspectives.

2 Statement of the problem

We consider a population of cells moving through a region of space that is filled with a porous embedding medium, e.g. the ECM. Mathematically, we identify such a region with a simply-connected spatial domain ${\cal D}\subset\mathbb{R}^{d}$, with $d=1,2,3$. Focussing on the biological scenario where the spatial domain is divided into two regions separated by a porous membrane, we let the domain ${\cal D}$ consist of three subdomains represented as the open sets ${\cal D}_{1}$, ${\cal D}_{2}$ and ${\cal D}_{3}$, as in the scheme depicted in Fig. 1 for a three-dimensional case. The subdomain ${\cal D}_{2}$ represents the porous membrane, and the interfaces between the membrane and the subdomains ${\cal D}_{1}$ and ${\cal D}_{3}$ are denoted by $\Sigma_{12}$ and $\Sigma_{23}$, respectively.

We model the cell volume fraction at position ${\bf x}\in{\cal D}$ and time $t\geq 0$ by means of the function $\rho(t,{\bf x})\geq 0$. The evolution of the cell volume fraction is governed by the mass balance equation

[TABLE]

complemented with the momentum-related equation for an elastic fluid, neglecting inertia,

[TABLE]

and a barotropic relation $p\equiv p(\rho)$ for the cell pressure $p$. If necessary, one can let the net growth rate $\Gamma$ depend also on the concentrations of some chemical factors, such as nutrients and growth factors, and couple Eq. (2.1) with the mass balance equations modelling their evolution. In analogy with the classical Darcy’s law for fluids, the function $\mu(t,{\bf x})\geq 0$ is the cell mobility coefficient and Eq. (2.2) models the tendency of cells to move towards regions where they feel less compressed [7].

Remark 2.1**.**

We remark that Eq. (2.2) is only an approximate representation of the far more complex process underlying the migration of cellular aggregates, which is governed by a multitude of sub-cellular pathways involving different proteins and chemical species [82, 81] and is influenced by the mechanical properties both of the single cells and of the sub-cellular elements of the aggregate [37], as well as by the conditions of the surrounding environment. However, when looking at cell migration at the tissue scale, the ensemble of cells that constitute a cellular aggregate can be described as a single phase material – or possibly a multi-phase material – with liquid [59, 18, 25, 33, 42, 22, 38], or elastic/hyperelastic [6, 10, 46], or elasto-viscoplastic [40] characteristics. In particular, the use of a liquid-like constitutive assumption is supported by experimental evidence [30, 5, 78, 74] indicating that cellular aggregates behave like elastic solids over short timescales (i.e. time scales of the order of a few minutes) but eventually display a fluid-like behaviour (i.e. over time scales of the order of cell division and apoptosis). For this reason, the representation of living materials as viscous/inviscid/elastic fluids is now commonly employed [59].

It is important to stress the fact that we let the cell mobility coefficient be a function of both $t$ and ${\bf x}$. This is to take into account the heterogeneous composition of the spatial domain ${\cal D}$ and the biological notion that the mobility of cells in the embedding medium, especially within the membrane, can vary considerably across space and time. Variability of the cell mobility can be due both to local variations in the micro-structure of the ECM and to spatio-temporal changes in the concentration of MMPs. Therefore, one may let the function $\mu$ depend explicitly on the local concentration of MMPs and then couple Eqs. (2.1) and (2.2) with a conservation equation for the MMP concentration, as we will do in Section 4.

From continuum mechanics, one has that mass flux and stresses must be continuous across the interfaces $\Sigma_{12}$ and $\Sigma_{23}$. Within the framework of Eqs. (2.1) and (2.2), such continuity conditions translate into the following interface conditions:

[TABLE]

and

[TABLE]

In Eqs. (2.3) and (2.4), the notation $[\![(\cdot)]\!]$ represents the jump across the interface $\Sigma_{ij}$, i.e. $[\![(\cdot)]\!]:=(\cdot)_{j}-(\cdot)_{i}$, with the subscript $i$ indicating that $(\cdot)$ is evaluated as the limit to a point of the interface coming from the subdomain ${\cal D}_{i}$. Moreover, as shown in Fig. 1, we denote by ${\bf n}_{ij}$ the unit vector normal to the interface $\Sigma_{ij}$ that points towards the subdomain ${\cal D}_{j}$. Substituting the expression (2.2) for the velocity field ${\bf v}$ into the flux-continuity condition (2.3) yields

[TABLE]

In order to close the transmission problem defined by Eqs. (2.1) and (2.2) complemented with the interface conditions (2.4) and (2.5), in addition to prescribing suitable boundary conditions on the outer boundaries (i.e. the non-interfacing parts of the boundaries of the three spatial subdomains) and suitable initial conditions, one should specify a barotropic relation $p(\rho)$.

In general, the three subdomains can differ in their biophysical properties. As a result, the mobility coefficient and the net growth rate can become discontinuous across the interfaces $\Sigma_{12}$ and $\Sigma_{23}$. In this case, denoting by $\rho_{i}(t,{\bf x})$, $\mu_{i}(t,{\bf x})$ and $\Gamma_{i}(\rho_{i})$ the restrictions to the subdomain ${\cal D}_{i}$ of the functions that represent the local cell volume fraction, the mobility coefficient and the net growth rate, respectively, we can rewrite the problem defined by Eqs. (2.1) and (2.2) subject to the interface conditions (2.4) and (2.5) as

[TABLE]

with $j=i+1$. We make the following assumptions:

Assumption 2.1**.**

The cell mobility coefficient $\mu_{i}$ is continuous in both arguments for all $i=1,2,3$.

Assumption 2.2**.**

The net growth rate $\Gamma_{i}$ is a continuously differentiable function of the cell volume fraction for all $i=1,2,3$.

Assumption 2.3**.**

The pressure $p$ is given by a barotropic relation $p\equiv f(\rho)$ where $f$ is a continuously differentiable and monotonically increasing function of the cell volume fraction.

Remark 2.2**.**

In the case where the pressure $p$ is a continuous function of the cell volume fraction $\rho$, the stress-continuity condition (2.4) implies that also $\rho$ is continuous across the interfaces $\Sigma_{12}$ and $\Sigma_{23}$, that is,

[TABLE]

Hence, the flux-continuity conditions (2.3) or (2.5) read as

[TABLE]

In most biologically relevant scenarios arising in the study of cell invasion through the basement membrane and other ECM barriers, the thickness of the membrane or the barrier is much smaller than the characteristic size $L>0$ of the spatial domain. In order to translate this biological observation into mathematical terms, we define the thickness of the membrane represented as the subdomain ${\cal D}_{2}$ as

[TABLE]

and we assume $\varepsilon\ll L$. In the biological scenarios corresponding to the assumption $\varepsilon\ll L$, one typically wishes to:

• i)

replace the subdomain ${\cal D}_{2}$ with an effective interface, which is obtained from the actual interfaces $\Sigma_{12}$ and $\Sigma_{23}$ by letting $\varepsilon\to 0$;

• ii)

find biophysically consistent transmission conditions to impose on the effective interface in this asymptotic regime.

With these goals in mind, we rewrite the transmission problem (2.6) as

[TABLE]

with $j=i+1$, while the limiting transmission problem whereby the subdomain ${\cal D}_{2\varepsilon}$ is replaced by an effective interface reads as

[TABLE]

where

[TABLE]

[TABLE]

Remark 2.3**.**

In the remainder of the paper, we will refer to the transmission problem ${\cal P}_{\varepsilon}$ defined by (2.8), or equivalently by (2.6), as the “thin layer problem”, and to the limiting transmission problem ${\cal P}_{0}$ defined by (2.9) along with the appropriate transmission conditions as the “effective interface problem”.

The next section will be devoted to derive the transmission conditions that are necessary to complete the effective interface problem ${\cal P}_{0}$. For the effective interface $\tilde{\Sigma}_{13}$ (i.e. an infinitesimal region) to have an effect on cell invasion analogous to that of the actual thin membrane represented as the subdomain ${\cal D}_{2\varepsilon}$ (i.e. a finite region), when letting $\varepsilon\to 0$ we will need to compact the membrane (vid. Remark 2.4). In other words, we will obtain the effective interface by virtually shrinking the pores of the membrane in such a way as to cause a reduction in the local permeability $\mu_{2\varepsilon}$ that is proportional to the local shrinkage. This ensures that the existing relationships between the structural characteristics of the thin membrane and the biophysical properties of the cells will remain intact across $\tilde{\Sigma}_{13}$. To this end, we will assume

[TABLE]

and

[TABLE]

where $\tilde{{\bf n}}_{13}$ is the unit vector normal to the interface $\tilde{\Sigma}_{13}$ that points towards the subdomain $\tilde{\cal D}_{3}$.

The positive bounded function $\tilde{\mu}_{13}$ can be seen as the “effective mobility coefficient” of the cells through the thin membrane represented as the effective interface $\tilde{\Sigma}_{13}$.

Remark 2.4**.**

By analogy, consider a liquid flowing through a layer of porous material with unitary cross-sectional area. The liquid flux $Q$ can be computed using the classical Darcy’s law as

[TABLE]

where $\Delta P$ is the pressure drop between the ends of the layer, $\Delta x$ is the thickness of the layer, $\kappa$ represents the hydraulic permeability of the material and $\nu$ is the dynamic viscosity of the liquid. We can draw a conceptual analogy between the biological problem at hand and the case of the liquid by noting that in order to preserve the flux $Q$ when taking the limit $\Delta x\to 0$ the key is to keep the pressure drop $\Delta P$ fixed. This can be achieved by letting

[TABLE]

where $\tilde{\kappa}$ represents an “effective permeability” of the porous layer in the case where the layer is thin. The latter assumption is analogous to assumption (2.12).

3 Formal derivation of the interface conditions for the effective transmission problem

In this section, we formally derive the transmission conditions required to complete the effective transmission problem ${\cal P}_{0}$ defined by (2.9). In summary, as established by Proposition 3.1, we show that the mass flux across the effective interface $\tilde{\Sigma}_{13}$ is continuous and we find an additional transmission condition that establishes a relationship between the mass flux across the effective interface $\tilde{\Sigma}_{13}$ and the effective cell mobility coefficient $\tilde{\mu}_{13}(t,{\bf x})$.

Proposition 3.1**.**

Under Assumptions 2.1-2.3, the following transmission condition formally applies to the effective interface problem (2.9)

[TABLE]

Moreover, under the additional assumptions (2.12) and (2.13),

[TABLE]

where the function $\Pi(\rho)$ is defined according to the equation

[TABLE]

Proof.

For ease of presentation, we formally derive the interface conditions (3.1) and (3.2) in the case where $\Sigma_{12\varepsilon}$ and $\Sigma_{23\varepsilon}$ are parallel planes, but there would be no additional difficulty in considering more general cases. We introduce the notation ${\cal D}_{2\varepsilon}\ni{\bf x}:=(x_{\perp},\mathbf{x}_{\Sigma})$, where $x_{\perp}:={\bf x}\cdot{\bf n}_{12}={\bf x}\cdot{\bf n}_{23}$. We also make the change of variables $x_{\perp}\mapsto x_{\perp}-\hat{x}_{12\perp}$, with $\hat{x}_{12\perp}$ given by $\hat{{\bf x}}_{12}=(\hat{x}_{12\perp},\hat{\bf x}_{12\Sigma})\in\Sigma_{12\varepsilon}$, and let $\eta:=\dfrac{x_{\perp}}{\varepsilon}\in(0,1)$, so that Eq. (2.8) for $\rho_{2\varepsilon}$ can be rewritten as

[TABLE]

and the related flux continuity conditions can be rewritten as

[TABLE]

[TABLE]

Rearranging terms in (3.4) yields

[TABLE]

We make the ansatz

[TABLE]

and compute the asymptotic expansions

[TABLE]

Substituting (3.8) and (3.9) into (3.7), and letting $\varepsilon\to 0$, under assumption (2.12) we formally obtain

[TABLE]

In a similar way, from the flux continuity conditions (3.5) and (3.6) we formally obtain

[TABLE]

[TABLE]

Using (3.10) along with (3.11) and (3.12) we find that for all $\eta\in(0,1)$ we have

[TABLE]

Hence, the transmission condition (3.1) is formally verified. Moreover, under the additional assumption (2.13), integrating both sides of (3.13) with respect to $\eta$ and noting that

[TABLE]

with $\Pi$ defined according to (3.3), we obtain

[TABLE]

Hence, the transmission condition (3.2) is formally verified as well. ∎

Remark 3.1**.**

If the cell pressure is given by the barotropic relation

[TABLE]

then Eq. (2.8) for $\rho_{i\varepsilon}$ becomes a nonlinear reaction-diffusion equation with a nonlinearity only in the reaction term, and $\Pi=P\,\rho+C$ with $C\in\mathbb{R}$. In this case, the interface condition (3.2) reduces to the classical Kedem-Katchalsky interface condition, i.e. $\tilde{\mu}_{13}\,(\tilde{\rho}_{3}-\tilde{\rho}_{1})=\tilde{\mu}_{1}\,\nabla\tilde{\rho}_{1}\cdot\tilde{{\bf n}}_{13}=\tilde{\mu}_{3}\,\nabla\tilde{\rho}_{3}\cdot\tilde{{\bf n}}_{13}$ on $\tilde{\Sigma}_{13}$.

Remark 3.2**.**

If $\tilde{\mu}_{13}\equiv 0$ then the thin membrane represented as the effective interface $\tilde{\Sigma}_{13}$ is impermeable and we recover no-flux boundary conditions on both sides of $\tilde{\Sigma}_{13}$, i.e. the cells in each subdomain are compartmentalised.

Taken together, the formal results established by Proposition 3.1 allow us to complete the effective interface problem ${\cal P}_{0}$ defined by the transmission problem (2.9) as follows

[TABLE]

In order to illustrate these formal results we constructed numerical solutions of a one-dimensional version of the thin layer problem ${\cal P}_{\varepsilon}$ for decreasing values of $\varepsilon$, and we compared the numerical solutions obtained with the numerical solutions of the corresponding effective interface problem ${\cal P}_{0}$. These results are reported in Section S.1 of the Supplementary Material and show that the relative error between the numerical solutions of the two transmission problems tends linearly to zero as $\varepsilon\to 0$.

Remark 3.3**.**

The results established by Proposition 3.1 can also be obtained using a control volume approach analogous to that typically used in continuum mechanics (i.e. considering a control volume that cuts across the subdomain ${\cal D}_{2\varepsilon}$).

4 Application of the effective interface conditions

The numerical solutions presented in this section show potential applications of the formal results established by Proposition 3.1. In Section 4.1, we construct numerical solutions for a two-dimensional model of cancer cell invasion through a basement membrane and the corresponding effective interface problem. The numerical results obtained indicate that the effective interface problem provides a good approximation of the original transmission problem for membranes of sufficiently small thickness. In Section 4.2, we construct numerical solutions for an effective interface problem modelling cell invasion dynamics in ovarian carcinoma. The numerical results obtained support the idea that the model can qualitatively reproduce the key steps of the complex process leading to the metastatic spread of ovarian cancer cells. All numerical simulations are carried out using the finite element software (FEM) COMSOL Multiphysics®, with the parallel sparse direct solver MUMPS. The method for constructing numerical solutions is based on the backward differentiation formula with an adaptive time-step, and a refined mesh is used in the region about the effective interface.

4.1 Numerical simulation of cancer cell invasion through the basement membrane

We compare the numerical solutions of a thin layer problem modelling a two-dimensional cell invasion process with the numerical solutions of the corresponding effective interface problem. We consider a biological scenario whereby cancer cells, which proliferate according to a logistic law with intrinsic growth rate $r>0$, invade a normal tissue composed of healthy cells in homeostatic equilibrium (i.e. cells for which proliferation is balanced by natural death) by squeezing through a damaged part of the basement membrane. Throughout this section, we use the notation ${\bf x}=(x/L,y/L)$ to denote the spatial position non-dimensionalised with respect to the thickness $L>0$ of the region represented as the subdomains ${\cal D}_{1\varepsilon}$ and $\tilde{\cal D}_{1}$, and we non-dimensionalise the time variable with respect to the intrinsic growth rate $r$.

We consider the net growth rate

[TABLE]

where $H(\cdot)$ denotes the Heaviside step function and the function $\varphi(t,{\bf x})$ is an auxiliary level set function that tracks the region of space occupied by cancer cells. Moreover, we use the barotropic relation

[TABLE]

where $(\cdot)_{+}$ is the positive part of $(\cdot)$. We remark that we consider a scenario whereby the cell volume fraction at $t=0$ is equal to or greater than $\rho_{0}$ for all ${\bf x}$. Since $\rho_{0}<1$, under definition (4.1) both the thin layer problem $\mathcal{P}_{\varepsilon}$ and the effective interface problem $\mathcal{P}_{0}$ are such that the cell volume fraction will be greater than or equal to $\rho_{0}$ for all $t\geq 0$. Under this scenario, the barotropic relation (4.2) is such that Assumption 2.3 is satisfied.

We choose the spatial domains schematised in Fig. 2 to carry out numerical simulations. For the thin layer problem [vid. Fig. 2], we let the subdomains ${\cal D}_{1\varepsilon}$ and ${\cal D}_{3\varepsilon}$ be separated by the basement membrane of thickness $\varepsilon$, which is represented as the subdomain ${\cal D}_{2\varepsilon}$ with boundaries $\Sigma_{12\varepsilon}$ and $\Sigma_{23\varepsilon}$. We identify the part of the membrane that is damaged, and thus permeable to cancer cells, with a subset $\mathcal{D}_{p}\subset\overline{{\cal D}}_{2\varepsilon}$. Similarly, for the effective interface problem [vid. Fig. 2], we let the subdomains $\tilde{\cal D}_{1}$ and $\tilde{\cal D}_{3}$ be separated by the effective interface $\tilde{\Sigma}_{13}$. In this case, the damaged part of the basement membrane is represented as a set $\tilde{\Sigma}_{p}\subset\tilde{\Sigma}_{13}$. For simplicity, we assume the cell mobility coefficients in the subdomains ${\cal D}_{1\varepsilon}$ and ${\cal D}_{3\varepsilon}$ to have the same constant value, i.e. $\mu_{1\varepsilon}~{}=~{}\mu_{3\varepsilon}~{}\equiv~{}\bar{\mu}\quad\text{with}\quad\bar{\mu}>0,$ and we define the mobility coefficient in the subdomain ${\cal D}_{2\varepsilon}$ as $\mu_{2\varepsilon}({\bf x})~{}:=~{}\varepsilon\,\bar{\mu}_{2}\,\mathbf{1}_{\mathcal{D}_{p}}({\bf x})\quad\text{with}\quad\bar{\mu}_{2}>0$ where $\mathbf{1}_{\mathcal{D}_{p}}({\bf x})$ is a mollification of the indicator function of the set $\mathcal{D}_{p}\subset\overline{{\cal D}}_{2\varepsilon}$. Accordingly, for the effective interface problem, we choose $\tilde{\mu}_{1}~{}=~{}\tilde{\mu}_{3}\equiv\bar{\mu},\quad\tilde{\mu}_{13}:=\bar{\mu}_{2}\,\mathbf{1}_{\tilde{\Sigma}_{p}}({\bf x})\,.$ We assume that cancer cells initially occupy only the region of space on the left of the membrane, while healthy cells reside in the remaining part of the spatial domain.

We describe the spatio-temporal evolution of the cell volume fraction $\rho_{i\varepsilon}(t,{\bf x})$ through the thin layer problem (2.8) with $f(\rho_{i\varepsilon})$ defined according to (4.2) and $\Gamma_{i\varepsilon}\equiv\Gamma(\varphi_{\varepsilon},\rho_{i\varepsilon})$ given by (4.1). The function $\varphi_{\varepsilon}(t,{\bf x})$ is the auxiliary level set function that tracks the region of space occupied by cancer cells – i.e. at any time instant $t\geq 0$, if $\varphi_{\varepsilon}(t,{\bf x})>0$ then the point ${\bf x}$ is occupied by cancer cells, whereas if $\varphi_{\varepsilon}(t,{\bf x})\leq 0$ then the point ${\bf x}$ is occupied by healthy cells. Hence, the zero level set of the function $\varphi_{\varepsilon}(t,{\bf x})$ corresponds to the boundary of the tumour region at time $t$. The evolution of the function $\varphi_{\varepsilon}(t,{\bf x})$ is governed by the following equation [69]

[TABLE]

for $i=1,2,3$, subject to the continuity conditions

[TABLE]

Notice that the transmission conditions (2.8)2 ensure the continuity of the normal velocity across the interfaces $\Sigma_{12\varepsilon}$ and $\Sigma_{23\varepsilon}$.

The corresponding effective interface problem is given by the transmission problem (3.14) with $\tilde{\Gamma}_{i}\equiv\Gamma(\tilde{\varphi},\tilde{\rho}_{i})$ defined according to (4.1) and with $f(\tilde{\rho}_{i})$ given by (4.2). As for the thin layer problem, the function $\tilde{\varphi}(t,{\bf x})$ is the level set function tracking the region of space occupied by cancer cells, the evolution of which is governed by the following equation [69]

[TABLE]

subject to the continuity condition

[TABLE]

A formal derivation of condition (4.6) is provided in Section S.2 of the Supplementary Material. Finally, we choose parameter values, boundary conditions and initial conditions corresponding to those of the thin layer problem.

The numerical results obtained are summarised by the plots in Figs. 3 and 4. The plots on the top line of Fig. 3 display the numerical solutions to the thin layer problem with $\varepsilon=0.1$ at different time instants. The numerical solutions to the effective interface problem at the same time instants are displayed in the plots on the bottom line. The discrepancy between the solutions to the thin layer problem and the solutions to the effective interface problem decays over time as the invasion front of cancer cells moves away from the basement membrane, which is represented either by the subdomain ${\cal D}_{2\varepsilon}$ or by the effective interface $\tilde{\Sigma}_{13}$. This is further clarified by the plots in Fig. 4. In particular, the curves reported in Fig. 4 indicate that the relative error between the numerical solution to the thin layer problem at the point $(\varepsilon,0)$ and the numerical solution to the effective interface problem at the point $(0^{+},0)$ decays over time. Moreover, in agreement with the formal results established by Proposition 3.1, the relative error decays as $\varepsilon\to 0$.

4.2 Numerical simulation of ovarian cancer invasion

In this section, we apply the formal results established by Proposition 3.1 to the mathematical modelling of cell invasion dynamics in ovarian carcinoma. In particular, we simulate the metastatic journey of a cancer multicellular mass, from the initial growth inside the ovary to the invasion of the healthy tissue adjacent to the peritoneum, using an effective interface problem.

For the sake of brevity, throughout this section we drop the tildes from all quantities and we work with dimensionless quantities, as specified in the previous subsection. In particular, we use the notation ${\bf x}=(x/L,y/L)$ to denote the spatial position non-dimensionalised with respect to the characteristic size $L>0$ of the region represented as the subdomain ${\cal D}_{1}$.

4.2.1 Biological background

Ovarian carcinoma originates either inside the ovary or in the fallopian tube. This type of cancer is known to invade the surrounding tissues and to metastasise both by direct extension and by cell detachment from the primary tumour [56]. The latter process of metastasis formation is peculiar to ovarian carcinoma and allows cancer cells to spread into the peritoneal cavity, to invade adjacent peritoneal tissues and, ultimately, to reach distant organs. Such a process encompasses multiple layers of complexity, which represents one of the main reasons why the metastatic behaviour of ovarian cancer cells remains poorly understood.

The detachment of ovarian cancer cells from the primary tumour starts with the destruction of the basement membrane underling the ovarian capsule (i.e. the ovarian surface epithelium) [3]. Cancer cells can subsequently break through the ovarian capsule as single cells or, more frequently, as spheroid-like aggregates. Such multicellular masses grow and passively move until they reach the walls of the peritoneal cavity – which represent the common site of disaggregation, dissemination and metastatic outgrowth for ovarian carcinoma [41].

The cancer cells that reach the walls of the cavity can attach to the mesothelial cells that constitute the peritoneal lining and, by secreting MMPs [56], they can degrade the basement membrane underling the mesothelium and cleave cell-cell adhesion molecules (e.g. N-cadherins) that hold mesothelial cells together [56]. This leads to the retraction of mesothelial cells at the cancer cells’ attachment sites and promotes the formation of foci of invasion, which enable the ovarian cancer cells to invade the healthy tissue adjacent to the peritoneum and form secondary tumours [41].

4.2.2 Mathematical model

In adult human females, the ovarian capsule consists of a single layer of epithelial cells and the peritoneal lining is constituted by a monolayer of mesothelial cells [3]. Hence, the thickness of the ovarian capsule and the peritoneal lining is small compared to the characteristic size of the ovary and of the peritoneal cavity. For this reason, we represent both the ovarian capsule and the peritoneal lining, along with the underling basement membranes, as two thin porous membranes. Moreover, using the formal results established by Proposition 3.1, we model each thin porous membrane as an effective interface.

On the basis of these observations, considering a two-dimensional scenario, we represent the ovary, the peritoneal cavity and the healthy tissue adjacent to the peritoneum as three distinct spatial subdomains ${\cal D}_{1}$, ${\cal D}_{2}$ and ${\cal D}_{3}$ separated by the effective interfaces $\Sigma_{12}$ (i.e. the ovarian capsule along with the underlying basement membrane) and $\Sigma_{23}$ (i.e. the peritoneal lining along with the underlying basement membrane) – cf. respectively, the blue curve and the red line in Fig. 5. We focus on the biological scenario whereby there is a part of the ovarian capsule that is damaged and thus permeable to cancer cells. We identify such a region with a subset $\Sigma_{p}$ of the effective interface $\Sigma_{12}$ (cf. the green line in Fig. 5).

Letting the function $\rho_{i}(t,{\bf x})$ model the cell volume fraction at position ${\bf x}\in{\cal D}_{i}$ and time $t\geq 0$, we describe the spatio-temporal evolution of the cells through the effective interface problem (3.14), with $i=\left\{1,2,3\right\}$, posed on the spatial domain illustrated in Fig. 5. Similarly to Section 4.1, we define $f(\rho_{i})$ according to (4.2) and we let ovarian cancer cells proliferate in all subdomains ${\cal D}_{i}$, with $i=\left\{1,2,3\right\}$, according to the net growth rate $\Gamma_{i}\equiv\Gamma(\varphi,\rho_{i})$ given by the logistic law (4.1). The evolution of the function $\varphi(t,{\bf x})$ is governed by Eq. (4.5) posed on the spatial domain illustrated in Fig. 5 and subject to the continuity condition (4.6) on $\Sigma_{12}$ and $\Sigma_{23}$.

We make the prima facie assumption that the effective mobility coefficient ${\mu}_{12}$ is a given function of ${\bf x}$ and does not depend on $t$. On the other hand, on the basis of the biological facts discussed in Section 4.2.1, we let the effective mobility coefficient ${\mu}_{23}$ be a function of the local concentration of MMPs $c(t,{\bf x})$, which can vary across space and time. In particular, using a modelling strategy similar to that proposed by Gallinato et al. [35] and Giverso et al. [37], we define ${\mu}_{23}$ as

[TABLE]

A detailed derivation of definition (S.3.3) is provided in Section S.3 of the Supplementary Material. Denoting the restriction of the function $c$ to the subdomain ${\cal D}_{i}$ by $c_{i}$, we describe the dynamics of the concentration of MMPs through the following transmission problem

[TABLE]

where the parameter $\gamma_{c}>0$ is the rate at which cancer cells release MMPs and the parameter $D_{c}>0$ is the diffusivity of MMPs. Notice that the transmission conditions in (4.8) are such that the MMP concentration $c(t,{\bf x})$ and its flux are continuous across the effective interfaces $\Sigma_{12}$ and $\Sigma_{23}$. This is because the size of the MMP molecules is much smaller than the size of the pores of the membranes (i.e. the membranes are permeable to the MMP molecules). Alternatively, one could impose the classical Kedem-Katchalsky interface conditions on $\Sigma_{12}$ and $\Sigma_{23}$.

4.2.3 Numerical solutions

The numerical results obtained are summarised by the plots in Fig. 6. As illustrated by these plots, which display the cell volume fraction in the different subdomains along with the boundaries of the cancer multicellular mass (white lines), the mathematical model defined by the effective interface problem (3.14) posed on the spatial domain of Fig. 5 and coupled with the transmission problem (4.8) can qualitatively reproduce the salient steps of the metastatic journey undertaken by an ovarian cancer multicellular mass.

In summary, cancer cells are initially confined to the ovary region ${\cal D}_{1}$ [vid. Fig. 6(a)], where they proliferate and grow into a multicellular mass. At later stages [vid. Figs. 6(b)-6(d)], cancer cells break through the damaged part of the ovarian capsule $\Sigma_{p}\subset\Sigma_{12}$ and spread across the peritoneal region ${\cal D}_{2}$, until they reach the peritoneal lining $\Sigma_{23}$. From there, secreting MMPs, cancer cells create one focus of invasion [vid. Fig. 6(e)], which enables the multicellular mass to squeeze through the peritoneal lining and form a secondary tumour in the healthy tissue adjacent to the peritoneum ${\cal D}_{3}$ – vid. Figs. 6(f)-6(h).

Note that the plots in Figs. 6(e)-6(h) indicate that the size of the focus of invasion grows over time. This is due to the diffusion of MMPs secreted by cancer cells, which increase the local value of the effective mobility coefficient $\mu_{23}(t,{\bf x})$ [cf. the expression given by Eq. (S.3.3)]. Moreover, throughout the simulations one can verify that the cell volume fractions can become discontinuous not only in the portions of the ovarian capsule and of the peritoneal lining that are impermeable, but also in the permeable part of the ovarian capsule and at the focus of invasion in the peritoneal lining.

5 Conclusions and research perspectives

We have developed a formal asymptotic method to mathematically address biological problems of cell invasion through thin membranes (i.e. the basement membrane and other ECM barriers of small thickness). We have showed how, starting from an original transmission problem in which thin membranes are represented as finite regions of small thickness, one can obtain a limiting transmission problem where each membrane is replaced by an effective interface, and we derived a set of biophysically consistent interface conditions to close the limiting problem.

The approximation of a thin porous layer with an effective interface and a set of suitable transmission conditions is a simplifying approach that has attracted attention in a wide range of application fields – e.g. heat transfer problems [68], flow simulations in porous media with immersed intersecting fractures [14], structural mechanical problems [15, 13], flows through thin membranes for biological applications [55] – as it brings considerable modelling and computational benefits. From the modelling point of view, the main benefit lies in the fact that, by using this approximation, one does not need to develop a detailed model of the phenomena that occur inside the thin layer. From the computational point of view, such an approximation ensures a stark reduction of simulation time in the case of very thin layers, since it makes it possible to avoid the computational cost associated with the fine mesh required to produce accurate numerical results in the proximity of a thin layer, where sharp variations of the dependent variables can lead to the emergence of numerical instabilities. The price to pay for having a simpler and more computationally efficient model is the introduction of effective interface parameters, such as our “effective mobility coefficient”, the estimation of which may require ad hoc experiments and extensive parameter fitting.

The formal results obtained have been validated via numerical simulations showing that the relative error between the solutions to the original transmission problem and the solutions to the limiting problem vanishes when the thickness of the membranes tends to zero. Moreover, in order to show potential applications of our effective interface conditions, we have employed the limiting transmission problem to model cancer cell invasion through the basement membrane and the metastatic spread of ovarian carcinoma.

Our work can be extended both from the analytical perspective and from the modelling point of view. From the analytical perspective, it would be interesting to provide a rigorous proof of the formal results established by Proposition 3.1. From the modelling point of view, we would like to generalise the results presented in this paper to the case of multiple cell populations. Moreover, it would be interesting to understand how to develop further our formal method for deriving effective interface conditions to consider momentum-related equations different from Eq. (2.2), in order to capture the visco-elasto-plastic behaviour of cellular aggregates, which is induced by the dynamic formation of bonds between cells and by the interaction between cells and the extracellular environment, and the active response of living aggregates.

Although the focus of this work has been on cancer invasion, cell penetration of thin membranes occurs also during development, immune surveillance and disease states other than cancer, such as fibrosis [76]. Hence, the effective interface conditions that we have derived can find fruitful application in a variety of research fields in the biological and medical sciences, including developmental biology and immunology.

Supplementary Material

Appendix S.1 Numerical solutions to a one-dimensional problem illustrating the results of Proposition 3.1

In order to illustrate the formal results established by Proposition 3.1, we construct numerical solutions to a one-dimensional thin layer problem $\mathcal{P}_{\varepsilon}$ of mobility $\mu_{2\varepsilon}=\varepsilon\,\bar{\mu}_{2}$, where $\bar{\mu}_{2}>0$. We compare the solutions obtained for decreasing values of $\varepsilon$ with the numerical solutions of the corresponding effective interface problem ${\cal P}_{0}$ of effective mobility $\tilde{\mu}_{13}=\mu_{2\varepsilon}/\varepsilon\equiv\bar{\mu}_{2}$. Throughout this section we make use of the notation ${\bf x}=x$.

For the solution of the thin layer problem ${\cal P}_{\varepsilon}$ to converge to a stationary profile that is non-constant in $x$, we consider a somehow artificial scenario whereby the cells proliferate according to a logistic law with intrinsic growth rate $r>0$ in the subdomain ${\cal D}_{1\varepsilon}$, whereas cell proliferation is balanced by natural death in the subdomains ${\cal D}_{2\varepsilon}$ and ${\cal D}_{3\varepsilon}$. Under these assumptions, letting $L>0$ be the thickness of the region represented as the subdomain ${\cal D}_{1\varepsilon}$, we introduce the non-dimensionalised independent variables $\hat{t}=r\,t$ and $\hat{x}=x/L$ so that, dropping the carets from the non-dimensionalised quantities, we have

[TABLE]

and

[TABLE]

We assume the cell mobility coefficients in the subdomains ${\cal D}_{1\varepsilon}$ and ${\cal D}_{3\varepsilon}$ to have the same constant value, i.e.

[TABLE]

Moreover, we use the following barotropic relation

[TABLE]

where $(\cdot)_{+}$ is the positive part of $(\cdot)$. We impose zero Neumann boundary condition on the left outer boundary, a Dirichlet boundary condition on the right outer boundary, and the following initial conditions

[TABLE]

These initial conditions model a biological scenario where the cells are initially uniformly distributed in space at the stress-free state.

Similarly, for the effective interface problem ${\cal P}_{0}$ we consider

[TABLE]

[TABLE]

Furthermore, we use the barotropic relation (S.1.1). We impose zero Neumann boundary condition on the left outer boundary, a Dirichlet boundary condition on the right outer boundary, and the following initial conditions

[TABLE]

We remark that we consider a scenario whereby the cell volume fraction at $t=0$ is equal to or greater than $\rho_{0}$ for all ${\bf x}$. Since $\rho_{0}<1$, under the above definitions of the growth rates $\Gamma_{i\varepsilon}$ and $\tilde{\Gamma}_{i}$ both the thin layer problem $\mathcal{P}_{\varepsilon}$ and the effective interface problem $\mathcal{P}_{0}$ are such that the cell volume fraction will be greater than or equal to $\rho_{0}$ for all $t\geq 0$. Under this scenario, the barotropic relation (S.1.1) is such that Assumption 2.3 is satisfied.

To construct numerical solutions we choose

[TABLE]

and we carry out computational simulations for $t\in[0,20]$, since numerical solutions appear to be stationary at $t=20$. The results obtained are summarised by the plots in Fig. 7, which display the numerical solutions to the effective interface problem $\mathcal{P}_{0}$ and the numerical solutions to the thin layer problem $\mathcal{P}_{\varepsilon}$ for decreasing values of $\varepsilon$ at $t=20$.

The curves in Fig. 7(a) indicate that the discrepancy between the numerical solutions to the thin layer problem and the numerical solutions to the effective interface problem decreases as $\varepsilon$ tends to zero. This is more precisely quantified by the curves in Fig. 7(b), which display: the relative error between the numerical solutions to the thin layer problem at $x=0$ and the numerical solutions to the effective interface problem at $x=0^{-}$ for $t=20$ as a function of $\varepsilon$ (blue line); the relative error between the numerical solutions to the thin layer problem at $x=\varepsilon$ and the numerical solutions to the effective interface problem at $x=0^{+}$ for $t=20$ as a function of $\varepsilon$ (red line). In agreement with the formal results established by Proposition 3.1, both relative errors tend to zero (linearly) as $\varepsilon\to 0$. Taken together, the numerical results presented in this section illustrate that there is a good match between the numerical solutions to the original transmission problem with membrane thickness $\varepsilon$ and mobility $\mu_{2\varepsilon}$ and the numerical solutions to the limiting transmission problem with effective mobility $\tilde{\mu}_{13}$. Hence, when the thickness of the membrane represented as the subdomain ${\cal D}_{2\varepsilon}$ is small, instead of solving the problem $\mathcal{P}_{\varepsilon}$ one can solve the approximate problem $\mathcal{P}_{0}$ with effective mobility coefficient $\tilde{\mu}_{13}$. The quality of the approximation is higher for membranes of smaller thickness.

Appendix S.2 Formal derivation of the continuity condition (4.6) for $\tilde{\varphi}$

Proceeding in the same way as at the beginning of the proof of Proposition 3.1, we rewrite Eq. (4.3) for $\varphi_{\varepsilon}$ in ${\cal D}_{2\varepsilon}$ as

[TABLE]

Multiplying both sides of the latter equation by $\rho_{2\varepsilon}$ and rearranging terms yields

[TABLE]

Substituting the ansatz

[TABLE]

and

[TABLE]

into the above equation, letting $\varepsilon\to 0$ and using the fact that, as shown in the proof of Proposition 3.1,

[TABLE]

we formally obtain

[TABLE]

from which we deduce the continuity condition (4.6) for $\tilde{\varphi}$. We remark that $\displaystyle{\tilde{\mu}_{13}\,\rho^{0}_{2}\,f^{\prime}(\rho^{0}_{2})\,\frac{\partial\rho^{0}_{2}}{\partial\eta}}=0$ only if \displaystyle{\tilde{\mu}_{1}\,\tilde{\rho}_{1}\,f^{\prime}(\tilde{\rho}_{1})\,\nabla\tilde{\rho}_{1}\cdot\tilde{{\bf n}}_{13}\big{|}_{\tilde{\Sigma}_{13}}}=\tilde{\mu}_{3}\,\tilde{\rho}_{3}\,f^{\prime}(\tilde{\rho}_{3})\,\nabla\tilde{\rho}_{3}\cdot\tilde{{\bf n}}_{13}\big{|}_{\tilde{\Sigma}_{13}}=0. In this case, the cell flux across ${\tilde{\Sigma}_{13}}$ is identically zero and, therefore, the level set does not cross the effective interface.

Appendix S.3 Derivation of the definition (4.7) of ${\mu}_{23}(t,{\bf x})$

Using a modelling strategy similar to that proposed by Gallinato et al. [35] and Giverso et al. [37], we define the effective mobility coefficient ${\mu}_{23}$ as

[TABLE]

where $\bar{\mu}_{23}>0$ is the maximum mobility of ovarian cancer cells through the interface that models the peritoneal lining, the function $A(t,{\bf x})>0$ represents the average cross-section of the pores of the membrane at position ${\bf x}\in\Sigma_{23}$ and time $t\geq 0$, and the parameter $A_{0}>0$ is the critical value of the average pores’ cross-section below which, according to “the physical limit of cell migration” [83], the membrane is completely impermeable to cancer cells. The evolution of the function $A(t,{\bf x})$ is governed by the following differential equation

[TABLE]

In Eq. (S.3.2), the parameter $A_{1}>0$, with $A_{1}<A_{0}$, models the non-dimensionalised average cross-section of the pores of the membrane in physiological conditions, $\alpha$ is the rate of remodelling of $A$ to the normal value $A_{1}$, and $\beta>0$ represents the rate at which MMPs increase the size of the pores of the membrane. Under the biologically realistic assumption that the process of pore cleavage and repairing is much faster than tumour expansion, we assume the average cross-section of the pores of the membrane to be in quasi-stationary equilibrium and rewrite Eq. (S.3.2) as

[TABLE]

Substituting the above expression for $A(t,{\bf x})$ into (S.3.1) and rearranging terms gives

[TABLE]

with

[TABLE]

In Section 4.2, we use (S.3.3) and, with a slight abuse of notation, we rename the rescaled concentration of MMPs $\hat{c}(t,{\bf x})$ to $c(t,{\bf x})$.