Numerical solution of the two-phase tumour growth model with moving boundary

TL;DR

This paper introduces a new numerical method for solving a one-dimensional two-phase tumour growth model that simplifies boundary handling and is validated against known solutions, showing consistency with existing results.

Contribution

A novel numerical technique for two-phase tumour growth models that avoids explicit boundary tracking and is validated against analytical solutions.

Findings

Method successfully solves tumour growth equations.

Results are consistent with existing literature.

Advantages over previous techniques demonstrated.

Abstract

A novel numerical technique has been proposed to solve a two-phase tumour growth model in one spatial dimension without needing to account for the boundary dynamics explicitly. The equivalence to the standard definition of a weak solution is proved. The method is tested against equations with analytically known solutions, to illustrate the advantages over the existing techniques. The tumour growth model is solved using the new procedure and showed to be consistent with results available in the literature.

Figures

Tables

Table 1. Table 1 : Δ ​ ℓ h Δ subscript ℓ ℎ \Delta\ell_{h} for case 1, method M.

$\Delta x$	${\alpha }_{\text{thr}}$
	$0.01$	$0.008$	$0.006$	$0.004$	$0.002$
0.01	1.67$E-3$	1.67$E-3$	1.67$E-3$	1.67$E-3$	5.00$E-3$
0.02	3.33$E-3$	3.33$E-3$	6.67$E-3$	1.33$E-2$	2.00$E-2$
0.04	6.67$E-3$	6.67$E-3$	2.00$E-2$	2.67$E-2$	4.00$E-2$
0.06	4.31$E-3$	1.58$E-2$	2.59$E-2$	4.60$E-2$	6.61$E-2$
0.08	2.10$E-2$	7.66$E-3$	1.92$E-2$	3.26$E-2$	5.93$E-2$
0.1	3.33$E-2$	1.67$E-2$	1.67$E-2$	5.00$E-2$	8.33$E-2$

Table 2. Table 2 : Δ ​ ℓ h Δ subscript ℓ ℎ \Delta\ell_{h} for case 1, method U.

	${\alpha }_{\text{thr}}$
$\Delta x$	$0.04$	$0.03$	$0.02$	$0.01$
0.01	3.33$E-3$	3.33$E-3$	1.66$E-2$	3.83$E-2$
0.02	3.33$E-2$	3.33$E-3$	1.33$E-2$	5.68$E-2$
0.04	1.20$E-1$	7.33$E-2$	6.66$E-3$	6.00$E-2$

Equations

$$\frac{\partial \alpha }{\partial t}+\frac{\partial }{\partial x}(u_c\alpha )$$

$$\frac{k{u}_c\alpha }{1-\alpha }-\mu \frac{\partial }{\partial x}\left(\alpha \frac{\partial u_c}{\partial x}\right)$$

$$\frac{\partial C}{\partial t}-\frac{{\partial }^{2}C}{\partial x^2}$$

$$\alpha (0,x)={\alpha }_0(x),C(0,x)=C_0(x)\forall x\in \Omega (0),$$

$$u_c(t,0)=0,\mu \frac{\partial u_c}{\partial x}(t,\ell (t))=\frac{\alpha (t,\ell (t))-{\alpha }_{\text{min}}}{(1-\alpha (t,\ell (t)){)}^2}H(\alpha (t,\ell (t))-{\alpha }_{\text{min}}),$$

$$\frac{\partial C}{\partial x}(t,0)=0,C(t,\ell (t))=1\forall t\in (0,T),$$

$${\ell }^{\prime }(t)=u_c(t,\ell (t))\forall t,\ell (0)={\ell }_0.$$

$${\mathbb{L}}^p\left(0,T;H^1(\Omega (t))\right)$$

$${\int }_{D_T}\left(\alpha f(\alpha ,C)\varphi +(\alpha ,u_c\alpha )\cdot {\nabla }_{t,x}\varphi \right)\mathrm{d}t\mathrm{d}x+{\int }_{\Omega (0)}\varphi (0,x){\alpha }_0(x)\mathrm{d}x=0.$$

