A tractable mathematical model for tissue growth

TL;DR

This paper develops a simplified mathematical model for tissue growth using free boundary problems and mean curvature flow, providing analytical and numerical insights into tumor dynamics.

Contribution

It introduces a tractable free boundary model for tissue growth incorporating a PDE-driven forcing term and offers linear stability analysis and diffuse-interface approximations.

Findings

Model captures key features of tissue growth and death.

Finite-element simulations validate the diffuse-interface approximation.

Stability analysis provides insights into interface evolution.

Abstract

Using formal asymptotic methods we derive a free boundary problem representing one of the simplest mathematical descriptions of the growth and death of a tumour or other biological tissue. The mathematical model takes the form of a closed interface evolving via forced mean curvature flow (together with a `kinetic under-cooling' regularisation) where the forcing depends on the solution of a PDE that holds in the domain enclosed by the interface. We perform linear stability analysis and derive a diffuse-interface approximation of the model. Finite-element discretisations of two closely related models are presented, together with computational results comparing the approximate solutions.

Tables

Table 1. Table 1 : Parameters for the diffuse–interface simulations for γ = 0 𝛾 0 \gamma=0 .

	$\epsilon$	$\Delta t$	$h_{max,f}$
$\gamma =0$, $\alpha =1$, $\beta =0.1$	$0.075$	$1.0\times {10}^{-4}$	0.0048
$\gamma =0$, $\alpha =0.1$, $\beta =0.1$	$0.01$	$1.0\times {10}^{-5}$	0.0028
$\gamma =0$, $\alpha =1$, $\beta =1$	$0.1$	$3.0\times {10}^{-6}$	0.0099
$\gamma =0$, $\alpha =0.1$, $\beta =1$	$0.1$	$1.0\times {10}^{-5}$	0.0042

Table 2. Table 2 : Parameters for the diffuse–interface simulations with β = 0 𝛽 0 \beta=0 .

	$\epsilon$	$\Delta t$	$h_{max,f}$
$\beta =0$, $\alpha =1$, $\gamma =0.1$	$0.04$	$1.0\times {10}^{-4}$	0.0048
$\beta =0$, $\alpha =0.1$, $\gamma =0.1$	$0.01$	$1.0\times {10}^{-5}$	0.0028
$\beta =0$, $\alpha =1$, $\gamma =1$	$0.1$	$3.0\times {10}^{-6}$	0.0099
$\beta =0$, $\alpha =0.1$, $\gamma =1$	$0.1$	$1.0\times {10}^{-5}$	0.0042

Equations

$$\Delta u=1$$

$$\displaystyle\nabla u\cdot{\mbox{\boldmath$\displaystyle\nu$\unboldmath}}+\frac{u}{\alpha}+\beta\kappa=Q$$

$$V=\frac{u}{\alpha }+\beta \kappa$$

$$\nabla u\cdot{\mbox{\boldmath$\displaystyle\nu$\unboldmath}}+V=Q,$$

$$u=\alpha (Q-s(t))+\frac{1}{2}(x^2-s^2(t)),$$

$$\frac{\mathrm{d}s}{\mathrm{d}t}=Q-s,s(t)=Q-(Q-s(0))e^{-t}.$$

$$\text{mbox}(lineargrowth)s(t)\sim s(0)+Qt$$

$$\text{mbox}(transitiontosaturation)s(t)\sim Q(1-e^{-t})$$

$$\left.\begin{array}[]{c}\frac{\partial n}{\partial t}+\nabla\cdot(\bm{v}n)=(k_{b}(c)-k_{d}(c))n,\quad\frac{\partial m}{\partial t}+\nabla\cdot(\bm{v}m)=k_{d}(c)n-\varepsilon\lambda m,\\ n+m=1,\\ \varepsilon^{2}\nabla\cdot(D(n)\nabla c)=K(c)n,\quad\bm{v}=-\nabla p/\mu(n;\varepsilon),\end{array}\right\}$$

$$c=1,\quad p=-\varepsilon^{2}\gamma(n)\kappa,\quad q_{{{\small\nu}}}=\bm{v}\cdot{\mbox{\boldmath$\displaystyle\nu$\unboldmath}},\quad\mbox{ on }\Gamma(t),$$

$$n(x,0)=1\text{mbox}in\Omega (0),$$

$$\mu (n;\epsilon )\sim {\mu }_0(n)\text{mbox}forn=O(1),n>0,\mu (0;\epsilon )\sim \epsilon {\mu }_1\text{mbox}as\epsilon \to 0,$$

$$\nu=\varepsilon\zeta,\quad\bm{v}\cdot{\mbox{\boldmath$\displaystyle\nu$\unboldmath}}=\varepsilon V_{{{\small\nu}}},\quad q_{{{\small\nu}}}=\varepsilon Q_{{{\small\nu}}},\quad p=\varepsilon^{2}P,\quad t=\varepsilon^{-1}T,$$

$$\frac{\mathrm{d}}{\mathrm{d}\zeta }(W_{\nu }n)=(k_b(c)-k_d(c))n,$$

$$n+m=1,$$

$$\frac{\mathrm{d}}{\mathrm{d}\zeta }\left(D(n)\frac{\mathrm{d}c}{\mathrm{d}\zeta }\right)=K(c)n,$$

$$c=1,P=-\gamma (n)\kappa ,W_{\nu }=0,\text{mbox}on\zeta =0,$$

$$\frac{\mathrm{d}{W}_{\nu }}{\mathrm{d}\zeta }=k_b(1)-k_d(1),m\frac{\mathrm{d}{W}_{\nu }}{\mathrm{d}\zeta }=k_d(1)n\text{mbox}at\zeta =0,$$

$$n=1-\frac{k_d(1)}{k_b(1)},m=\frac{k_d(1)}{k_b(1)}\text{mbox}at\zeta =0,$$

$${\gamma }_0=\gamma \left(1-\frac{k_d(1)}{k_b(1)}\right),$$

$$n\to 0,m\to 1,c\to c_{\infty },W_{\nu }\to -Q,$$

$$V_{\nu }\sim -\frac{1}{\mu (n;\epsilon )}\frac{\mathrm{d}P}{\mathrm{d}\zeta },$$

$$P\sim \epsilon {\mu }_1(Q-Q_{\nu })\zeta +Q_{\nu }{\int }_{\zeta }^0\left(\mu (n({\zeta }^{\prime });\epsilon )-\epsilon {\mu }_1\right)\mathrm{d}{\zeta }^{\prime }+{\int }_{\zeta }^0\left(\mu (n({\zeta }^{\prime });\epsilon )W_{\nu }({\zeta }^{\prime })+\epsilon {\mu }_{1}Q\right)\mathrm{d}{\zeta }^{\prime }-{\gamma }_0\kappa$$

$$\frac{\mathrm{d}{W}_{\nu }}{\mathrm{d}\zeta }=k_b(c)n.$$

$$P\sim \epsilon {\mu }_1(Q-Q_{\nu })\zeta +a{Q}_{\nu }-b-{\gamma }_0\kappa ,$$

$$a\equiv {\int }_{-\infty }^0{\mu }_0(n(\zeta ))\mathrm{d}\zeta ,b\equiv -{\int }_{-\infty }^0{\mu }_0(n(\zeta ))W_{\nu }(\zeta )\mathrm{d}\zeta .$$

$$\nabla \cdot \mathbf{V}=-\lambda ,\mathbf{V}=-\nabla P/{\mu }_1,$$

$$\Delta c=0.$$

$$V_{\nu }=Q_{\nu }-Q,P=a{Q}_{\nu }-b-{\gamma }_0\kappa ,$$

$${\mbox{\boldmath$\displaystyle\nu$\unboldmath}}\cdot\nabla P=\mu_{1}V_{{{\small\nu}}},$$

$$\Delta P=\lambda {\mu }_1,\text{mbox}in\Omega (t),$$

$$Q_{\nu }=Q-\frac{1}{{\mu }_1}\frac{\partial P}{\partial \nu },P=a{Q}_{\nu }-{\gamma }_0\kappa ,\text{mbox}on\Gamma (t),$$

$${\dot{s}}_c=Q-\frac{1}{2}{s}_c,u_c=\frac{1}{4}{r}^2-\frac{1}{4}{s}_c^2-\frac{\alpha \beta }{s_c}+\alpha {\dot{s}}_c,$$

Full text

1 Introduction

The aim of this paper is to derive and simulate numerically one of the simplest mathematical descriptions of the growth and death of a tumour (or, more generally, other tissue in a multicellular organism or a bacterial biofilm) comprising a nutrient–deprived necrotic (dead) core surrounded by a thin rim of nutrient–rich dividing cells. We emphasise we are not seeking biological realism here, but rather to clarify the properties of simple formulations, onto which any level of additional complexity can be built; we also seek to provide proof of concept for the numerical approaches adopted. Despite the simplicity of the model, the numerical simulations in Section 6 show that the model produces many features that are seen in significantly more complicated ones. Applying formal asymptotic methods, the thin rim collapses onto the interface that separates the tumour from its exterior; this generates two non–standard terms in the moving–boundary conditions for the tumour (the model otherwise corresponding to a Hele–Shaw reverse squeeze film with surface tension). We assume that the nutrient density is constant on the tumour surface and Darcy flow is used, this being the simplest constitutive assumption appropriate to modelling in a porous tissue–engineering scaffold, and having historically been widely adopted (at least since [14, 15]) in the modelling of tumour growth. This results (see below, where we clarify the role of nutrient in this formulation) in the following free boundary problem for the tissue pressure $u$,

[TABLE]

Here $\Omega(t)$ denotes the tumour, $\Gamma(t)=\partial\Omega(t)$, and $V$, $\kappa$ and $\displaystyle\nu$ are respectively the normal velocity, mean curvature and outer unit normal of $\Gamma(t)$, and $\alpha,\beta\in\mathbb{R}$ are positive constants. The key ingredients of the model are as follows: the right–hand side of $\Delta u=1$ is a volumetric sink that reflects loss of necrotic material and destabilises the free boundary; $Q$ is a surface source, capturing the viable rim around the tumour edge; and $\alpha$ and $\beta$ are regularising parameters, the unregularised case being associated with the limit $\alpha\rightarrow 0$ with $\beta$ fixed - the exact implication

[TABLE]

of (1.1) – (1.3) is instructive for this limit case.

The novelty of this free boundary problem is that it is derived via a thin–rim limit, resulting in a formulation that does not require additional equations to be solved for the nutrient density in the interior of the tumour. For a brief discussion of the framework of the model see [17], pp. 305.

The model is both simpler and more tractable than many found in the literature, compare for example [6, 12] and the references therein, and see [5] and [18] for additional background to moving–boundary approaches to tumour growth related to that which is analysed here; indeed, it is perhaps the simplest such model that is able to capture the transition from the linear growth phase to growth saturation: this can be exemplified by the one–dimensional case: imposing symmetry about $x=0$ and writing $\Gamma(t)$ as $x=s(t)$ gives

[TABLE]

with

[TABLE]

The case $s(0)\ll Q$ is then instructive, giving

[TABLE]

Two–component mixture models of tumour and healthy cells are considered in [6, 12], these models taking the form of Cahn–Hilliard type diffuse–interface models in which the interfacial region between the tumour and the healthy cells has a non–zero thickness. In [18, 19] both Darcy and Stokes constitutive assumptions are considered, allowing the tumour to be modelled as a highly viscous fluid moving either through a porous medium or unconstrained; this is also the case in [6], where movement of cells along chemical gradients (chemotaxis) is introduced via a nutrient diffusion equation. In [12] a Cahn–Hilliard–Darcy model is considered that models chemotaxis and active transport, and formal matched asymptotic expansions are performed to yield several sharp–interface models. In [14] the authors consider a cylindrical or spherical tumour and focus on the diffusion of nutrients and waste products in order to model tumour growth patterns. The geometrically constrained model in [14] was extended in [15] whereby, among other things, the geometric constraints are dropped: here the surface normal stress is taken to be proportional to the mean curvature, and Darcy’s law is adopted. In [5] the authors consider a single–nutrient sharp–interface model with two growth inhibitors, one external to the tumour, e.g. an anti–cancer drug or immune response, and one internal to the tumour, e.g. a by–product of the degradation of necrotic cells. The diffusion of the nutrient is considered, and its description on the boundary depends upon the curvature. Although we do not consider anti–cancer drugs or immune responses here, a simplified version of these could potentially be modelled by reducing the size of our surface source.

In Section 2 we derive the model using formal–asymptotic methods then, in Section 3, we present calculations pertaining to the linear stability of radially symmetric solutions. A diffuse–interface formulation is derived in Section 4 and in Section 5 we present finite–element approximations of the model and its diffuse–interface approximation. Numerical simulations are presented and discussed in Section 6. In the appendix we perform analysis on the thin–film limit of the model.

2 Derivation of the model

The model corresponds to a distinguished limit of the following dimensionless formulation (here we elaborate, with some minor differences, on analysis briefly outlined in [19]). We have

[TABLE]

wherein $n$ and $m$ are the volume fractions of live and dead cells, $c$ is the nutrient concentration, $k_{b}$ and $k_{d}$ are the cellular birth and death rates (the former being an increasing function of $c$ and the latter a decreasing one), $\varepsilon\lambda$ specifies the (slow) degradation rate of the necrotic material, $K(c)/\varepsilon^{2}$ expresses the nutrient consumption rate, which will be taken to be large, $D(n)$ is the nutrient diffusivity (nutrient transport can be treated as quasi–steady), $\bm{v}$ is the velocity field (the two phases being treated as a single continuum, i.e. the model can be characterised as being one–and–a–half phase rather than fully two phase), $p$ is the pressure (i.e. Darcy’s law for a Newtonian fluid is adopted as the constitutive assumption, in keeping with many existing models) and $\mu(n;\varepsilon)$ is proportional to the tissue viscosity; $\varepsilon$ is a small parameter that appears in a number of places, associted with a key distinguished limit.

The system (2.1) is taken to hold in a finite domain $\Omega(t)$, with boundary $\Gamma(t)$ and dimensionless boundary conditions (the sign convention being such that $\kappa<0$ for a sphere)

[TABLE]

where $\varepsilon^{2}\gamma(n)$ expresses cell–cell adhesion, $\kappa$ is the mean curvature of $\Gamma(t)$, $\displaystyle\nu$ is its unit outward normal, and $q_{{{\small\nu}}}$ its normal velocity. Since we shall be concerned with what is in effect the large–time behaviour, initial conditions are not important, though

[TABLE]

would represent a plausible assumption.

In addition to the powers of $\varepsilon$ already introduced in (2.1), (2.2), the final condition that leads to the distinguished limit in question involves

[TABLE]

for constant $\mu_{1}$, associated with the physically reasonable assumption that the necrotic material is much less viscous than living tissue.

We now derive the asymptotic structure of the problem in the limit $\varepsilon\rightarrow 0$, which comprises two regions, namely a boundary–layer (thin–rim) around the tissue edge in which the living cells are concentrated, resulting from nutrient consumption therein, and a necrotic core. The former (thin–rim) is governed by a one–dimensional travelling–wave balance, described next; the latter generates, on matching to the former, the moving–boundary problem that is the subject of the rest of the paper.

The boundary–layer scalings are as follows

[TABLE]

$\nu$ denoting the outward normal distance from $\Gamma(T)$; the tangential velocity components are of $O(\varepsilon^{2})$ but are not needed in the sequel. We note that the relevant timescale is $t=O(\varepsilon^{-1})$ due to growth being confined to the thin rim (i.e. only a small minority of the volume is experiencing cell division) and dictating the velocity and pressure scalings above. Thus to leading order the travelling–wave balance

[TABLE]

holds in $\zeta<0$ with

[TABLE]

where $W_{{{\small\nu}}}=V_{{{\small\nu}}}-Q_{{{\small\nu}}}$; we note that the leading–order solution depends only on $\zeta$, so that ordinary derivative notation is indeed appropriate, this being a deduction from the boundary value problem (2.4) – (2.6) rather than an assumption; in particular $W_{\nu}$ depends only on $\zeta$ while $Q_{\nu}$ depends on $T$ and the tangential coordinates, and $V_{\nu}$, and hence $P$, depends on all the variables, this property of $W_{\nu}$ being a consequence of the Galilean invariance of the problem when $p$ is disregarded. Indeed, $P$ (and hence the mechanics) decouples, the remaining system being independent of the Darcy constitutive assumption and depending on $V_{{{\small\nu}}}$ and $Q_{{{\small\nu}}}$ only through $W_{{{\small\nu}}}$. It follows from (2.5) that (2.4) implies

[TABLE]

so that

[TABLE]

and (2.5) in consequence involves the constant

[TABLE]

which serves as a surface–tension parameter in what follows.

As $\zeta\rightarrow-\infty$ we have

[TABLE]

where the positive constants $c_{\infty}$ and $Q$ are to be determined as part of the solution to (2.4) – (2.6); $Q$ plays a crucial role in what follows, providing a surface source of material that constrasts with the volumetric sink $\lambda$. To complete the matching we must turn to the pressure $P$; the assumption (2.3) implies a non–uniformity here, but it is a consequence of $P$ decoupling that

[TABLE]

holds uniformly in the boundary layer regions (it will be clear that there is an ad–hoc flavour to this statement, but it can readily be fully justified in the sense of formal asymptotics). In consequence

[TABLE]

applies, with (2.4) implying that $W_{{{\small\nu}}}(\zeta)$ satisfies

[TABLE]

Hence as $\zeta\rightarrow-\infty$, $\varepsilon\rightarrow 0$ with $\zeta=O(1/\varepsilon)$,

[TABLE]

for positive constants

[TABLE]

We can now turn to the necrotic core region, wherein $n$ is exponentially small in $\varepsilon$, so (2.1) reduces in $\Omega(t)$ at leading order to, on setting $\bm{v}=\varepsilon\bm{V}$,

[TABLE]

with $m=1$ and

[TABLE]

The boundary conditions on (2.8) – (2.9) result from matching to (2.6) – (2.7) (i.e. by the usual matching arguments from singular–perturbation theory, (2.1) does not apply to the core). Thus (2.8) is subject on $\Gamma(t)$ to

[TABLE]

(the creation of extra volume through cell division necessitating the deformation of the surrounding tissue), while the regularising $aQ_{\nu}$ term (mathematically equivalent to kinetic undercooling, but arising for reasons distinct from those relevant in widely studied phase change problems) is associated with the mechanical interactions between the large low–viscosity core and the high–viscosity rim (i.e. the former cannot readily deform the latter despite their relative sizes, (2.3) reflecting the associated distinguished limit). We stress that the ‘kinetic–undercooling’ regulariastion arises from the thin–rim limit, whereas the ‘surface–tension’ term $\gamma_{0}\kappa$ is imposed through the original boundary condition (2.2). $b$ can be viewed as a proliferation–driven contribution to the pressure, while (2.6) – (2.9) simply give $c=c_{\infty}$; note that the first term on the right–hand side of (2.7) gives the matching condition

[TABLE]

which is automatically satisfied given (2.8) and (2.10). Finally we obtain, on translating $P$ by $b$, the moving–boundary problem

[TABLE]

where $\frac{\partial}{\partial\nu}$ denotes the outward–normal derivative. For $a=\gamma_{0}=Q=0$, this is the classical Hele–Shaw reverse (negative) squeeze–film problem (see, e.g., [22]), which is ill–posed; the $a$ term corresponds to a kinetic–undercooling regularisation and $\gamma_{0}$ to a surface–energy regularisation (cf. [23], for example, for a discussion of such regularisations in the Hele–Shaw context) - both are stabilising, but for small enough values (2.11) will be susceptible to fingering instabilities. The surface source term $Q$ is novel, making (2.11) distinct from previous Hele–Shaw formulations; since the instabilities are associated with shrinking fluid domains (due to $\lambda\mu_{1}>0$), it might be expected to play a stabilising role, though the linear–stability analysis in Section 3 quantifies the extent to which this intuition is misleading (indeed, the growth associated with $Q$ has the potential to stabilise by filling inward–facing fingers but also allows more room for the destabilising mechanism to manifest itself). Setting $P=u$, $a=\alpha$, $Q_{{{\small\nu}}}=V$, $\beta=\gamma_{0}/a$ and scaling such that $\mu_{1}=\lambda=1$ yields (1.1) – (1.3), with $t$ reinstated as the time variable, in place of (and equivalent to) $T$ above.

We regard the rigorous justification of the above ‘thin–rim’ limiting process for (1.1) – (1.3), from (2.1) as a worthwhile open problem. The limit has features in common with both sharp–interface and thin–film limits, without being exactly of either type (for example an interface is present both before and after the limit is taken, but the conditions that apply thereon are modified by the limit process).

3 Model analysis

3.1 Radially symmetric solutions and their stability

3.1.1 Circular solutions

We limit ourselves here, for brevity, to the two–dimensional case, though the three–dimensional one is of course also tractable. Taking the cylindrically symmetric solution to be $u=u_{c}(r,t)$ in $\Omega$, with $r=s_{c}(t)$ corresponding to $\Gamma$, (1.1) – (1.3), reduces to

[TABLE]

so that

[TABLE]

for initial data $s_{c}(0)=s_{0}$.

3.1.2 Linear stability

Setting

[TABLE]

where $r=s(\theta,t)$ parametrises $\Gamma$, and retaining only those terms that are linear in $\delta$ implies

[TABLE]

[TABLE]

on projecting $\Gamma$ onto $r=s_{c}(t)$ in the usual way. Fourier decomposing in the form

[TABLE]

for non–negative integer $n$ then implies that

[TABLE]

so that

[TABLE]

Hence

[TABLE]

the former corresponding to perturbing $s_{0}$ in (3.2), and the latter arising because the $n=1$ modes simply correspond to small shifts in the coordinate origin. $\dot{B}_{n}$ is negative for $n\geq 2$, and the circular solution is hence linearly stable, if

[TABLE]

i.e. the steady state $s_{c}=2Q$ is linearly stable if

[TABLE]

In view of (3.3), it is instructive to determine the faster growing mode, i.e. value of $n$ that maximises

[TABLE]

whenever the latter is positive, and this may change as $s_{c}$ evolves; for $s_{0}$ sufficiently small (and less than $2Q$), all modes will initially decay, but when (3.5) is violated the $n=2$ mode will grow for sufficiently large $t$, and other modes may subsequently become unstable. For a wide range of parameters the $n=2$ mode will be the most unstable one for all $t$, but an instructive regime in which this is not the case arises for $\beta$ small. Thus for $\beta=0$ all modes with $n\geq 2$ are unstable, with the limit $n\rightarrow\infty$ giving the fastest growth rate, namely

[TABLE]

(this expression clarifies the role of $\alpha$ as a regularising parameter). For $\beta\ll 1$ and $n\gg 1$, (3.6) is given asymptotically by

[TABLE]

for $n=O(\beta^{-\frac{1}{3}})$, so the fastest growing mode at a given $t$ has

[TABLE]

Because this quantity changes with $t$, the dominant wavelength observed at large time will presumably be dependent on how the relative magnitudes of the modes in the initial data feed through into the ultimate nonlinear evolution, and it is important to note that the leading term in (3.7) is independent of $n$. While this instability mechanism is closely related to the well known Saffman–Taylor instability in Hele–Shaw and Darcy flows, its relevance to tissue growth applications is less widely recognised, a key point being that (in contrast to the Hele–Shaw case) instabilities can arise for growing domains (due to the surface source $Q$) as well as shrinking ones; see [18, 20] for a detailed classification of instability mechanisms for such tissue–growth moving–boundary problems.

3.2 Other properties

We briefly note some other properties here (additionally the appendix records a thin–film limit). The result

[TABLE]

which readily follows on applying the divergence theorem, expresses the overall mass balance arising from the surface source and the volumetric sink. Thus tissue can enhance its growth by maximising its surface area for a given volume, making spheres the most sub–optimal shapes in this sense, and providing an interpretation of interface fingering (this does not characterise a fingering mechanism but does clarify why such instabilities may be favoured).

Secondly, steady state solutions satisfy

[TABLE]

[TABLE]

and this free–boundary problem is itself of interest. The circular case is addressed in the linear–stability section, Section 3.1.2, but there are presumably also non–trivial (but probably unstable) steady configurations. Applying the linear stability analysis may, by identifying bifurcation points, be instructive in classifying these solutions and in facilitating weakly nonlinear analyses but, given that they seem likely to be unstable, we shall not pursue their analysis here. Nevertheless, we highlight their potential interest from the free–boundary problem perspective.

Finally, we record the interfacial dynamics limit that arises for small domains (though of course nevertheless much larger than the rim thickness) with $\alpha\gg 1$, whereby the leading–order problem involves the evolution of the free boundary only. Setting $\alpha=1/\varepsilon$, $\beta=\varepsilon\gamma$, and rescaling $\bm{x}\rightarrow\varepsilon\bm{x}$, $V\rightarrow\varepsilon V$, $Q\rightarrow\varepsilon Q$ to give

[TABLE]

[TABLE]

so that

[TABLE]

and the divergence theorem (or, equivalently, the solvability condition at $O(\varepsilon^{2})$) gives (as above)

[TABLE]

implying the novel non–local interfacial dynamics law of forced mean–curvature–flow type

[TABLE]

that involves the total mean curvature; this becomes more explicit in two dimensions, whereby

[TABLE]

We note that the exact result (3.8) is readily recovered from the limit case (3.9).

4 A diffuse-interface model

In this section we derive a diffuse–interface formulation of (1.1) – (1.3). In the diffuse–interface paradigm we approximate the interface $\Gamma(t)$ by an interfacial region $\Gamma_{\varepsilon}(t)$ with width $C\varepsilon\ll 1$. To this end a phase–field parameter $\varphi$ is introduced such that, within the thin interfacial region, $\varphi$ varies smoothly from $\varphi=-1$ to $\varphi=1$, while outside it $|\varphi|=1$, thus we can define $\Gamma_{\varepsilon}$ by

[TABLE]

We set $\varphi=1$ in the region inside $\Gamma_{\varepsilon}$ and $\varphi=-1$ exterior to $\Gamma_{\varepsilon}$. In addition we let $D$ denote a bounded domain in $\mathbb{R}^{d}$, with a polygonal boundary, such that $D$ is chosen large enough to contain $\Gamma_{\varepsilon}(t)$ for all $t$. We denote the outward unit normal to $\partial D$ by \tilde{\mbox{\boldmath\displaystyle\nu\unboldmath}}.

To derive a diffuse–interface approximation of (1.1), (1.2) we first introduce the variational form:

[TABLE]

We now follow the techniques described in [1] to obtain a diffuse–interface approximation of (4.1). Setting $\tilde{u}\in D$ to be a diffuse–interface approximation to $u\in\Omega$ we have

[TABLE]

where

[TABLE]

Using (4.2) we arrive at the following diffuse–interface approximation of (4.1):

[TABLE]

Here $\tilde{\kappa}$ approximates an extension of $\kappa$ to $\Gamma_{\varepsilon}$ such that at $t=0$, $\tilde{\kappa}$ is obtained by constantly extending $\kappa$ in the normal direction via the signed distance function to $\Gamma(0)$, while for $t>0$ we use the Allen–Cahn approximation of forced–mean–curvature flow, presented below in (4.5), to set

[TABLE]

For $\beta=0$, the strong convergence of the solution $\tilde{u}$ of (4.3) to the solution $u$ of (4.1), as the approximation parameter $\varepsilon$ tends to zero, is shown in [1], while in [4] a convergence result of the form

[TABLE]

for a constant $C>0$ independent of $\varepsilon$, and some $p>0$ is shown for a slightly different diffuse interface approximation in which $\delta(\varphi)=|\nabla\varphi|$.

We now present a diffuse–interface approximation or, as it is more commonly known, a phase–field approximation, of (1.3). Since (1.3) describes forced–mean–curvature flow, we approximate it using an Allen–Cahn equation with forcing,

[TABLE]

with boundary data

[TABLE]

Here $\tilde{u}$ is the solution to (4.3), $W(\varphi)$ is the double–obstacle potential, see [3],

[TABLE]

where $I_{[-1,1]}(u)$ is the indicator function

[TABLE]

and the factor of $\pi/4$ on the right hand size of (4.5) is chosen such that as $\varepsilon\rightarrow 0$, we recover the velocity law (1.3).

We define the initial data for the phase–field parameter $\varphi$ by

[TABLE]

where $d(x)$ is the signed distance function to $\Gamma(0)$. This initial condition corresponds to a diffuse–interface of width $\varepsilon\pi$, with the level set $\varphi=0$ approximating $\Gamma(0)$.

Combining (4.3) with the weak formulation of (4.5), (4.6) yields the following diffuse–interface approximation of (1.1) – (1.3):

Given $\Gamma(0)\in\mathbb{R}^{d}$, find $(\tilde{u},\varphi)$ such that $\varphi$ satisfies the initial data (4.7) and for all $t\in(0,T)$, $\varphi(t)\in K$,

[TABLE]

and

[TABLE]

where

[TABLE]

The well-posedness of the variational problem (4.8), (4.9) is an open problem. Setting $\beta=0$ in (4.8) and adding a time regularisation term, $\varepsilon\tilde{u}_{t}$, to the left hand side of the resulting equation would result in a diffuse–interface system similar to the one studied in [8] in which a parabolic variational inequality is coupled to a degenerate diffusion equation. In [8] the authors use a regularisation technique to prove an existence and uniqueness result for the coupled system.

5 Finite–element approximations

In this section we derive finite–element approximations to (1.1) – (1.3), and to the diffuse–interface model (4.8), (4.9).

5.1 Parametric approach

We partition the time interval $[0,T]$ into $N$ equidistant steps, $0=t_{0}<t_{1}<\ldots<t_{N-1}<t_{N}=T$, so that for $i=1,\ldots,N$, $\Delta t:=t_{n}-t_{n-1}$. We let $\Gamma_{h}^{n}$ be a polyhedral surface approximating $\Gamma(t_{n})$, with ${\Gamma_{h}^{n}=\bigcup_{i=1}^{I}\overline{\sigma}_{i}^{n}}$, where $\{\sigma_{i}^{n}\}^{I}_{i=1}$ is a family of disjoint open simplices (straight line segments in $\mathbb{R}^{2}$ and triangles in $\mathbb{R}^{3}$). We let $\mathcal{J}$ be the set of nodes of $\{\sigma_{i}^{n}\}^{I}_{i=1}$ and let $\{\bm{q}_{j}^{n}\}_{j\in\mathcal{J}}$ be the coordinates of these nodes.

We denote the interior of $\Gamma_{h}^{n}$ by $\Omega_{h}^{n}$ and we partition $\Omega_{h}^{n}$ into a family of disjoint open simplices $\{\mu_{k}^{n}\}^{K}_{k=1}$ (triangles in $\mathbb{R}^{2}$ and tetrahedra in $\mathbb{R}^{3}$) such that $\Omega_{h}^{n}=\bigcup_{k=1}^{K}\overline{\mu}_{k}^{n}$. Here $\{\mu_{k}^{n}\}^{K}_{k=1}$ is chosen such that the nodes and the faces (straight line segments in $\mathbb{R}^{2}$ and triangles in $\mathbb{R}^{3}$) of the elements of $\{\mu_{k}^{n}\}^{K}_{k=1}$ that make up the boundary of $\Omega_{h}^{n}$ coincide with nodes and simplices of $\{\sigma_{i}^{n}\}^{I}_{i=1}$. We let $\mathcal{L}$ be the set of nodes of $\{\mu_{k}^{n}\}^{K}_{k=1}$ and let $\{\bm{p}_{l}^{n}\}_{l\in\mathcal{L}}$ be the coordinates of these nodes.

We define {\mbox{\boldmath\displaystyle\nu\unboldmath}}^{n} to be the outward unit normal to $\Gamma_{h}^{n}$, such that {\mbox{\boldmath\displaystyle\nu\unboldmath}}_{i}^{n}:={\mbox{\boldmath\displaystyle\nu\unboldmath}}^{n}|_{\sigma_{i}}, for $i\in 1,\ldots,I$ and for $n\in 0,\ldots,N-1$ we introduce the finite–element spaces

[TABLE]

and

[TABLE]

with $\displaystyle{\{\chi_{j}^{n}\}_{j\in\mathcal{J}}}$ denoting the standard basis for $V_{h}(\Gamma_{h}^{n})$.

We describe $\Gamma_{h}^{n}$ by the piecewise linear vector function $\bm{X}^{n}_{h}$ such that $\bm{X}_{h}^{n}\in[V_{h}(\Gamma_{h}^{n-1})]^{d}$.

We denote the $L^{2}$–inner product over ${\Omega_{h}^{n}}$ by ${(f,g)_{\Omega_{h}^{n}}:=\int_{{\Omega_{h}^{n}}}fgdx}$ and the discrete $L^{2}$–inner product by

[TABLE]

with $\Pi^{h}:C(\overline{{\Omega_{h}^{n}}})\to S_{h}$ denoting the interpolation operator, such that $(\pi^{h}\eta)(\bm{p}_{l})=\eta(\bm{p}_{l})$ for all $l\in\mathcal{L}$. We adopt similar notation for $(u_{h},v_{h})_{\Gamma_{h}^{n}}^{h}$.

We now present a finite element approximation of (1.1) – (1.3) in which we used (4.1) in order to approximate (1.1) – (1.2) and we followed the authors in [9] to approximate (1.3).

Given $\Gamma^{0}$ and the identity function $\bm{X}^{0}_{h}\in[V_{h}(\Gamma_{h}^{0})]^{d}$ on $\Gamma^{0}$, for $n=0\to N-1$ find $\{u_{h}^{n},\bm{X}_{h}^{n+1}\}\in V_{h}(\Omega_{h}^{n})\times[V_{h}(\Gamma_{h}^{n})]^{d}$ such that

[TABLE]

and

[TABLE]

where \displaystyle{{\mbox{\boldmath\displaystyle\omega\unboldmath}}^{n}_{h}:=\sum_{j\in\mathcal{J}}{\mbox{\boldmath\displaystyle\omega\unboldmath}}_{j}^{n}\chi_{j}^{n}(x)\in[V_{h}(\Gamma_{h}^{n})]^{d}} with

[TABLE]

with $|\sigma_{i}|$ denoting the measure of $\sigma_{i}$, $\mathcal{T}_{j}^{n}:=\{\sigma_{i}^{n}:\bm{q}_{j}^{n}\in\overline{\sigma_{i}^{n}}\}$ and $\Lambda_{j}^{n}:=\cup_{\sigma_{i}\in\mathcal{T}_{j}^{n}}\overline{\sigma_{i}^{n}}$, such that {\mbox{\boldmath\displaystyle\omega\unboldmath}}_{j}^{n} can interpreted as a weighted normal defined at the node $\bm{X}_{h}^{n}(\bm{q}_{j}^{n})=\bm{q}_{j}^{n}$ of the surface $\Gamma_{h}^{n}$.

Noting from (1.3) that $\beta\kappa=V-u/\alpha$, the approximation, $\displaystyle{\beta\kappa_{h}^{n}=\beta\sum_{j\in\mathcal{J}}\kappa_{j}^{n}\chi_{j}^{n}(x)\in V_{h}(\Gamma_{h}^{n})}$, of the curvature is given by

[TABLE]

and we set $\kappa_{h}^{0}\in V_{h}(\Gamma_{h}^{0}):=\Pi^{h}(\kappa_{0})$, where $\kappa_{0}$ is the curvature of $\Gamma(0)$.

5.2 Diffuse–interface approach

We partition the time interval $[0,T]$ into $N$ equidistant steps, $0=t_{0}<t_{1}<\ldots<t_{N-1}<t_{N}=T$, so that for $i=1,\ldots,N$, $\Delta t:=t_{n}-t_{n-1}$. Let $\{{\cal T}^{h}\}_{h>0}$ be a family of partitionings of $D$ into disjoint open simplices $\sigma$ with $h_{\sigma}:={\rm diam}(\sigma)$ and $h:=\max_{\sigma\in{\cal T}^{h}}h_{\sigma}$.

Associated with ${\cal T}^{h}$ is the finite element space

[TABLE]

We denote the standard basis for $S_{h}(D)$ by $\displaystyle{\{\psi_{k}^{n}\}_{k=1}^{M}}$ and we introduce

[TABLE]

Using standard finite element approximations of (4.8) and (4.9) we obtain the following.

Given an approximation $\varphi_{h}^{0}\in K_{h}$ of $\varphi_{0}\in K$, for $n=0\to N-1$ find $\{\tilde{u}_{h}^{n},\varphi_{h}^{n+1}\}\in S_{h}\times K_{h}$ such that for all $\eta_{h}\in S_{h}$,

[TABLE]

and

[TABLE]

where $\tilde{u}_{h}^{n}(x)$ is set to zero if $\varphi_{h}^{n}(x)=-1$, $(\cdot,\cdot)_{D}^{h}$ is discrete inner product over $D$ and noting (4.4) the approximation $\displaystyle{\beta\tilde{\kappa}_{h}^{n}=\beta\sum_{j=1}^{M}\tilde{\kappa}_{j}^{n}\psi_{j}\in S_{h}}$, of the curvature is given by

[TABLE]

with $\tilde{\kappa}_{h}^{0}\in S_{h}=\Pi^{h}(\tilde{\kappa}_{0})$, where $\kappa_{0}$ is the curvature of $\Gamma(0)$ and $\tilde{\kappa}_{0}$ is such that $\tilde{\kappa}_{0}(x):=\kappa_{0}(\hat{p}(x))$, where $\hat{p}(x)$ denotes the closest point projection of a point $x\in D$ onto $\Gamma(0)$.

6 Numerical Results

In this section we first briefly describe the implementation of the finite element schemes presented in Section 5, and then we present a number of numerical simulations. The results displayed are visualised in Paraview [16].

6.1 Implementation of the finite element schemes

6.1.1 Implementation of the parametric scheme in $\mathbb{R}^{2}$

For $\Omega_{h}^{n}\in\mathbb{R}^{2}$ the parametric finite–element scheme (5.1), (5.2) was implemented with the finite–element toolbox ALBERTA 2.0, see [21]. As we are solving on an evolving domain, care has to be taken to ensure good mesh properties for the bulk mesh $\Omega_{h}^{n}$. To achieve this we implemented Algorithm 1, as introduced in [10] in which the authors present a variant of the so-called the DeTurck trick using the the harmonic map heat flow that leads to computational meshes of high quality for both for the bulk mesh $\Omega_{h}^{n}$ and its boundary $\Gamma_{h}^{n}$.

A simple mesh refinement algorithm is also presented in [10], in which we seek to refine element $\mu^{n}\in\Omega_{h}^{n}$ if $|\mu^{n}|\geq 1.5\overline{|\mu^{0}|}$ and coarsen element $\mu^{n}\in\Omega_{h}^{n}$ if $|\mu^{n}|\leq 0.5\overline{|\mu^{0}|}$, where $\mu^{0}\in\Omega_{h}^{0}$, $|\mu^{n}|$ is the area of element $\mu^{n}$, and $\overline{|\mu^{0}|}$ is the mean area of all elements at time $t=0$.

When the mesh reached a point at which the algorithm presented in [10] was unable to maintain the required mesh quality, we used the re–meshing software GMSH, see [13], to re–mesh $\Omega_{h}^{n}$, whilst keeping the nodes on $\Gamma_{h}^{n}$ fixed. We performed this re–meshing procedure when

[TABLE]

where $H(\mu^{n})$ is the length of the longest edge of $\mu^{n}$ and $h(\mu^{n})$ is the length of the shortest edge.

6.1.2 Implementation of the parametric scheme in $\mathbb{R}^{3}$

For $\Omega_{h}^{n}\in\mathbb{R}^{3}$ the parametric finite–element scheme (5.1), (5.2) was implemented with the finite–element toolbox ALBERTA 3.1; the quality of the mesh was maintained by moving the nodes with a velocity $\bm{u}$, using the harmonic extension method. We compute $\bm{u}$ by solving a standard finite element approximation to

[TABLE]

where $\bm{v}$ is the velocity of $\Gamma$. This extends the velocity on $\Gamma$ to $\Omega$ in a ‘smooth’ manner.

6.1.3 Implementation of the diffuse–interface scheme

The diffuse–interface finite–element scheme (5.3), (5.4) was implemented with the finite–element toolbox ALBERTA 2.0. Since the interfacial thickness is proportional to $\varepsilon$, in order to resolve the interfacial layer, we need to choose the mesh size $h$ such that $h\ll\varepsilon$, see [7] for details. Away from the interface, $h$ can be chosen larger and hence adaptivity in space can be used to speed up the computations. We use the built in algorithms from ALBERTA 2.0 for adaptivity to implement a mesh refinement strategy in which the following meshes are constructed: a fine mesh in the interfacial region where $|\varphi_{h}^{n}|<1$, a coarse mesh exterior to the tumour where $\varphi_{h}^{n}=-1$ and a standard sized mesh in the interior of the tumour where $\varphi_{h}^{n}=1$, see the left hand plot in Figure 2 in which an enlarged section of the diffuse–interface mesh associated with the results in Figure 10 at $t=45$ is displayed, in addition an enlarged section of the corresponding parametric mesh is displayed in the right hand plot. We denote the maximum diameter of the triangles in the three meshes by $h_{max,f}=\max_{\sigma\in\mathcal{T}_{f}^{n}}h_{\sigma},h_{max,m}=\max_{\sigma\in\mathcal{T}_{m}^{n}}h_{\sigma}$ and $h_{max,c}=\max_{\sigma\in\mathcal{T}_{c}^{n}}h_{\sigma}$, where

[TABLE]

and

[TABLE]

For more details on a similar mesh refinement strategy see [2].

6.2 Generalisation of the model

In the numerical results presented below we consider the following generalisation of (1.1) – (1.3)

[TABLE]

Here $\gamma\in\mathbb{R}$ is a nonnegative regularisation parameter, such that setting $\gamma=0$ recovers the original model (1.1) – (1.3). We introduce this generalised model for numerical purposes as it enables us to run computations with $\beta=0$, thus simplifying the numerical schemes by removing the curvature term from the Robin boundary conditions for $u$, whilst maintaing a velocity law of forced mean curvature flow.

6.3 Parameter values

Unless otherwise specified, in all simulations in $\mathbb{R}^{2}$ the following hold: we set $Q=1.0$ and the initial geometry $\Gamma(0)$ was set to be an ellipse with length $0.5$ and height $1.0$, in the parametric examples we set $\Delta t=10^{-3}$ and the mesh size at $t=0$ was taken to be $h\approx 0.022$. In the diffuse–interface simulations we set $h_{max,c}\approx 1.77$, $h_{max,m}\approx 0.02$, and used the values of $\varepsilon,~{}\Delta t$ and $h_{max,f}$ given in Tables 1 and 2.

6.4 Simulations for $\gamma=0$

In the first set of simulations, Figures 3 – 8, we take $\gamma=0$ in the generalised model to obtain results for the original model (1.1) – (1.3).

6.4.1 Radially symmetric solutions

To test the accuracy of the numerical schemes, we set $\Gamma(t)\subset\mathbb{R}^{2}$ to be a circle with radius $R(t)$ and we express $u$ in polar coordinates such that $u(r,\theta)=u(r)$. In this setting (1.1) – (1.3) reduces to

[TABLE]

[TABLE]

so that

[TABLE]

For some choices of parameters, typically $\alpha,\,\beta\ll 1$ and $Q=O(1)$, instabilities occurred in the parametric scheme that gave rise to non-radially symmetric solutions, these instabilities seem likely to have stemmed from the initial triangulation of $\Omega(0)$. In the following results we avoided such parameter combinations.

In Figure 3 we compare the radii obtained from the finite element approximations of the parametric and the diffuse–interface formulations for $\gamma=0$ to the corresponding analytical solutions. Similarly we compare the pressure solutions obtained from the parametric scheme and the diffuse–interface scheme to the analytical solution.

We consider four combinations of $\alpha$ and $\beta$ and set $Q=R_{0}=1.5$. In the upper plots we display the radius of the circle for the parametric model (red line), the radius of the circle for the diffuse–interface model, obtained from the zero–level line of $\varphi$, (black line), and the radius of the analytical solution of (6.6) (blue dots). In the lower plots we display $u_{h}$, $\tilde{u}_{h}$, and $u$ at $t=0.5$, where $u$ is the solution of (6.4), again, we use a red line for the parametric scheme, a black line for the diffuse–interface scheme, and blue dots for the analytical solution. We display the results for $\gamma=0$, $\alpha=1$, $\beta=0.1$ in the first column, $\gamma=0$, $\alpha=\beta=0.1$ in the second column, $\gamma=0$, $\alpha=\beta=1$ in the third column and $\gamma=0$, $\alpha=0.1$, $\beta=1$ in the fourth column. In the parametric simulations the mesh size at $t=0$ was taken to be $h\approx 0.1$.

From Figure 3 we see that for all four combinations of $\alpha$ and $\beta$ there is very good agreement between the solutions of the parametric scheme and the analytical solution, while the solutions of the diffuse–interface scheme compare reasonably well to the analytical solution for $\gamma=0$, $\alpha=1$, $\beta=0.1$ and $\gamma=0$, $\alpha=\beta=0.1$ but not well for $\gamma=0$, $\alpha=\beta=1$ and $\gamma=0$, $\alpha=0.1$, $\beta=1$. We believe the poor approximation of the diffuse–interface scheme for $\beta=O(1)$ is due to the fact that the scheme is only an approximation to forced mean curvature flow and hence larger errors occur if the curvature term in the evolution becomes more dominant. Based on these results in the subsequent simulations we solved the diffuse–interface model only with $\beta=0.1$.

6.4.2 Linear stability

In Figure 4 we display results from the parametric scheme (5.1), (5.2) relating to the linear stability result, (3.5). We set $Q=1$ and $\beta=0.1$ and set the initial geometry to be a perturbation of a circle with radius $R(0)=2Q=2$ such that $R(\theta,0)=2Q+0.25\sin(\theta)\cos(9\theta)$, (left plot) and we set $h\approx 0.1$. We consider two values of $\alpha$; $\alpha=0.68$ (centre plot) such that $3\alpha\beta>2$ and $\alpha=0.65$ (right plot) such that $3\alpha\beta<2$, in each case we display the solution $u_{h}$ at $t=700$. From the figure we see that when $3\alpha\beta>2$ the stable state in which $\Gamma(t)$ is a circle of radius $R=2Q=2$ is obtained, whereas for $3\alpha\beta<2$ the radial symmetry of the problem is lost. We note that for this parameter set with $\alpha=0.65$, $n=2$, recall (3.3), is the only unstable mode, the simulations leading to the same (reflectional) symmetries as those of an ellipse.

6.4.3 Dependence on $\alpha$ and $\beta$

We now investigate the morphologies that arise using differing combinations of $\alpha$ and $\beta$.

We first note that, for the initial data chosen, an ellipse with length $0.5$ and height $1.0$, our simulations show that for $\gamma=0$, $\alpha=\beta=1.0$ the geometry tends to a radially symmetric steady state, with the radius of the resulting circle decreasing as $\beta$ increases; this radius can be calculated by setting $R^{\prime}(t)=0$ in (6.6).

In Figure 5 we display the solution $u_{h}$ at $t=0,7,10,12.6$, computed using the parametric scheme (5.1), (5.2) with $\gamma=0$, $\alpha=0.1$ and $\beta=1$

In Figure 6 we compare the parametric and diffuse–interface numerical solutions obtained by setting $\gamma=0$, $\alpha=1$ and $\beta=0.1$. The solution $u_{h}$ given by the parametric scheme (5.1), (5.2) is displayed in the top row, while the solution $\tilde{u}_{h}$ given by the diffuse–interface scheme (5.3), (5.4) is displayed in the middle row. The bottom row displays the order parameter $\varphi_{h}$ (in red and blue) given by the diffuse–interface scheme and $\mathbf{X}_{h}$ (in white) given by the parametric scheme. The solutions are displayed at times $t=0$ (left) $t=11$ (centre) and $t=23$ (right). At $t=11$ there is good agreement between the parametric and diffuse–interface schemes, though for $t=23$ the difference becomes more pronounced.

The results obtained by setting $\gamma=0$, $\alpha=0.1$ and $\beta=0.1$ are displayed in Figure 7. For this choice of parameters the agreement between the parametric and diffuse–interface schemes is quite pronounced; in particular, the evolution of the diffuse–interface solution is slower than that of the parametric solution, however, the geometries of the two solutions compare well if the diffuse–interface solution is displayed at a different time to the parametric solution. To this end in Figure 7, the diffuse–interface results are displayed at $t=0,4,5$ while the parametric results are displayed at $t=0,3.5,4.4$. We believe the disparity between the two sets of solutions is predominantly caused by the diffuse–interface approximation of the pressure that is used in the Allen-Cahn approximation of the velocity law, since for $\alpha\ll 1$, a small error in the approximation of $u$ on $\Gamma(t)$ can result in a significant error in the velocity of $\Gamma(t)$. When we solved the model with $\beta=0$ with the same set of parameters we saw good agreement between the parametric and diffuse–interface schemes, see Figure 12 below, and thus we conclude that the diffuse interface approximation of the curvature term in the Robin boundary conditions is the main reason for the differences.

6.4.4 Thin–film limit (see the appendix)

In Figure 8 we display results from the parametric scheme (5.1), (5.2), relating to the thin–film limit from the appendix. Here the initial geometry $\Gamma(0)$ is given by an ellipse of length $1.0$ and height $0.1$. We set $Q=0.1$, $\gamma=0$ and $\alpha=\beta=0.1$. We chose $\Delta t=10^{-4}$ and $h\approx 0.005$ in order to maintain a stable evolution.We present the solutions at $t=0,0.2,0.3$, demonstrating how the large–aspect–ratio of the domain is lost, consistent with the analysis of the appendix.

6.5 Simulations for $\beta=0$ and $\gamma\neq 0$

We conclude our numerical results with some simulations with $\beta=0$ and $\gamma\neq 0$, such that the curvature term from the Robin boundary conditions for $u$ is removed, whilst a velocity law of forced mean curvature flow is maintained.

6.5.1 Radially symmetric solutions in $\mathbb{R}^{2}$

Figure 9 presents radially symmetric results with $\beta=0$ and $Q=R_{0}=1.5$ displayed in the same format as the results in Figure 3. In the parametric simulations the mesh size at $t=0$ was taken to be $h\approx 0.1$. As in Figure 3 we see that for all four combinations of $\alpha$ and $\gamma$ there is very good agreement between the solutions of the parametric scheme and its analytical solution. In addition we see that for $\beta=0$ , $\alpha=1$, $\beta=0.1$ and $\alpha=\gamma=0.1$ the diffuse–interface solution is much closer to the analytical solution than it was for the original model, we also note that for $\beta=0$ , $\alpha=1$, $\gamma=0.1$ we set $\varepsilon=0.04$ rather than the larger value $\varepsilon=0.075$ that was used for the original model. Similar to the original model, the results of the diffuse–interface model with $\beta=0$, $\alpha=\gamma=1$ and $\alpha=0.1$, $\gamma=1$ are not good and so in the subsequent simulations we will not include any diffuse–interface results using these parameter combinations.

6.5.2 $\mathbb{R}^{2}:$ $\beta=0,$ $\alpha=1.0,~{}\gamma=1.0$

As in the case for $\gamma=0$, with the initial data chosen, an ellipse with length $0.5$ and height $1.0$, our simulations show that for $\beta=0$, $\alpha=\gamma=1.0$ the geometry tends to a radially symmetric steady state, with the radius of the resulting circle decreasing as $\beta$ increases.

6.5.3 $\mathbb{R}^{2}:$ $\beta=0,$ $\alpha=0.1,~{}\gamma=1.0$

The parameters $\gamma=0$, $\alpha=0.1$ and $\beta=1.0$ gave rise to the geometries displayed in Figure 5, however for the parameters $\beta=0$, $\alpha=0.1$ and $\gamma=1.0$, such that the curvature term in the Robin boundary conditions for $u$ is removed, the initial ellipse with length $0.5$ and height $1.0$ shrinks to a point.

6.5.4 $\mathbb{R}^{2}:$ $\beta=0,$ $\alpha=1.0,~{}\gamma=0.1$

Figure 10 presents results with $\beta=0$, $\alpha=1$ and $\gamma=0.1$, displayed in the same format as the results in Figure 6. The solutions are displayed at times $t=0,15,30,45$ (columns one to four respectively). As in Figure 6, we see good agreement between the parametric and diffuse–interface schemes. When comparing Figure 10 with Figure 6 we see that the additional curvature term in the Robin boundary conditions for $u$ gives rise to more a rounded structure in the geometry.

Figure 11 shows the merging of two circular tumours. The results are obtained from the diffuse–interface scheme with $\beta=0$, $\alpha=1.0$, $\gamma=0.1$ and $\varepsilon=0.04$ and the solutions are displayed at $t=0,1,20,30$. The initial geometry is given by two circles of radius $1.0$, with centres at $(1.3,0)$ and $(-1.3,0)$. This example highlights the power of the diffuse–interface model as the parametric model would not be able to handle directly the topology change that takes place when two tumours merge.

6.5.5 $\mathbb{R}^{2}:$ $\beta=0,$ $\alpha=0.1,~{}\gamma=0.1$

Figure 12 presents results with $\beta=0$, $\alpha=\gamma=0.1$ displayed at $t=0,3,7$. We see good agreement between the parametric and diffuse–interface schemes. This is in contrast to the corresponding simulations with $\gamma=0$, $\alpha=\beta=0.1$, displayed in Figure 6, where the results from the parametric and diffuse–interface schemes are displayed at different times.

In Figure 13 we display results using the parametric scheme with $\beta=0$, $\alpha=\gamma=0.1$. Here the initial geometry $\Gamma(0)$ is purposefully chosen to reduce symmetry. To create the initial geometry we perturb an ellipse of length $0.5$ and height $1.0$ by a distance of $0.2\sin(6\theta)$ in the normal direction, where $\theta$ is the polar angle. We present the solutions at $t=0,2.3,4.6,6.9$.

6.5.6 $\mathbb{R}^{2}:$ $\beta=0,$ $\alpha=1.0,~{}\gamma=0.05,~{}0.025,~{}0.01$

In Figure 14 we show the effect that reducing $\gamma$ has on the evolution of the tumour. We set $\beta=0$, $\alpha=1.0$ and consider three values of $\gamma$: $\gamma=0.05$ at $t=20$ (left plot), $\gamma=0.025$ at $t=14$ (centre plot) and $\gamma=0.01$ at $t=11.6$ (right plot). From this figure we see that as $\gamma$ is reduced the extent to which the tumour extends in the horizontal direction is also reduced.

6.5.7 $\mathbb{R}^{3}:$ $\beta=0,$ $\alpha=1.0,~{}\gamma=0.1$

We set $\beta=0$, $\alpha=1.0$, $\gamma=0.1$ and $Q=1.25$ and the initial geometry $\Gamma(0)$ is given by the oblate spheroid with equation $\frac{x^{2}}{1.0^{2}}+\frac{y^{2}}{0.5^{2}}+\frac{z^{2}}{1.0^{2}}=1$ .

In the parametric examples the mesh size was taken to be $h\approx 0.15$ at $t=0$, and we set $\Delta t=5\times 10^{-3}$. In the diffuse–interface examples we set $\varepsilon=0.1$, $h_{max,f}\approx 0.011$, $h_{max,m}\approx 0.045$, $h_{max,c}\approx 8.485$ and we set $\Delta t=10^{-3}$.

Figure 15 displays results from the parametric scheme at $t=0,7,13$. In this figure we display the surface $\Gamma_{h}$ in three different orientations: looking down the $x$ axis (first row), down the $y$ axis (second row), and a cross section in the plane $z=0$ (third row). We also display the solution $u_{h}$ on a cross section in the plane $z=0$ (fourth row).

In Figure 16 we display results from the diffuse–interface scheme. We display the zero level surface of $\varphi_{h}$ in three orientations: looking down the $x$ axis (first row), down the $y$ axis (second row), and a cross section in the plane $z=0$ (third row). We also display plots of the solution $\tilde{u}_{h}$ on the zero level surface of $\varphi_{h}$ on a cross section in the plane $z=0$ (fourth row).

6.6 Summary of the numerical simulations

From the above results we see that the model gives rise to complicated structures developing from simple elliptical geometries. Comparing Figure 12 with Figure 13 we see the complexity of the structure increases with the complexity of the initial geometry. The additional curvature term in the Robin boundary condition for the pressure that is present when $\beta\neq 0$, results in a more compact rounded geometry, compare Figures 6 and 7 with the corresponding Figures 10 and 12. From Figures 6, 7, 10 and 12 we see that significantly larger values of the pressure, $u$, on the boundary are achieved with $\alpha=1$ compared to $\alpha=0.1$, whereas the difference between the minimum values in the interior is much smaller.

In Figure 11 simulations from the diffuse–interface model show two tumours merging to form a single tumour. The change in topology that occurs in this example is naturally handled by the diffuse–interface model and demonstrates the power of this model, in addition it illustrates the complicated pathologies that can result even for the current simple tissue–growth model

7 Discussion

We conclude with some observations, noting possible extensions to the work and highlighting open problems.

The formulation (1.1) – (1.3) derived and analysed above represents arguably the simplest macroscale spatially structural model for tissue growth (another of comparable complexity replaces Darcy flow with Stokes flow) and hence bucks the current trend in which, for good reason, numerous additional mechanisms are incorporated with the goal of achieving better biological realism. It nevertheless exhibits growth morphologies sharing features of those much more complicated models, allowing the associated mechanisms much more readily to be identified and mitigating against possible pitfalls of assuming responsibility for behaviour that is in fact rather generic to specific biological processes.

The biological aspects implicit in (2.1) pertain primarily to the constitutive assumptions made on $k_{b}(c)$, $k_{d}(c)$, $D(n)$ and $K(c)$, on the interpretation of $\lambda$ and on the choice of Darcy flow. Obvious generalisations would include the inclusion of active cell–cell interactions (necessitating the adoption of genuinely multiphase models), the adaption of more complex and realistic constitutive laws for the mechanics and the proper treatment of the exterior domain - notably, replacing the Dirichlet condition on $c$ in (2.2) by describing exterior nutrient transport would allow the inclusion of Mullins–Sekerka instabilities as well as those of Saffron–Taylor type captured above (whereby the exterior is in effect inviscid and fingers are prone to develop due to a self-reinforcing mechanism whereby for the sink it is easier to pull in tissue in directions for which there is less to pull). We believe the current framework provides a firm and well–characterised basis on which to incorporate such additional effects.

From a mathematical perspective, open problems for rigorous analysis abound in this context, including the derivation of the thin–rim limit, the analysis of the travelling–wave boundary value problem and the derivation and analysis of the thin–film evolution equation from the appendix, the last of these having been successfully applied in many contexts, including in biological ones such as the growth of bacteria biofilms. From both the formal and rigorous point of view more could be done by way of weakly nonlinear analysis on the basis of the linear–stability results and plausibly exploiting bifurcatia–with–symmetry approaches. Changes of topology result from evolution beyond that demonstrated in Figure 7, for example, with the parts of the boundary that become isolated from the exterior requiring specific modelling. We record here the simplest such scenario, namely a concentric circular annulus $s_{1}(t)<r<s_{2}(t)$ containing passive material in $0<r<s_{1}(t)$ that is devoid of nutrient, so that

[TABLE]

[TABLE]

(since no kinetic–undercooling regularisation is present at $r=s_{1}(t)$, the limit $\alpha\to 0$ with $\alpha\beta=O(1)$ needs to be taken there). Hence

[TABLE]

with $s_{1},s_{2},A$ and $B$ given by the system

[TABLE]

[TABLE]

the most revealing aspect of which is its large–time behaviour

[TABLE]

which is equivalent to the corresponding one–dimensional solution that provides the dominant balance in the growing but asymptotically constant thickness, annulus. The stability or otherwise of this solution highlights more general, and challenging, open questions regarding the full bifurcation structure of solutions in order to characterise the variety of possible large–time outcomes. Finally we note that for $\beta=0$ time–reversal ($t\to-t$, $u\to-u$, $V\to-V$) of (1.1) – (1.3) leads to a meaningful problem comprising a volumetric source and surface sink of possible reference to drug treatment from the exterior of nutrient–rich tissue, though the assumption on the viscosity ratio would need reversing if the kinetic–undercooling term is to be retained.

**Acknowledgements

**The authors would like to thank the Isaac Newton Institute for Mathematical Sciences for its hospitality during the programme Coupling Geometric PDEs with Physics for Cell Morphology, Motility and Pattern Formation supported by EPSRC Grant Number EP/K032208/1. JE gratefully acknowledges the support of the EPSRC grant 1507261, VS gratefully acknowledges the support of the Leverhulme Trust Research Project Grant (RPG-2014-149) and JK and VS would like to thank the Mathematisches Forschungsinstitut Oberwolfach (MFO) for its hospitality during the programme 1704: Emerging Developments in Interfaces and Free Boundaries.

Appendix: The thin–film limit

In this appendix we apply thin–film methods to (1.1) – (1.3) in a traditional way (indeed, when the radially symmetric solution is unstable, this is perhaps the most natural next regime to study, motivation being provided by numerical results such as those in Figure 4) to derive the novel evolution equation (A.1), but supplement this derivation with a warning about its limited applicability. We consider here a tumour that is initially thin and, for brevity, two–dimensional and symmetric about the $x$–axis. Thus we take $\Gamma$ to be given in $y>0$ by

[TABLE]

with $0<\varepsilon\ll 1$, and rescale via $y=\varepsilon Y$ to give

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where to obtain a distinguished limit we have also rescaled $Q$ and $\alpha$, replacing $Q$ by $\varepsilon Q$ and $\alpha$ by $\alpha/\varepsilon$.

Integrating with respect to $Y$ from [math] to $h$ gives

[TABLE]

where the final boundary condition has yet to be used; taking $\varepsilon\rightarrow 0$ gives at leading order that

[TABLE]

[TABLE]

so an evolution equation with both a pseudoparabolic and a fourth–order regularisation, namely

[TABLE]

ensues. While interesting in its own right, (A.2) contains no destabilising terms, contrary to the conclusions above about the prevalence of instabilities. To explain this conundrum, which amounts to a warning against the naive application of thin–film approaches ((A.2) does not in and of itself provide such a warning, being a composite of well–established types), we now complement the current analysis by characterising the stability properties of the one–dimensional solution (the results being of more general local relevance for slender initial data), retaining for subsequent transparency the scalings in (A.1); $\varepsilon$ amounts to an artificial small parameter that can be removed from (A.3) by reverting the previous scalings via $\alpha\rightarrow\varepsilon\alpha$ and $h_{p}\rightarrow h_{p}/\varepsilon$.

Thus the planar solution $u=u_{p}(Y,t)$, $h=h_{p}(t)$, in which surface tension of course plays no role, is given by

[TABLE]

Perturbing in the form

[TABLE]

and retaining terms that are linear in $\delta$ only gives

[TABLE]

[TABLE]

[TABLE]

Hence for

[TABLE]

we obtain

[TABLE]

so that

[TABLE]

In the limit $\varepsilon\rightarrow 0$, $\lambda=O(1)$, this recovers the linear stability result for (A.2), namely

[TABLE]

but for short wavelengths, $\lambda=\mu/\varepsilon$, and small $\beta$, $\beta=\varepsilon^{4}b$, we obtain a different distinguished limit which is in no way foreshadowed by the thin–film analysis, namely

[TABLE]

wherein $T=\varepsilon t$, in which the first term on the right–hand side will lead to instability for sufficiently large $h_{p}$, limiting the range of applicability of the thin–film approach and reflecting a more general property of (A.3) that can also be exhibited by other distinguished limits. Thus the linear–stability analysis represents an essential adjunct to the thin–film description in characterising the evolution from slender initial data. With regard to the interpretation of the numerical results above, it is also important to note that when $Q$ in (A.2) is large the domain may in any case cease to be slender simply due to the associated growth in the $y$ direction; moreover, a thin–film analysis is not of course applicable close to the tips.