Numerical solution of the twophase tumour growth model with moving boundary
Abstract
A novel numerical technique has been proposed to solve a twophase 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.
IITBMonash Research Academy, IIT Bombay, Powai, Mumbai 400076, INDIA \mailtogopikrishnan.chirappurathuremesan@monash.edu \myorcid0000000245074463
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 timedependent onedimensional 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 , and the fluid phase is inviscid.
Let be the interval representing the tumour where is the tumour radius at time . Set and (Figure 1). is assumed to be of class [4, p. 627]. The model seeks the variables and that denote volume fraction of the tumour cells, velocity of the tumour cells and oxygen tension, respectively such that the following hold in :
(1.1a)  
(1.1b)  
(1.1c)  
where . The positive constants (drag coefficient) controls the drag between the phases; and , the stress in the cellular phase; , , and , birth and death rates. The Heaviside function if and zero otherwise. The twophase model in Breward et al. [1] uses a quasisteady state assumption for oxygen tension which is relaxed in this study. This means the explicit temporal variation of oxygen tension is considered which makes it a parabolic equation (1.1c). The initial and boundary conditions are  
(1.1d)  
(1.1e)  
(1.1f)  
(1.1g)  
Here, satisfies for every . 
The standard method to solve the system of equations of the form (1.1a)(1.1g) is to transform the domain 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 in 1D case, the discretisation error is proportional to . An alternative choice is to discretise to apply numerical schemes. In this method remeshing 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 timedependent domains for a finite time. This domain, referred to as the extended domain, is timeindependent and requires only one initial spatial discretisation; thereby avoiding the need to remesh. Also, the discretisation error becomes free from the dependence on .
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 and a family of domains , define
Multiply (1.1a) by a test function and apply integration by parts. A use of (1.1g) and (1.1d) yields
(2.1) 
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 .
[Weak solution I] By a weak solution of the system (1.1) in we mean a 4tuple such that , and
The definition of the weak solution in the domain is given next. is the extended domain given by (Figure 1) where is chosen such that for every .
[Weak solution II] By a weak solution of the system (1.1) in the extended domain (Figure 1) we mean a 4tuple such that , and

such that for every
(2.2) 
For a fixed , , , on . Define .
If is a weak solution I, then defined by and in and, in with is weak solution II. Conversely, if is a weak solution II, then with and and is a weak solution I.
Proof.
Let be a weak solution I and . Since , (2.1) holds true. Let in and in . Then a use of the definitions of and yields
(2.3)  
(2.4) 
(2.5) 
Therefore (2.2) holds true. The conditions on and follow naturally from the definition 2. Since in and in , for every . Therefore is a weak solution II.
Conversely, assume that is a weak solution II. Let . Define such that in . Since for every , in . Using this in (2.2) we obtain (2.1). We shall recover (1.1g) next. For this define a vector field by We set and . Since the weak divergence of the vector field is , the flux of is continuous across . Since in , . Therefore, where is the normal to given by . This gives . Since , . The conditions on and follows directly from the definitions. Therefore is a weak solution I.
This completes the proof of the equivalence between the solutions. ∎
3 Numerical experiments
By Theorem 2 it is enough to solve (1.1a) in the extended domain . Equation (1.1a) is solved using cellcentred 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 , with and . The right hand side boundary is approximated by where is a small positive number. Define for to eliminate the error caused by small positive values of (created by numerical diffusion) in there. Equations (1.1b) and (1.1c) are solved using conforming finite element method (FEM) in space and forward finite difference in time, in the reconstructed domain . This procedure, referred to as scheme A in the rest of the paper, is outlined below.

Start at . Solve using (initial condition).
For to , .

.

Find in ( conforming FEM).

Extrapolate to .

Find in (method U or M).
The complete elimination of remeshing 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 [1]. Two test cases are considered in the numerical experiments.
In the first case, the cell velocity and the oxygen tension are assumed to be unity. In this case (1.1a) reduces to a semilinear 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 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 , , , and to preserve conformity with Breward et al. [1].
Case 1
The analytical solution to (1.1a) in the case where is:
(3.1) 
where and . The initial data considered are (i) (ii) and (iii) , where in and otherwise. Here (method U) and (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 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 for the MUSCL method. The relative error, , at is used as a quantification of the error in the recovered radius. Two sets of experiments are conducted; (a) varying at a fixed (b) varying at a fixed . Table 1 shows that there exists a wide range of and for which the error remains below 1. This assures the accuracy of the method while the selection of remains a pertinent problem.
0.01  1.67  1.67  1.67  1.67  5.00 

0.02  3.33  3.33  6.67  1.33  2.00 
0.04  6.67  6.67  2.00  2.67  4.00 
0.06  4.31  1.58  2.59  4.60  6.61 
0.08  2.10  7.66  1.92  3.26  5.93 
0.1  3.33  1.67  1.67  5.00  8.33 
0.01  3.33  3.33  1.66  3.83 

0.02  3.33  3.33  1.33  5.68 
0.04  1.20  7.33  6.66  6.00 
The range of and 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 , , , , (method M) and (method U) based on Tables 1, 2. The initial condition is for and otherwise. The moving boundary is well captured by methods U and M. Since the exact value of is not available, the error is quantified as the relative difference between scheme B and scheme A. The difference for method U is and method M is . 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.
References
 [1] C. J. W. Breward, H. M. Byrne, and C. E. Lewis. The role of cellcell interactions in a twophase model for avascular tumour growth. Journal of Mathematical Biology, 45(2):125–152, 2002.
 [2] C. J. W. Breward, H. M. Byrne, and C. E. Lewis. A multiphase model describing vascular tumour growth. Bulletin of Mathematical Biology, 65(4):609–640, 2003.
 [3] H. M. Byrne, J. R. King, D. L. Sean McElwain, and L. Preziosi. A twophase model of solid tumour growth. Appl. Math. Lett., 16:567–573, 2003.
 [4] L. Evans. Partial Differential Equations (Graduate Studies in Mathematics, V. 19) GSM/19. American Mathematical Society, 1998.
 [5] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In Solution of Equation in ℝn (Part 3), Techniques of Scientific Computing (Part 3), volume 7 of Handbook of Numerical Analysis, pages 713 – 1018. Elsevier, 2000.
 [6] J. Ward and J. R. King. Mathematical modelling of avasculartumour growth. IMA Journal of Mathematics Applied in Medicine and Biology, 14:39–69, 04 1997.
 [7] J. Ward and J. R. King. Mathematical modelling of avasculartumor growth ii: Modelling growth saturation. IMA Journal of Mathematics Applied in Medicine and Biology, 16:171–211, 06 1999.