$${\int }_{\mathcal{D}}\left(\tilde{\alpha }f(\tilde{\alpha },\tilde{C})\tilde{\varphi }+(\tilde{\alpha },{\tilde{u}}_c\tilde{\alpha }).{\nabla }_{t,x}\tilde{\varphi }\right)\mathrm{d}t\mathrm{d}x+{\int }_{\Omega (0)}\tilde{\varphi }(0,x){\alpha }_0(x)\mathrm{d}x=0.$$

$${\int }_{D_T}\left(\tilde{\alpha }f(\tilde{\alpha },\tilde{C})\tilde{\varphi }+(\tilde{\alpha },{\tilde{u}}_c\tilde{\alpha })\cdot {\nabla }_{t,x}\tilde{\varphi }\right)\mathrm{d}t\mathrm{d}x+{\int }_{\Omega (0)}\tilde{\varphi }(0,x){\alpha }_0(x)\mathrm{d}x$$

$${\int }_{\mathcal{D}\backslash D_T}(\tilde{\alpha },{\tilde{u}}_c\tilde{\alpha })\cdot {\nabla }_{t,x}\tilde{\varphi }\mathrm{d}t\mathrm{d}x+{\int }_{\mathcal{D}\backslash D_T}\tilde{\alpha }f(\tilde{\alpha },\tilde{C})\tilde{\varphi }\mathrm{d}t\mathrm{d}x$$

$${\int }_{\mathcal{D}}\left(\tilde{\alpha }f(\tilde{\alpha },\tilde{C})\tilde{\varphi }+(\tilde{\alpha },{\tilde{u}}_c\tilde{\alpha })\cdot {\nabla }_{t,x}\tilde{\varphi }\right)\mathrm{d}t\mathrm{d}x+{\int }_{\Omega (0)}\tilde{\varphi }(0,x){\alpha }_0(x)\mathrm{d}x=0.$$

$$\alpha (t,x)=\frac{(c_2-c_1){\alpha }_0(x-t)\exp ((c_1-c_2)t)}{c_1{\alpha }_0(x-t)\left(1-\exp ((c_1-c_2)t)\right)+c_2-c_1}$$

Full text

1. 1 Introduction
2. 2 Extended model
3. 3 Numerical experiments
4. 4 Conclusion

1 Introduction

We consider the tumour growth model presented in the seminal paper by Breward et al. [1]. The partial differential equations are defined in a time-dependent one-dimensional spatial domain. Such systems generally account for higher spatial dimensional models reduced to a single spatial dimension by symmetry arguments [3]. In the current model, tumour cells and surrounding fluid medium are considered as two distinct, actively interacting phases. The cell phase is viscous with viscosity $\mu$, and the fluid phase is inviscid.

Let $\Omega(t)=(0,\ell(t))$ be the interval representing the tumour where $\ell(t)$ is the tumour radius at time $t$. Set $D_{T}:=\cup_{0<t<T}\{t\}\times\Omega(t)$ and $B_{T}:=\partial D_{T}\backslash\left(\{T\}\times\Omega(T)\cup\{0\}\times\Omega(0)\cup[0,T]\times\{0\}\right)$ (Figure 1). $B_{T}$ is assumed to be of class $\mathscr{C}^{1}$ [4, p. 627]. The model seeks the variables $\alpha,u_{c}$ and $C$ that denote volume fraction of the tumour cells, velocity of the tumour cells and oxygen tension, respectively such that the following hold in $D_{T}$:

[TABLE]

The standard method to solve the system of equations of the form (1.1a)-(1.1g) is to transform the domain $\Omega(t)$ into a fixed interval using suitable change of variables [2, 6, 7]. An inverse transform is then applied to obtain the solution in the moving domain. Even though this method is commonly adopted, it comes with significant drawbacks.

Firstly, the change of variable is computable only in a few cases where the geometry of the problem is simple enough. This is even harder in 2D and 3D domains. Secondly, for the clear choice of $x\rightarrow\xi:=x/\ell(t)$ in 1D case, the discretisation error is proportional to $\ell(t)\Delta\xi$. An alternative choice is to discretise $(0,\ell(t))$ to apply numerical schemes. In this method re-meshing needs to be done at each time step which may become computationally expensive.

In this article, a new numerical technique that overcomes theses disadvantages is introduced. These are done by presenting the notion of solutions on a larger domain which contains all the time-dependent domains $\Omega(t)$ for a finite time. This domain, referred to as the extended domain, is time-independent and requires only one initial spatial discretisation; thereby avoiding the need to re-mesh. Also, the discretisation error becomes free from the dependence on $\ell(t)$.

This paper is organised as follows. In Section 2, a novel method that is referred to as extended model in the rest of the paper is introduced and its equivalence to the standard model is proved. In Section 3, the numerical technique is developed and the results are presented. The results are compared with a model problem for which analytic solutions are known. The effect of parameters in the new method is also investigated. The extended version is solved using the numerical technique developed and compared with results available from the literature. The paper ends with a conclusion in Section 4.

2 Extended model

The notion of weak solutions in the given domain and extended domain are presented in this section. The solutions for the extended and original version are proved to be equivalent. For $p\in[1,\infty]$ and a family of domains $\{\Omega(t)\}_{0\leq t\leq T}$, define

[TABLE]

Multiply (1.1a) by a test function $\phi\in\mathscr{C}_{c}^{\infty}\left(\overline{D_{T}}\backslash\left(\{T\}\times\Omega(T)\right)\right)$ and apply integration by parts. A use of (1.1g) and (1.1d) yields

[TABLE]

This constitutes the weak formulation of the hyperbolic conservation law. The weak solutions of the problem in 1.1 and in the extended model are defined next. Firstly, we give the definition of the solution in the domain $D_{T}$.

Definition 1** (Weak solution I).**

By a weak solution of the system (1.1) in $D_{T}$ we mean a 4-tuple $(\alpha,u_{c},C,\Omega)$ such that $0<\overline{m}_{\alpha}\leq\alpha\leq\overline{M}_{\alpha}<1$, $C\geq 0$ and

$\alpha\in L^{\infty}(D_{T})$* satisfies (2.1) for every $\phi\in\mathscr{C}_{c}^{\infty}\left(\overline{D_{T}}\backslash\left(\{T\}\times\Omega(T)\right)\right)$.*

• 2.

$u_{c}\in\mathbb{L}^{\infty}(0,T;H^{1}(\Omega(t)))$* with $u|_{x=0}=0$ and $C\in\mathbb{L}^{2}(0,T;H^{1}(\Omega(t)))$ with $C|_{x=\ell(t)}=1$ are solutions of (1.1b) and (1.1c) in the sense of distributions.*

• 3.

The domain $\Omega(t)$ is the open interval $(0,\ell(t))$ where $\ell(t)$ is governed by (1.1g).

The definition of the weak solution in the domain $\mathscr{D}$ is given next. $\mathscr{D}$ is the extended domain given by $(0,T)\times(0,L)$ (Figure 1) where $L$ is chosen such that $\ell(t)<L$ for every $t\leq T$.

Definition 2** (Weak solution II).**

By a weak solution of the system (1.1) in the extended domain $\mathscr{D}$ (Figure 1) we mean a 4-tuple $(\tilde{\alpha},\tilde{u}_{c},\tilde{C},\tilde{\Omega})$ such that $0<\overline{m}_{\alpha}\leq\tilde{\alpha}|_{\tilde{\Omega}}\leq\overline{M}_{\alpha}<1$, $\tilde{C}\geq 0$ and

$\tilde{\alpha}\in L^{\infty}(\mathscr{D})$* such that for every $\tilde{\phi}\in\mathscr{C}_{c}^{\infty}([0,T)\times(0,L))$*

[TABLE]

• 2.

For a fixed $t$, $\tilde{\Omega}(t):=\{x:\tilde{\alpha}(t,x)>0\}$, $\tilde{u}_{c}=0$, $\tilde{C}=1$ on $(0,L)\backslash\tilde{\Omega}(t)$. Define $\tilde{D}_{T}:=\cup_{0<t<T}\{t\}\times\tilde{\Omega}(t)$.

• 3.

$\tilde{u}_{c}\in L^{\infty}(\mathscr{D})$* with $u_{c}:=\tilde{u}_{c}|_{\tilde{D}_{T}}$, $u_{c}\in\mathbb{L}^{\infty}(0,T;H^{1}(\tilde{\Omega}(t)))$ and, $\tilde{C}\in L^{2}(\mathscr{D})$ such that $C:=\tilde{C}|_{\tilde{D}_{T}}\in\mathbb{L}^{2}(0,T;H^{1}(\tilde{\Omega}(t)))$ are solutions of (1.1b) and (1.1c) in the sense of distributions.*

Theorem 3**.**

If $(\alpha,u_{c},C,\Omega)$ is a weak solution I, then $(\tilde{\alpha},\tilde{u}_{c},\tilde{C},\tilde{\Omega})$ defined by $\tilde{\alpha}:=\alpha,\tilde{u}_{c}:=u_{c}$ and $\tilde{C}:=C$ in $D_{T}$ and, $\tilde{\alpha}:=0,\tilde{u}_{c}:=0,\tilde{C}:=1$ in $\mathscr{D}\backslash D_{T}$ with $\tilde{\Omega}(t):=\Omega(t)$ is weak solution II. Conversely, if $(\tilde{\alpha},\tilde{u}_{c},\tilde{C},\tilde{\Omega})$ is a weak solution II, then $(\alpha,u_{c},C,\Omega)$ with $\Omega=\tilde{\Omega}$ and $\alpha:=\tilde{\alpha}|_{\tilde{D}_{T}},u_{c}:=\tilde{u}_{c}|_{\tilde{D}_{T}}$ and $C:=\tilde{C}|_{\tilde{D}_{T}}$ is a weak solution I.

Proof:

Let $(\alpha,u_{c},C,\Omega)$ be a weak solution I and $\tilde{\phi}\in\mathscr{C}_{c}^{\infty}([0,T)\times(0,L))$. Since $\tilde{\phi}|_{D_{T}}\in\mathscr{C}_{c}^{\infty}\left(\overline{D_{T}}\backslash\left(\{T\}\times\Omega(T)\right)\right)$,  (2.1) holds true. Let $\tilde{\alpha}=\alpha$ in $D_{T}$ and $\tilde{\alpha}=0$ in $\mathscr{D}\backslash D_{T}$. Then a use of the definitions of $\tilde{u}_{c}$ and $\tilde{C}$ yields

[TABLE]

Add (2.3) and (2.4) to obtain

[TABLE]

Therefore (2.2) holds true. The conditions on $\tilde{u}_{c}$ and $\tilde{C}$ follow naturally from the definition 2. Since $\tilde{\alpha}>0$ in $D_{T}$ and $\tilde{\alpha}=0$ in $\mathscr{D}\backslash D_{T}$, $\tilde{\Omega}(t)=\Omega(t)$ for every $t\in[0,T)$. Therefore $(\tilde{\alpha},\tilde{u}_{c},\tilde{C},\tilde{\Omega})$ is a weak solution II.

Conversely, assume that $(\tilde{\alpha},\tilde{u}_{c},\tilde{C},\tilde{\Omega})$ is a weak solution II. Let $\phi\in\mathscr{C}_{c}^{\infty}\left(\overline{D_{T}}\backslash\left(\{T\}\times\Omega(T)\right)\right)$. Define $\tilde{\phi}\in\mathscr{C}_{c}^{\infty}([0,T)\times(0,L))$ such that $\tilde{\phi}=\phi$ in $D_{T}$. Since $\Omega(t)=\tilde{\Omega}(t)$ for every $t$, $\tilde{\alpha}=0$ in $\mathscr{D}\backslash D_{T}$. Using this in (2.2) we obtain (2.1). We shall recover (1.1g) next. For this define a vector field $\mathbf{F}:\mathscr{D}\rightarrow\mathbb{R}^{2}$ by $\mathbf{F}(t,x):=(\tilde{\alpha},\tilde{u}_{c}\tilde{\alpha}).$ We set ${\bf F}|_{B_{T}^{+}}=({\bf F}|_{D_{T}})|_{B_{T}}$ and ${\bf F}|_{B_{T}^{-}}=({\bf F}|_{\mathscr{D}\backslash D_{T}})|_{B_{T}}$. Since the weak divergence of the vector field ${\bf F}$ is -$\tilde{\alpha}f(\tilde{\alpha},\tilde{C})\in L^{2}(\mathscr{D})$, the flux of ${\bf F}$ is continuous across $B_{T}$. Since $\tilde{\alpha}=0$ in $\mathscr{D}\backslash D_{T}$, ${\bf F}|_{B_{T}^{-}}={\bf 0}$. Therefore, $({\bf F}|_{B_{T}^{+}}-{\bf F}|_{B_{T}^{-}})\cdot{\bf n}_{B_{T}}=(\alpha,u_{c}\alpha)\cdot{\bf n}_{B_{T}}=0$ where ${\bf n}_{B_{T}}$ is the normal to $B_{T}$ given by $\left(\left|\ell^{\prime}(t)\right|^{2}+1\right)^{-1/2}\left(-\ell^{\prime}(t),1\right)$. This gives $(\alpha,u_{c}\alpha)\cdot{\bf n}_{B_{T}}=0$. Since $\alpha>0$, $\ell^{\prime}(t)=u_{c}(t,\ell(t))$. The conditions on $u_{c}$ and $C$ follows directly from the definitions. Therefore $(\alpha,u_{c},C,\Omega)$ is a weak solution I.

This completes the proof of the equivalence between the solutions.   $\spadesuit$

3 Numerical experiments

By Theorem 3 it is enough to solve (1.1a) in the extended domain $(0,L)$. Equation (1.1a) is solved using cell-centred finite volume methods. In particular, we use two methods to solve the volume fraction equation: upwinding with Godunov flux [5, p. 135] (method U), and MUSCL with Godunov flux [5, p. 146] (method M). The uniform space and time discretisations are $0=x_{0}<x_{1}<\cdots<x_{i}<\cdots<x_{M}=L$, $0=t_{0}<t_{2}<\cdots<t_{j}<\cdots<t_{N}=T$ with $h=x_{i+1}-x_{i}$ and $\Delta t=t_{j+1}-t_{j}$. The right hand side boundary $\ell(t)$ is approximated by $\ell_{h}(t)=\min_{x}\{x:\tilde{\alpha}<\alpha_{\text{thr}}\;\text{ on }(x,L)\}$ where $\alpha_{\text{thr}}$ is a small positive number. Define $\tilde{\alpha}:=0$ for $x\geq\ell_{h}(t)$ to eliminate the error caused by small positive values of $\tilde{\alpha}$ (created by numerical diffusion) in there. Equations (1.1b) and (1.1c) are solved using conforming $P_{1}$ finite element method (FEM) in space and forward finite difference in time, in the reconstructed domain $(0,\ell_{h}(t))$. This procedure, referred to as scheme A in the rest of the paper, is outlined below.

Start at $t_{0}=0$. Solve $\tilde{u}_{ch}^{0}$ using $\tilde{\alpha}_{h}^{0}$ (initial condition).

For $j=1$ to $N$, $t_{j}=t_{j-1}+\Delta t$.

$\ell_{h}^{j}=\displaystyle\min\limits_{x_{i}}\{x_{i}:\tilde{\alpha}_{h}^{j-1}<\alpha_{\text{thr}}\;\text{ on }(x_{i},L)\}$.

• 3.

Find $\tilde{u}_{ch}^{j}\;\&\;\tilde{C}_{h}^{j}$ in $(0,\ell_{h}^{j})$ ($P_{1}$ conforming FEM).

• 4.

Extrapolate $\tilde{u}_{ch}^{j}=0\;\&\;\tilde{C}_{h}^{j}=1$ to $(\ell_{h}^{j},L)$.

• 5.

Find $\tilde{\alpha}_{h}^{j}$ in $(0,L)$ (method U or M).

The complete elimination of re-meshing and applicability in higher dimensions are the major advantages of this scheme. Scheme B denotes the procedure of obtaining numerical solution in the scaled domain $(0,1)$ [1]. Two test cases are considered in the numerical experiments.

In the first case, the cell velocity $u_{c}$ and the oxygen tension $C$ are assumed to be unity. In this case (1.1a) reduces to a semi-linear advection equation which can be solved analytically by the method of characteristics. The analytical solution is compared with the numerical solutions. We also study the influence of $\alpha_{\text{thr}}$ on locating the tumour frontier. In the second case, we compare the approximate solutions of the full system (1.1a)-(1.1g).

In all numerical tests the values of the parameters are set to be $s_{1}=10=s_{4},\,s_{2}=0.5=s_{3}$, $k=1=\mu$, $Q=0.5$, $\hat{Q}_{1}=0$ and $\ell(0)=1$ to preserve conformity with Breward et al. [1].

Case 1

The analytical solution to (1.1a) in the case where $C=u_{c}=1$ is:

[TABLE]

where $c_{1}=1$ and $c_{2}=\frac{s_{2}+s_{3}}{1+s_{4}}$. The initial data considered are (i) $\alpha_{0}(t,x)=0.5\left(0.02+\cos^{2}\left(x\right)\right)\chi_{[0,1]}$ (ii) $\alpha_{0}(t,x)=0.5\left(0.02+\sin^{2}\left(x\right)\right)\chi_{[0,1]}$ and (iii) $\alpha_{0}(t,x)=\frac{\chi_{[0,1]}}{2}\frac{1+\exp(x-0.5)^{2}}{1+\exp(2(x-0.5)^{2})}$, where $\chi_{[0,1]}=1$ in $[0,1]$ and [math] otherwise. Here $T=5,\,L=6,\,\Delta t=0.01,\,\Delta x=0.02,\,\alpha_{\text{thr}}=0.04$ (method U) and $\alpha_{\text{thr}}=0.004$ (method M). The reduction in the numerical diffusion from former to latter method explains the reduction in the threshold value.

Approximate solution obtained in the extended domain captures the properties of the analytical solution better than the one obtained in the scaled domain (Figure 2), though it is less accurate towards the discontinuity at $\ell_{h}^{j}$ in method U owing to high diffusion. But method M overcomes this disadvantage; the extended solution agrees well with the scaled solution towards the discontinuity, and remarkably better in the interior region (Figure 2). The recovered radius, on the other hand, is in excellent agreement with the exact radius for both method M and U with a proper choice of the threshold value (Table 1).

We conclude this section by analysing the dependency of the recovered radius on the threshold value $\alpha_{\text{thr}}$ for the MUSCL method. The relative error, $\Delta\ell_{h}=\frac{|\ell(T)-\ell_{h}^{N}|}{\ell(T)}$, at $T=5$ is used as a quantification of the error in the recovered radius. Two sets of experiments are conducted; (a) varying $\alpha_{\text{thr}}$ at a fixed $\Delta x$ (b) varying $\Delta x$ at a fixed $\alpha_{\text{thr}}$. Table 1 shows that there exists a wide range of $\alpha_{\text{thr}}$ and $\Delta x$ for which the error remains below 1$\operatorname{E}-2$. This assures the accuracy of the method while the selection of $\alpha_{\text{thr}}$ remains a pertinent problem.

The range of $\alpha_{\text{thr}}$ and $\Delta x$ for the error remains low is thin for method U, which is expected considering the high numerical diffusion associated with it.

Case 2

The two phase model with all the system variables treated as unknowns is considered in this case. The parameters are chosen as $\Delta t=0.01$, $\Delta x=0.01$, $T=228$, $L=25$, $\alpha_{\text{thr}}=0.004$ (method M) and $\alpha_{\text{thr}}=0.01$ (method U) based on Tables 1, 2. The initial condition is $\alpha_{0}(x)=0.8$ for $0\leq x\leq 1$ and [math] otherwise. The moving boundary $\ell(t)$ is well captured by methods U and M. Since the exact value of $\ell(T)$ is not available, the error is quantified as the relative difference between scheme B and scheme A. The difference for method U is $6.18\operatorname{E}-3$ and method M is $5.69\operatorname{E}-3$. The numerical solution in the extended domain is in good agreement with the solution obtained from the scaled domain [1].

4 Conclusion

A novel numerical technique is developed to solve the two phase tumour growth problem and is tested against problems for which analytical solutions are known. For a fixed spatial mesh size the new method gives better solution than the standard method of solving in a scaled domain. The moving boundary is recovered from the numerical solution by comparing with a threshold value. It is found that an appreciable range of threshold values can be used along with higher order methods like MUSCL so that the error in the recovered radius can be kept low. The solution obtained by applying this technique shows very good agreement with solutions obtained using standard methods. This emphasises the reliability of the new method in extending it to solve tumour growth problems in higher dimensions while not solving for the boundary explicitly.

Acknowledgement

The author expresses gratitude towards A/Prof. Jérôme Droniou (Monash Univeristy), Dr. Jennifer Anne Flegg (Melbourne University) and Prof. Neela Nataraj (I.I.T. Bombay) for the valuable suggestions and help.

Author address