On Finite Volume Discretization of Infiltration Dynamics in Tumor Growth Models
OON FINITE VOLUME DISCRETIZATION OF INFILTRATION DYNAMICS INTUMOR GROWTH MODELS
XIANYI ZENG, MASHRIQ AHMED SALEH, AND JIANJUN PAUL TIANA bstract . We address numerical challenges in solving hyperbolic free boundary problemsdescribed by spherically symmetric conservation laws that arise in the modeling of tumorgrowth due to immune cell infiltrations. In this work, we normalize the radial coordinateto transform the free boundary problem to a fixed boundary one, and utilize finite volumemethods to discretize the resulting equations. We show that the conventional finite volumemethods fail to preserve constant solutions and the incompressibility condition, and theytypically lead to inaccurate, if not wrong, solutions even for very simple tests. Theseissues are addressed in a new finite volume framework with segregated flux computationsthat satisfy su ffi cient conditions for ensuring the so-called totality conservation law and thegeometric conservation law. Classical first-order and second-order finite volume methodsare enhanced in this framework. Their performance is assessed by various benchmark teststo show that the enhanced methods are able to preserve the incompressibility constraintand produce much more accurate results than the conventional ones.
1. I ntroduction
Modeling the tumor growth due to immune cell infiltration using partial di ff erentialequations (PDE) has been an active research area in recent years. One of the earliest papersaddressing this phenomenon from a mathematical point of view is by Evelyn F. Keller andLee A. Segel [17], who model the cell movements by Brownian motion and concludethat they generally move towards a region with high chemoattractant concentration. ThePatlak-Keller-Segel (PKS) chemotaxis system, which describes the interaction between thecell and the chemoattractant, is then studied both theoretically and numerically by variousauthors [1, 3, 11, 23, 18, 9, 6]. Existing literature focuses on solving the PKS system on afixed domain; hence they are suitable for describing the cell movements inside the tumorbut not for modeling how the tumor grows. Recently, B. Niu and the authors of the currentpaper propose a free boundary model that extends the PKS system to describe the growthof tumor due to immune cell infiltration [21]. In this model, the immune cells are attractedby the chemoattractant that usually has higher concentration inside the tumor and enterthe tumor boundary; the mean cell movement velocity is derived by assuming the cellsare incompressible, i.e., the total cell number per unit volume is assumed to be constant.The incompressibility is a crucial assumption – because the cells have fixed volume, whenimmune cells enter through the tumor boundary they need to compete with native ones forspace and eventually promote tumor growth.It should be noted that treating biological systems as free boundary problems is by nomeans new. In the literature, there are numerous successful studies addressing the existence Date : May 9, 2019.2010
Mathematics Subject Classification.
Key words and phrases.
Finite volume methods; Cell incompressibility; Free boundary problems; Patlak-Keller-Segel system; Tumor growth modeling. a r X i v : . [ m a t h . NA ] S e p X. ZENG, M. SALEH, AND J. TIAN and uniqueness of solutions to such PDE systems [5, 8] as well as conducting well-behavednumerical simulations [16]. We would like to emphasize, however, that these studies relyon the fact that the same velocity field is used for the advection of all cell species; hencea characteristic method (in the analytical approach) or a Lagrangian strategy (in the com-putational approach) can be applied. This is not the case with infiltration dynamics, asby nature the invading species and the native ones are carried by di ff erent velocities. Itis worth mentioning that in a recent work by A. Friedman et al. [13], the authors provethe global existence and uniqueness of solutions to a free boundary problem that containsinfiltrating species; however, the governing equations therein are of parabolic type, whichis very di ff erent from what we’re considering here – because there is no di ff usion term forcell species, the model considered in this paper does not contain regularization and shocksdo occur in the solution process.Studies on free boundary problems of hyperbolic type that involves infiltration dynam-ics, to our best knowledge, remain scarce in the literature. In this work, we attempt toclose this gap by proposing a new finite volume framework for the discretization of a gen-eral class of equations. Particularly, we consider the migration of two categories of cellspecies – the cell species belonging to the first category move inside the tumor and willnever cross the boundary, whereas the second category involve all the infiltrating species;the motion of both types are governed by hyperbolic equations. The methodology is de-scribed in a very general setting, in the sense that it is not restricted to any particular cellproliferation, apoptosis, and interaction models. To this end, the method we propose issuitable for the investigation of any similar systems, such as the plaque development andthe wound healing processes [13, 15, 12]. However, for the ease of statement we set ourcontext in tumor growth modeling and use the term “cell” to refer to any entities that playa part in the incompressibility constraint, see Section 2.In previous work [21], the spherically symmetric free-boundary problem is consideredand solved numerically by first mapping the physical coordinates onto a fixed logical oneand then discretization using the conventional finite volume methods, see also Section 2 fora brief review of this model. Although the shocks are captured nicely, clear violation of theincompressibility assumption is observed, especially near the tumor center. A major causeis that incompressibility is not enforced directly by the model; instead, it is assumed in thederivation of the velocity equation. In addition, geometrical source terms appear when wechange from the physical coordinate to the logical one; and existing finite volume methodscannot balance them well, even when the solutions are constants.To resolve these issues, we investigate a simplified model that easily extends to thefull tumor growth model of [21]. The totality conservation law (TCL) and the geometricconservation law (GCL) are defined and justified as the criterion for any numerical methodto maintain constant solutions and satisfy the incompressibility condition. The new finitevolume methods are developed in three stages. First, we design a general finite volumeframework for solving the model system, and extend the TCL and GCL to the discrete level,called DTCL and DGCL, respectively, where the first letter “D” stands for “discrete”. Next,we propose several consistency properties, so that for any numerical flux that satisfies theseproperties, the resulting method will satisfy both DTCL and DGCL. Finally, the classicalfirst-order upwind method and the second-order MUSCL flux [25] are enhanced accordingto these conditions.The remainder of the paper is organized as follows. In Section 2 we briefly review theoriginal tumor growth model as well as the incorporation of the incompressibility assump-tion. Then, a simplified model that captures the most important features is described in VM OF INFILTRATION IN TUMOR GROWTH 3
Section 3. The main results and the proposed methods are derived in Section 4, where wepropose the DTCL and DGCL conditions and prove su ffi cient conditions for the numericalmethod to satisfy these conditions. Extensive numerical tests are provided in Section 5 toassess the performance of the enhanced methods, which is compared to the existing finitevolume methods. Finally, Section 6 concludes this paper.2. A R eview of the T umor G rowth M odel and I ts F inite V olume D iscretization In the tumor growth model proposed earlier [21], we consider the movement of glioma(or cancer) cells, necrotic cells, and immune cells, whose number densities are denotedby G ( r , t ), N ( r , t ), and M ( r , t ), respectively. Here r is the distance from a point insidethe (spherically symmetric) tumor to the center and t is the time ordinate. The cells aresupposed to be incompressible, in the sense that one expects:(2.1) G ( r , t ) + N ( r , t ) + M ( r , t ) = θ , for some constant θ that designates the total number of cells per unit volume.The velocities of the cell movements are determined by two aspects. First, because ofthe incompressibility assumption each cell takes a fixed volume; hence when the cells aresqueezed they tend to move to the nearby region and eventually cause the tumor to growor shrink. The velocity due to the cell-volume-preserving mechanism is denoted by V ( r , t ),and it is the solely velocity that is responsible for the movement of glioma cells and necroticcells. Second, in addition to V , the immune cells are also guided by the chemoattractantconcentration, as discussed by Evelyn F. Keller and Lee A. Segal [17] in the early 1970s.The corresponding velocity is denoted by U ( r , t ), and it is positive related to the gradient ∂ A ( r , t ) /∂ r , where A ( r , t ) is the chemoattractant concentration.In the spherical coordinates, the equations that govern the cell movements are thuslygiven by: ∂ G ∂ t + r ∂∂ r (cid:104) r GV (cid:105) = f ( r , t , G , N , M ) , (2.2a) ∂ N ∂ t + r ∂∂ r (cid:104) r NV (cid:105) = g ( r , t , G , N , M ) , (2.2b) ∂ M ∂ t + r ∂∂ r (cid:104) r M ( V + U ) (cid:105) = h ( r , t , G , N , M ) , (2.2c) U = α ∂ A ∂ r . (2.2d)Here α > f , g , h are sourceterms that describe the production and diminishing of the cells. In this paper, we follow theconvention that a single upper case letter denotes a dependent variable to be solved, and asingle lower case letter designates an independent variable or a prescribed function. Ournumerical method will not depend on the particular forms of the source functions; from amodeling point of view, however, examples of these functions are given below. Let λ and µ be the self-production and transformation rates of the cancer cells, we have:(2.2e) f ( G , N , M ) = λ G − µ G ;here µ is the rate at which the cancer cells convert to necrotic cells, which are removedfrom the tumor by the rate δ , hence one can model:(2.2f) g ( G , N , M ) = µ G − δ N ; X. ZENG, M. SALEH, AND J. TIAN and finally if the only way that the immune cells are gone is through their own death, whichhappens at the rate ρ , then the source term h can be modeled as:(2.2g) h ( G , N , M ) = − ρ M . For more details about the rationale behind these source functions, the readers are referredto [21] and the references therein.The equations (2.2a)–(2.2d) are valid for all ( r , t ) : 0 ≤ r ≤ R ( t ), where R ( t ) > R (cid:48) ( t ) = V ( R ( t ) , t ) . The equation for the velocity field V ( r , t ) is derived by summing up (2.2a)–(2.2c) andinvoking the incompressibility assumption (2.1):(2.2i) 1 r ∂∂ r (cid:104) r θ V + r U M (cid:105) = f + g + h . To complete the system, the chemoattractant A is generally secreted by the glioma cellsand subject to the di ff usion rate ν and diminishing rate γ :(2.2j) ∂ A ∂ t = ν r ∂∂ r (cid:34) r ∂ A ∂ r (cid:35) + χ mG β + G − γ A , ≤ r < + ∞ . Note that this equation is valid on the entire domain since the chemoattractant exists inthe entire body, which is supposed to be much larger than the tumor. The indicator func-tion χ in the second term of the right hand side equals 1 when 0 ≤ r ≤ R ( t ) and equals 0otherwisely.Finally, the governing equation (2.2) is complemented by appropriate initial conditionsfor G , N , M , and A such that (2.1) is satisfied, and the following boundary conditions: ∂ G (0 , t ) ∂ r = ∂ N (0 , t ) ∂ r = , (2.3a) ∂ M (0 , t ) ∂ r = , M ( R ( t ) , t ) = M bc ( t ) if U ( R ( t ) , t ) < , (2.3b) ∂ A (0 , t ) ∂ r = , lim r → + ∞ A ( r , t ) = , (2.3c) V (0 , t ) = . (2.3d)Here the second part of (2.3b) is known as the incoming boundary condition and M bc isthe prescribed embient number density of immune cells.2.1. Conservation form in normalized coordinates.
To avoid the di ffi culty of dealingwith a time-varying domain, we cast the equations to the normalized coordinates ( η, τ ) = VM OF INFILTRATION IN TUMOR GROWTH 5 ( R ( t ) / t , t ) and rescale the equations to obtain a conservation system: ∂ ( η R G ) ∂τ + ∂∂η (cid:34)(cid:32) VR − η R (cid:48) R (cid:33) η R G (cid:35) = η R f ( G , N , M ) − η R (cid:48) RG , (2.4a) ∂ ( η R N ) ∂τ + ∂∂η (cid:34)(cid:32) VR − η R (cid:48) R (cid:33) η R N (cid:35) = η R g ( G , N , M ) − η R (cid:48) RN , (2.4b) ∂ ( η R M ) ∂τ + ∂∂η (cid:34)(cid:32) VR − η R (cid:48) R + UR (cid:33) η R N (cid:35) = η R h ( G , N , M ) − η R (cid:48) RM , (2.4c) U = α R ∂ A ∂η , (2.4d) 1 η R ∂∂η (cid:104) η ( θ V + U M ) (cid:105) = f + g + h , (2.4e)for all 0 ≤ η ≤ τ ≥
0; and(2.4f) ∂ ( η R A ) ∂τ + ∂∂η (cid:34)(cid:32) − η R (cid:48) R (cid:33) η R A (cid:35) = ν ∂∂η (cid:32) η ∂ A ∂η (cid:33) + χ m η R G β + G − γη R A − η R (cid:48) RA , on the domain 0 ≤ η < + ∞ , τ ≥
0, and χ is the indicator function that equals 1 when0 ≤ η ≤ R (cid:48) ( τ ) = V (1 , τ ) . Note that the − η R (cid:48) R terms are new, and they appear because of our change of coordinates.The boundary conditions are given by: ∂ G (0 , τ ) ∂η = ∂ N (0 , τ ) ∂η = , (2.5a) ∂ M (0 , τ ) ∂η = , M (1 , τ ) = M bc ( τ ) if U (1 , τ ) < , (2.5b) ∂ A (0 , τ ) ∂η = , lim η → + ∞ A ( η, τ ) = , (2.5c) V (0 , τ ) = . (2.5d)2.2. Finite volume discretization.
In previous work, the conservative equations (2.4a)–(2.4c) and (2.4f) are discretized by the standard finite volume methods, see for exam-ple [25, 19]. We briefly review the first-order upwind method here as well as introducesome notations that will be used throughout the paper.The logical domain η ∈ [0 ,
1] is divided into N η uniform intervals , each of which haslength ∆ η = / N η ; and we denote the interval faces by η j = j ∆ η and interval centers by η j − / = ( j − / ∆ η . For easy reading, we use the integer subscripts to denote nodal vari-ables, whereas the half-integer subscripts to denote the variables that are associated withintervals, such as the interval-averages.In particular, because the cell numbers are conserved quantities, in the general finitevolume discretization these variables are defined for each interval, and they’re denoted by G j − / , N j − / , and M j − / , where 1 ≤ j ≤ N η . Considering in addition the forward-Euler To avoid confusion, we reserve the word “cell” exclusively for denoting the cell species; whereas the com-monly used “cell” in finite volume discretization is referred to as “interval” throughout the paper.
X. ZENG, M. SALEH, AND J. TIAN time integrator and designating the discrete solutions at time step τ n by the superscript n ,the general finite volume discretization reads:( R n + ) G n + j − / − ( R n ) G nj − / ∆ τ + F G , nj − F G , nj − η j − / ∆ η = ( R n ) f nj − / − R (cid:48) n R n G nj − / , (2.6a) ( R n + ) N n + j − / − ( R n ) N nj − / ∆ τ + F N , nj − F N , nj − η j − / ∆ η = ( R n ) g nj − / − R (cid:48) n R n N nj − / , (2.6b) ( R n + ) M n + j − / − ( R n ) M nj − / ∆ τ + F M , nj − F M , nj − η j − / ∆ η = ( R n ) h nj − / − R (cid:48) n R n M nj − / , (2.6c)where f nj − / = f ( G nj − / , N nj − / , M nj − / ) and h nj − / and g nj − / are similarly computed; theradius related quantities are:(2.6d) R (cid:48) n = V nN η , R n + = R n + ∆ τ V nN η . We define the velocity at the nodes, and V nN η is the numerical approaximation to V (1 , τ n ),see also the discussion below (2.8).The numerical flux F X , nj , where X stands for G , N , or M , is an approximation to thecorresponding flux for X at η j . If we apply the existing finite volume methods to computethese numerical fluxes, for example, by using the first-order upwind flux, we have: F G , nj = F upw V nj R n − η j R (cid:48) n R n ; η j − / ( R n ) G nj − / , η j + / ( R n ) G nj + / , (2.7a) F N , nj = F upw V nj R n − η j R (cid:48) n R n ; η j − / ( R n ) N nj − / , η j + / ( R n ) N nj + / , (2.7b) F M , nj = F upw V nj R n − η j R (cid:48) n R n + U nj R n ; η j − / ( R n ) M nj − / , η j + / ( R n ) M nj + / . (2.7c)Here the upwind flux is defined as:(2.8) F upw ( W ; X l , X r ) = WX l , if W ≥ WX r , if W < , where the subscripts l and r mean “left” and “right”, respectively, and W is the local ad-vection velocity at the interval face between the two interval values X l and X r .In (2.7), the velocity variables V nj and U nj , where 0 ≤ j ≤ N η , are collocated at the intervalface η j . Here the velocity V is computed using the integral form of (2.4e): V n = , (2.9a) V nj = θη j j (cid:88) k = η k − / R n ( f nk − / + g nk − / + h nk − / ) − θ U nj M nj , ≤ j ≤ N η , (2.9b)where M j is computed as the mean of surrounding interval-averaged values: M j = ( M nj − / + M nj + / ) /
2, except for the last node, in which case M N η = M N η − / . As for the velocity U , we have U n = U nj = α ( A nj + / − A nj − / ) / ( ∆ η R n ), where A nj − / is the averagedchemoattractant concentration on [ η j − , η j ]. The chemoattractant concentration is com-puted by approximating the convective term (2.4f) by straightforward finite volume dis-cretization and the di ff usion term by central di ff erence approximation. Because the only VM OF INFILTRATION IN TUMOR GROWTH 7
Time: t T u m o r r ad i u s : R ( t ) ( a ) Radius growth history. Radial coordinate: r C e ll nu m be r s Glioma cells: GImmune cells: MTotal: G+M ( b ) Cell numbers at T = . F igure A is to compute the nodal velocities U nj , the method we will propose later is inde-pendent of how A is computed, as long as the nodal U nj is computable; more details areprovided in the next section.2.3. A simple case study.
Whether (2.1) can be maintained by the solutions to (2.2) re-mains an open problem, since analytical approach to solve these equations remain di ffi cult.Nevertheless, one may justify that (2.1) should be respected by adding (2.2a) to (2.2c) toobtain: ∂ Θ ∂ t + r ∂∂ r (cid:104) r Θ V + r MU (cid:105) = f + g + h , where Θ = G + N + M ; and then incorporating (2.2i):(2.10) ∂ Θ ∂ t + r ∂∂ r (cid:104) r ( Θ − θ ) V (cid:105) = . Clearly, the equation (2.1), or equivalently Θ ≡ θ is a solution to the latest equation.In this section, we consider a simple case whose parameters and initial / boundary con-ditions are given as follows: • Most diminishing rates are set to zero except for λ , which models the self-productionof the glioma cells: λ = . , µ = δ = ρ = . . • We normalize the cell number by setting θ = .
0, and in the chemoattractant equa-tion: m = . , β = . , γ = . , ν = . , α = . . • The initial radius is R (0) =
1, and the initial cell numbers are: G ( r , = . , N ( r , = . , M ( r , = . , for all 0 ≤ r ≤ A ( r , = − r ≤ r ≤ e − ( r − r ≥ • The boundary condition for M is M bc = . T = . N η =
50 uniform intervals and fixed time step size ∆ t = . N is always zero (and so is our numerical solutions), hencewe clearly observe the violation of incompressibility in the numerical solutions at T = . X. ZENG, M. SALEH, AND J. TIAN
Time: t L - no r m o f G + M - F igure d θ .especially near the tumor center. To make this point clearer, the L -norm of G + M − d θ ( t n ) def == (cid:90) R ( t n )0 (cid:12)(cid:12)(cid:12) G ( r , t n ) + M ( r , t n ) − (cid:12)(cid:12)(cid:12) dr ≈ R n N η N η (cid:88) j = (cid:12)(cid:12)(cid:12)(cid:12) G nj − / + M nj − / − (cid:12)(cid:12)(cid:12)(cid:12) , where R n ≈ R ( t n ) is the numerical solution of the radius at t n . The history of d θ is providedin Figure 2.2, where we observe violation of the incompressibility constraint in increasingmagnitude as t grows. In the rest of the paper, we try to address this issue and investigateenhancement of existing finite volume methods to improve the numerical results.3. A M odel P roblem and T he T otality C onservation L aw To make the idea clear, we consider a simplified model instead of the original one. Firstof all, only two cell species are considered, namely the glioma cells G and the immunecells M . Second, noticing that the chemoattractant A is only used to compute the velocityfield U , in this simplified model we treat U as a given velocity field and denote it by u since it is prescribed; A is thusly ignored altogether. To this end, the governing equationsin spherical coordinate are given by: ∂ G ∂ t + r ∂∂ r (cid:104) r GV (cid:105) = f , (3.1a) ∂ M ∂ t + r ∂∂ r (cid:104) r M ( V + u ) (cid:105) = h , (3.1b) 1 r ∂∂ r (cid:104) r ( V + uM ) (cid:105) = f + h , V (0 , t ) = ≤ r ≤ R ( t ) and 0 ≤ t ≤ T ; the radius function R : [0 , T ] (cid:55)→ R + satisfies:(3.1d) R (cid:48) ( t ) = V ( R ( t ) , t ) . As before, the lower case letters in (3.1) represent prescribed functions:(3.1e) u = u ( r , t , G , M ) , f = f ( r , t , G , M ) , g = g ( r , t , G , M ) . We further require that u = r = t ∈ [0 , T ].Note that in (3.1c) there is no θ , comparing to the previous model; indeed, we havesupposed that θ = G ( r , + N ( r , = , ∀ ≤ r ≤ R (0) . Similar as before, we can define the total number Θ ( r , t ) = G ( r , t ) + H ( r , t ); then theincompressibility assumption requires Θ ≡
1. If this holds, we actually have a very con-venient way to estimate the growth of the tumor. In particular, let C ( t ) denote the total VM OF INFILTRATION IN TUMOR GROWTH 9 number of cells in the tumor; then on the one hand the assumption Θ ≡ C ( t ) = (cid:90) R ( t )0 π r Θ ( r , t ) dr = (cid:90) R ( t )0 π r dr = π R ( t ) , hence the rate of change in C ( t ) is:(3.4) C (cid:48) ( t ) = π R (cid:48) ( t ) R ( t ) . On the other hand, the only mechanism such that the new cells can enter the tumor isthrough the boundary condition for M at r = R ( t ):(3.5) C (cid:48) ( t ) = − π R ( t ) u ( R ( t ) , t ) ˇ M ( t ) , where u is the prescribed infiltration velocity and ˇ M ( t ) is the flow out of / into the tumor:(3.6) ˇ M ( t ) = (cid:40) M ( R ( t ) , t ) , u ( R ( t ) , t ) ≥ M bc ( t ) , u ( R ( t ) , t ) < , where as before M bc ( t ) is the prescribed ambient number of immune cells. Equating (3.4)and (3.5), we obtain an ODE for R ( t ):(3.7) R (cid:48) ( t ) = − u ( R ( t ) , t ) ˇ M ( t ) , which will help us design numerical tests for which the exact tumor growth curve can becalculated.3.1. The totality conservation law.
Adding (3.1a) and (3.1b) then equating the right handside with that of (3.1c), we obtain an analogy of (2.10): ∂ ( r Θ ) ∂ t + ∂∂ r (cid:104) r Θ V + r Mu (cid:105) = r ( f + g ) = ∂∂ r (cid:104) r ( V + uM ) (cid:105) , or equivalently:(3.8) ∂ ( r Θ ) ∂ t + ∂∂ r (cid:104) r ( Θ − V (cid:105) = , which admits the solution Θ ( r , t ) ≡ G and M by their sum Θ , an equivalent PDE system is obtained byreplacing either (3.1a) or (3.1b) by (3.8) without changing the solutions. Hence we expectthe incompressibility constraint G ( r , t ) + M ( r , t ) = Θ ( r , t ) ≡ totality conservation law or TCL in the context of current work and expect the numericalmethod satisfies a discrete version to be specified later.3.2.
The model, TCL, and GCL in normalized coordinate system.
Similar as before,after the coordinate transformation ( r , t ) (cid:55)→ ( η, τ ) = ( r / R ( t ) , t ), we obtain the model in thenormalized coordinates: ∂ ( η R G ) ∂τ + ∂∂η (cid:34)(cid:32) VR − η R (cid:48) R (cid:33) η R G (cid:35) = η R f − η R (cid:48) RG , (3.9a) ∂ ( η R M ) ∂τ + ∂∂η (cid:34)(cid:32) VR − η R (cid:48) R + uR (cid:33) η R M (cid:35) = η R g − η R (cid:48) RM , (3.9b) 1 η ∂∂η (cid:20) η (cid:18) VR + uR M (cid:19)(cid:21) = f + g , V (0 , t ) = the computational domain is ( η, τ ) ∈ [0 , × [0 , T ] for some positive T >
0; and R : [0 , T ] (cid:55)→ R + denotes the radious of the spherical domain, which satisfies:(3.9d) R (cid:48) ( τ ) = V (1 , τ ) . The lower case letters in (3.9) represent prescribed source terms.Correspondingly, the equation (3.8) is converted to: ∂ ( η R Θ ) ∂τ − η R (cid:48) R ∂∂η (cid:104) η R Θ (cid:105) + R ∂∂η (cid:104) η R ( Θ − V (cid:105) = , or equivalently:(3.10) ∂ ( η R Θ ) ∂τ + ∂∂η (cid:104) η R ( Θ − V (cid:105) − ∂∂η (cid:104) η R (cid:48) R Θ (cid:105) = − η R (cid:48) R Θ , which is the TCL in the normalized coordinates. In (3.8), both terms vanishes if we set Θ =
1; whereas in (3.10) setting
Θ = ∂ ( η R ) ∂τ + ∂∂η (cid:104) − η R (cid:48) R (cid:105) = − η R (cid:48) R . This equation only involve geometric quantities and it is rooted in using a mesh coordinate(normalized coordinate system) that is di ff erent from the physical one (the radial coordi-nates). Similar identities are studied in other contexts, especially the arbitrary Lagrangian-Eulerian (ALE) methods, see for example [10, 22], where it is called the geometric con-servation law or GCL . In this work, we follow this convention and call (3.11) the GCLfor the free-boundary problem in radial coordinates. At the continuous level, (3.11) holdsnaturally; but we will see in a moment that it may not hold at the discrete level. Existingliterature has demonstrated that violating GCL at the discrete level lead to unstable solu-tions; in this work, we thusly require the proposed method to satisfy a discrete version ofGCL, called the discrete geometric conservation law (DGCL), which will also be specifiedin the next section. 4. A n E nhanced F inite V olume M ethod Neither TCL nor GCL is automatically satisfied by classical finite volume discretiza-tions. For example to see why GCL could be violated, let us consider a numerical dis-cretization of (3.11), so that GCL is satisfied discretely, and look at what this discretizationmay look like. Using a mesh with N η uniform intervals and nodal velocities V nj where0 ≤ j ≤ N η , if the straightforward forward Euler time-integrator is used (see for example,Section 2.2), we have the following formula to update the solutions from t n to t n + :1 ∆ τ n (cid:104) η j − / ( R n + ) − η j − / ( R n ) (cid:105) + R (cid:48) n R n D j − / (cid:104) − η (cid:105) = − η j − / R (cid:48) n R n , (4.1a) R (cid:48) n = V nN η , (4.1b) R n + = R n + ∆ τ n V nN η , (4.1c)where (4.1a) collocates at the interval center η j − / and D j − / is the spatial discretizationfor ∂ η at η j − / as a result of the finite volume discretizations of (3.1). Rearranging (4.1a)there is: D j − / (cid:104) − η (cid:105) = − η j − / − η j − / ∆ τ n ( R n + ∆ τ n V nN η ) − ( R n ) V nN η R n = − η j − / − η j − / + ∆ τ n V nN η R n . VM OF INFILTRATION IN TUMOR GROWTH 11
This is a highly undesirable property, because it means that when we apply the chosennumerical discretization D j − / to a purely geometric quantity − η , the result needs todepend on the solutions of both V and R .An easy way to fix the issue is to make sure that the radius update satisfies:(4.2) R (cid:48) n = ( R n + ) − ( R n ) ∆ τ n R n , then (4.1a) reduces to:(4.3) D j − / ( − η ) = − η j − / , which is independent of V and R as desired.For example, if one wish to update the radius as R n + = R n + ∆ τ n R (cid:48) n , c.f., (4.1c), then itrequires R (cid:48) n to be computed as R (cid:48) n = V nN η (1 + (1 / ∆ τ n V nN η / R n ) rather than (4.1b). In thispaper, however, we propose to compute R (cid:48) n as:(4.4) R (cid:48) n = (cid:32) − ∆ η (cid:33) − V nN η , and then compute R n + according to (4.2). The motivation is to make sure that our numer-ical method is compatible with the no-flux biological condition at the moving boundary,see the discussion after the proof of (4.6).The preceding case study indicates that we must design the time-integrator carefully;furthermore, the spatial discretization D j − / needs to compute the derivative of third-degree polynomials exactly, as required by (4.3).The rest of this section focuses on constructing finite volume methods that satisfy boththe GCL and TCL in a discrete sense, which is yet to be made precise. To this end, wefollow the same notations as before and denote discrete cell numbers by G nj − / and M nj − / ,where 1 ≤ j ≤ N η , and they represent: G nj − / ≈ ∆ ηη j − / ( R n ) (cid:90) η j η j − η ( R n ) G ( η, τ n ) d η , (4.5) M nj − / ≈ ∆ ηη j − / ( R n ) (cid:90) η j η j − η ( R n ) M ( η, τ n ) d η , ≤ j ≤ N η ;(4.6)the discrete velocities are given by: V nj ≈ V ( η j , τ n ) , ≤ j ≤ N η , (4.7) u nj = u ( η j , τ n ) , ≤ j ≤ N η , (4.8)where no special approximation is needed for u since it can be evaluated explicitly, c.f. (3.1e).The remainder of this section is organized as follows. A general finite volume formu-lation is provided in Section 4.1 and our main result is given in Section 4.2, where bothDGCL and TGCL are defined and su ffi cient conditions for numerical methods to satisfythese conditions are provided. The subsequent sections then focus on various numericalfluxes that obey these conditions. A general finite volume formulation.
The explicit first-order time-accurate finitevolume formulation of (3.9a) and (3.9b) is obtained by integrating these equations overeach interval [ η j − , η j ] and then discretizing the time-derivative by forward-Euler method: η j − / [( R n + ) G n + j − / − ( R n ) G nj − / ] ∆ τ n + ∆ η (cid:104) F G , nj − F G , nj − (cid:105) (4.9) = η j − / ( R n ) f nj − / − η j − / R (cid:48) n R n G nj − / ,η j − / [( R n + ) M n + j − / − ( R n ) M nj − / ] ∆ τ n + ∆ η (cid:104) F M , nj − F M , nj − (cid:105) (4.10) = η j − / ( R n ) g nj − / − η j − / R (cid:48) n R n M nj − / , ≤ i ≤ N η . F G , nj and F M , nj are numerical fluxes for G and M at η j , respectively: F G , nj ≈ (cid:32) VR − η R (cid:48) R (cid:33) η R G (cid:12)(cid:12)(cid:12)(cid:12) η = η j , τ = τ n , (4.11) F M , nj ≈ (cid:32) VR − η R (cid:48) R + uR (cid:33) η R M (cid:12)(cid:12)(cid:12)(cid:12) η = η j , τ = τ n . (4.12)In Section 2.2, the velocities V / R − η R (cid:48) / R and V / R − η R (cid:48) / R + u / R are used to compute thetwo fluxes F Gj and F Mj , respectively. For our problems, however, it is advantageous toconsider each component of the velocity separately; namely, we segregate the numericalfluxes as: F G , nj = F G , nV , j + F G , nR (cid:48) , j , (4.13) F G , nV , j ≈ V η RG , F G , nR (cid:48) , j ≈ − η R (cid:48) RG ; F M , nj = F M , nV , j + F M , nR (cid:48) , j + F M , nu , j , (4.14) F M , nV , j ≈ V η RM , F M , nR (cid:48) , j ≈ − η R (cid:48) RM , F M , nu , j ≈ u η RM . The velocity equation is obtained similarly as before, but we keep the approximation to u η R M as unspecified:(4.15) η j R n V nj + F M , nu , j = j (cid:88) k = ∆ η (cid:16) η k − / ( R n ) f nk − / + η k − / ( R n ) g nk − / (cid:17) . Here F M , nu , j approximates u η R M at ( η j , τ n ); and the source terms on the right hand sideare computed the same way as those in (4.9) and (4.10). In Section 2.2, F M , nu , j is approxi-mated by averaging M nj − / and M nj + / ; as we will see soon, this is a good choice for ourproblem.4.2. Su ffi cient conditions for DTCL and DGCL. It is fair to assume that we use thesame flux function to compute the numerical fluxes associated with the same velocity,such as F G , nV , j and F M , nV , j ; to this end we suppose: F X , nV , j = F nj ( { X nj − / + k : − l ≤ k ≤ r } , P ) , (4.16) F X , nR (cid:48) , j = ˆ F nj ( { X nj − / + k : − l ≤ k ≤ r } , ˆ P ) , (4.17)where l ≥ r ≥ X repre-sents either species, and the parameter sets P and ˆ P are placeholders for high-resolutionfluxes that are described later. VM OF INFILTRATION IN TUMOR GROWTH 13
We distinguish the flux functions F nj and ˆ F nj because the former approximates the fluxesdue to a spatially varying velocity V whereas the latter can be interpreted as fluxes due to aspatially constant velocity R (cid:48) ; furthermore, we maintain the subscript j and the superscript n in these generic functions to indicate their dependence on the spatial coordinates η , do-main size R n , as well as R (cid:48) n , which are determined independently from the finite volumediscretizations.For our next purpose, we note that both flux functions are in the form F ( { X nj − / + k : − l ≤ k ≤ r } , · · · ), where the omitted quantities represent the parameters that are the samewhen the flux function is applied to compute fluxes for di ff erent species, such as G and M ,respectively. Definition 4.1.
The flux function F ( { X nj − / + k : − l ≤ k ≤ r } , · · · ) is called additive if forall X , Y and Z = X + Y : F ( { X nj − / + k : − l ≤ k ≤ r } , · · · ) + F ( { Y nj − / + k : − l ≤ k ≤ r } , · · · ) = F ( { Z nj − / + k : − l ≤ k ≤ r } , · · · ) . (4.18)where the omitted inputs are kept the same in all the three function evaluations.Furthermore, we define the V -consistency for the flux function F nj of (4.16) and cubic-preserving for the flux function ˆ F nj of (4.17) as follows. Definition 4.2.
The numerical flux function F nj of (4.16) is V-consistent if for all V nj :(4.19) F , nV , j = η j R n V nj , that is, setting X nj − / + k = , ∀ k in the right hand side of (4.16) yields η j R n V nj . Definition 4.3.
The numerical flux function ˆ F nj of (4.17) is cubic-preserving if(4.20) 1 ∆ η (cid:16) F , nR (cid:48) , j − F , nR (cid:48) , j − (cid:17) = − η j − / R (cid:48) n R n . Note that F X , nR (cid:48) , j can be treated as the flux for an advection equation with the spatiallyconstant velocity − R (cid:48) n R n and convected variable η X , this will be our basis to construct acubic-preserving flux function, see the further discussions in Section 4.5.The purpose of this section is to derive su ffi cient conditions such that our method satis-fies GCL and TCL discretely. To this end, we have the following definitions: Definition 4.4.
The method given by (4.9), (4.10) and (4.15) satisfies the DGCL providedthat: Suppose M mj − / + G mj − / = j and m = n , n +
1, then we can derive (4.15) from(4.9) and (4.10).
Definition 4.5.
The method given by (4.9), (4.10) and (4.15) satisfies the DTCL if theylead to a conservative discretization of (3.8).Now we state the main theorem that will eventually guide us in the construction of theenhanced numerical methods.
Theorem 4.6.
The numerical method given by (4.9), (4.10) and (4.15) satisfies both DGCLand DTCL if: (1) F nj is additive and V-consistent, (2) ˆ F nj is additive and cubic-preserving,(3) F M , nu , j = F M , nu , j , and (4) R (cid:48) n equals the right hand side of (4.2). Proof.
Adding (4.9) and (4.10) then incorporating (4.15), we have: η j − / ∆ τ n (cid:104) ( R n + ) ( G n + j − / + M n + j − / ) − ( R n ) ( G nj − / + M nj − / ) (cid:105) + ∆ η (cid:16) F M , nu , j − F M , nu , j − (cid:17) ∆ η (cid:16) F G , nV , j + F M , nV , j − F G , nV , j − − F M , nV , j − (cid:17) + ∆ η (cid:16) F G , nR (cid:48) , j + F M , nR (cid:48) , j − F G , nR (cid:48) , j − − F M , nR (cid:48) , j − (cid:17) = η j − / ( R n ) f nj − / + η j − / ( R n ) g nj − / − η j − / R (cid:48) n R n ( G nj − / + M nj − / )(4.21) = R n ( η j V nj − η j − V nj − ) ∆ η + ∆ η (cid:16) F M , nu , j − F M , nu , j − (cid:17) − η j − / R (cid:48) n R n ( G nj − / + M nj − / ) . Define Θ nj − / = G nj − / + M nj − / and Θ n + j − / = G n + j − / + M n + j − / as before; following theadditivity of the fluxes F nj and ˆ F nj we obtain: F G , nV , j + F M , nV , j = F Θ , nV , j and F G , nR (cid:48) , j + F M , nR (cid:48) , j = F Θ , nR (cid:48) , j . Invoking in addition the assumption that F M , nu , j = F M , nu , j , we obtain from (4.21): η j − / ∆ τ n (cid:104) ( R n + ) Θ n + j − / − ( R n ) Θ nj − / (cid:105) + ∆ η (cid:16) F Θ , nV , j − F Θ , nV , j − (cid:17) + ∆ η (cid:16) F Θ , nR (cid:48) , j − F Θ , nR (cid:48) , j − (cid:17) (4.22) = R n ( η j V nj − η j − V nj − ) ∆ η − η j − / R (cid:48) n R n Θ nj − / . Clearly (4.22) represents a conservative finite volume discretization of the continuous to-tality conservation law (3.8) using the same flux functions F nj and ˆ F nj ; hence the methodsatisfies DTCL.Now we move on to show DGCL and to this end assume Θ nj ≡ Θ n + j ≡
1, then(4.22) reduce to: η j − / ∆ τ n (cid:104) ( R n + ) − ( R n ) (cid:105) + ∆ η (cid:16) F , nV , j − F , nV , j − (cid:17) + ∆ η (cid:16) F , nR (cid:48) , j − F , nR (cid:48) , j − (cid:17) (4.23) = R n ( η j V nj − η j − V nj − ) ∆ η − η j − / R (cid:48) n R n . Since R (cid:48) n = (( R n + ) − ( R n ) ) / (2 ∆ τ n R n ), (4.23) is equivalent to:(4.24) 1 ∆ η (cid:16) F , nV , j − F , nV , j − (cid:17) + ∆ η (cid:16) F , nR (cid:48) , j − F , nR (cid:48) , j − (cid:17) = R n ( η j V nj − η j − V nj − ) ∆ η − η j − / R (cid:48) n R n . This equality is trivial to prove following the V -consistency of F nj and the cubic-preservingof ˆ F nj . Hence we conclude that given all the assumptions as stated, and that M mj + G mj =Θ mj = , ∀ j and m = n , n +
1, (4.9) and (4.10) gives rise to (4.15). Thus the method satisfiesDGCL. (cid:3)
In the theorem and its proof, we only considered the radius update condition (4.2). Onthe one hand, the theorem only requires R n , R (cid:48) n , and R n + to be related by (4.2); and it doesnot pose any restriction on how R (cid:48) n is to be computed. On the other hand, biologicallypeople do not expect any G to flow across the moving boundary, which translates to:(4.25) F , nV , N η + F , nR (cid:48) , N η = , VM OF INFILTRATION IN TUMOR GROWTH 15 and no geometrical flux at η = F , nR (cid:48) , = . However, the V -consistency condition requires that: F , nV , N η = η N η R n V nN η = R n V nN η , and incorporating (4.26), the cubic-preserving condition requires: F , nR (cid:48) , N η = F , nR (cid:48) , + N η (cid:88) j = ( F , nR (cid:48) , j − F , nR (cid:48) , j − ) = − ∆ η N η (cid:88) j = η j − / R (cid:48) n R n = (cid:32) − ∆ η (cid:33) R (cid:48) n R n . Hence the no-flux condition (4.25) indicates V nN η = (cid:16) − ∆ η (cid:17) R (cid:48) n R n , or equivalently (4.4)as proposed before.4.3. A review of the conventional flux functions.
We briefly review the conventionalfinite volume method in the context of (3.9); particularly we consider the spherically sym-metric conservation law for a generic species X in spherical coordinates and radial advec-tive velocity W :(4.27) ∂ ( η R X ) ∂τ + ∂∂η (cid:104) W ( η R X ) (cid:105) = , where we omitted any source terms on the right hand side since their approximation isgenerally independent of the finite volume discretizations.The conservative variable of (4.27) is ˜ X def == η R X rather than X , particularly the variablefor the interval [ η j − , η j ] is ˜ X j − / = η j − / R X j − / . Hence (4.27) is simply the conserva-tive advection equation for ˜ X by the velocity W :(4.28) ∂ ˜ X ∂τ + ∂∂η (cid:16) W ˜ X (cid:17) = , whose finite volume discretization (at the semi-discretized level) reads: d ˜ X j − / d τ + ∆ η (cid:16) F j − F j − (cid:17) = , where F j ≈ W ˜ X (cid:12)(cid:12)(cid:12) η = η j and F j − ≈ W ˜ W (cid:12)(cid:12)(cid:12) η = η j − .If the conventional first-order upwind flux is used (see Section 2.2), there is:(4.29) F j = F upw ( W j ; ˜ X j − / , ˜ X j + / ) , where W j is the nodal velocity at η j and F upw is given by (2.8).Extension to higher accuracy is achieved by the limited polynomial reconstruction. Oneof the most widely used second-order extension is given by the high-resolution MUSCLmethod [25]:(4.30) F j = F muscl ( W j ; ˜ X j − / , ˜ X j − / , ˜ X j + / , ˜ X j + / , φ j − / , φ j + / ) , where φ j − / and φ j + / are slope limiters and the MUSCL flux function is: F muscl ( W j ; Z j − / , Z j − / , Z j + / , Z j + / , φ j − / , φ j + / )(4.31) def == F upw (cid:32) W j ; Z j − / + φ j − / ∆ Z j , Z j + / − φ j + / ∆ Z j + (cid:33) , where ∆ Z k = Z k + / − Z k − / , k = j , j + Z is a generic variable that equals ˜ X in thecase of (4.30). The slope limiter φ j − / ∈ [0 ,
1] usually depends on the solutions, but only weakly in the following sense. Slope limiters are introduced to reduce the magnitude ofthe slope such that the reconstruction will not create any new local extremum – a propertycalled monotone preserving. Hence if setting φ j − / = c satisfies the monotone preservingproperty for some particular value c , so is all slope limiters φ j − / ∈ [0 , c ]. For this reasonthe slope limiters are introduced as free (or more precisely semi-free) parameters.The bounds for slope limiters are nonlinear functions of the discrete solutions, for ex-ample, the minmod limiter computes:(4.32) φ j − / = ϕ minmod ( ∆ Z j − , ∆ Z j ) = , ∆ Z j − ∆ Z j ≤ , min (cid:18) ∆ Z j − ∆ Z j , (cid:19) , ∆ Z j − ∆ Z j > . Other widely used limiter functions can be found in [20, 19, 26].Meanwhile, we show that these conventional flux functions are neither V -consistent norcubic-preserving. The latter is easy to verify; indeed, the upwind flux is only first-orderaccurate and the MUSCL flux is at most second-order; whereas cubic-preserving requiresa third-order flux for advection equations.Let us focus on the V -consistency and consider, for example, the upwind flux F upw .Then V -consistency requires that F upw ( V j / R ; η j − / R , η j + / R ) = η j RV j ; however, thisequality does not hold either when V j ≥
0, in which case according to (2.8): F upw (cid:32) V j R ; η j − / R , η j + / R (cid:33) = η j − / RV j (cid:44) η j RV j ;or when V j >
0, in which case: F upw (cid:32) V j R ; η j − / R , η j + / R (cid:33) = η j + / RV j (cid:44) η j RV j . In the next sub-sections, we focus on designing numerical methods such that they leadto a method that satisfies both DGCL and DTCL, following the results of Section 4.2.4.4.
Modified fluxes: Part I.
In this section, we construct V -consistent fluxes F nj bymodifying the conventional upwind or MUSCL fluxes; in the latter case a synchronizedlimiter is introduced to ensure the additivity property as required by Theorem 4.6. Thefluxes F X , nV , j and F X , nu , j will subsequently be constructed accordingly.To construct a V -consistent flux F nj , instead of applying the conventional flux functionsto the conservative variables ˜ X , we consider the primitive ones X . Particularly, a first-orderupwind method for (4.16) can be constructed by setting l = r =
1, and P = ∅ :(4.33) F nj ( { X nj − / , X nj + / } ) = η j ( R n ) F upw V nj R n ; X nj − / , X nj + / . Because F upw ( V nj / R n ; 1 , ≡ V nj / R n , the flux (4.33) is V -consistent.Similarly, extension to higher-order accuracy can make use of the MUSCL flux (4.31): F nj ( { X nj − / , X nj − / , X nj + / , X nj + / } , { φ X , nj − / , φ X , nj + / } )(4.34) = η j ( R n ) F muscl V nj R n ; X nj − / , X nj − / , X nj + / , X nj + / , φ X , nj − / , φ X , nj + / ,φ X , nk − / = ϕ minmod ( ∆ X nk − , ∆ X nk ) , k = j , j + . (4.35)Here P = { φ X , nj − / , φ X , nj + / } and the minmod limiter can be replaced by any other limiterof choice. It is not di ffi cult to verify that if X nk − / ≡
1, the MUSCL flux F muscl gives VM OF INFILTRATION IN TUMOR GROWTH 17 rise to V nj / R n regardless of the values of the limiters; hence the flux function (4.34) is V -consistent, no matter what limiter we will choose.Next the additivity of these fluxes is considered, which is essentially requiring that thefluxes are linear in the inputs X j − / . Hence the upwind fluxes are by nature additive; forexample let us consider F nj given by (4.33) and suppose V nj ≥
0, then: F G , nV , j + F M , nV , j = F nj ( { G nj − / , G nj + / } ) + F nj ( { M nj − / , M nj + / } ) = η j ( R n ) F upw V nj R n , G nj − / , G nj + / + η j ( R n ) F upw V nj R n , M nj − / , M nj + / = η j ( R n ) · V nj R n G nj − / + η j ( R n ) · V nj R n M nj − / = η j ( R n ) · V nj R n Θ nj − / = η j ( R n ) F upw V nj R n , Θ nj − / , Θ nj + / = F Θ , nV , j . The argument for the case V nj < F nj is additive.Extension to the MUSCL-based fluxes (4.34) is not straightforward, as the limiter func-tion is generally nonlinear. Following the discussion below Equation (4.31), we can cir-cumvent this di ffi culty by synchronizing the limiters for G and M , that is: F G , nV , j = η j ( R n ) F muscl V nj R n , G nj − / , G nj − / , G nj + / , G nj + / , φ nj − / , φ nj + / (4.36) F M , nV , j = η j ( R n ) F muscl V nj R n , M nj − / , M nj − / , M nj + / , M nj + / , φ nj − / , φ nj + / (4.37) where φ nk − / = min (cid:16) φ G , nk − / , φ M , nk − / (cid:17) , k = j , j + , (4.38)here φ G , nk − / and φ M , nk − / are obtained by applying (4.35) to X = G and X = M , respectively.Note that the same limiters are used to compute the two fluxes. To show the additivity, weassume again without loss of generality that V nj ≥
0, then: F G , nV , j = η j ( R n ) · V nj R n (cid:32) G nj − / + φ nj − / ( G nj + / − G nj − / ) (cid:33) , F M , nV , j = η j ( R n ) · V nj R n (cid:32) M nj − / + φ nj − / ( M nj + / − M nj − / ) (cid:33) ⇒ F G , nV , j + F M , nV , j = η j ( R n ) · V nj R n (cid:32) Θ nj − / + φ nj − / ( Θ nj + / − Θ nj − / ) (cid:33) = η j ( R n ) F muscl V nj R n ; Θ nj − / , Θ nj − / , Θ nj + / , Θ nj + / , φ nj − / , φ nj + / . The latest flux is in general not a monotone flux for Θ , since the selected limiter may betoo large. Nevertheless, this is not an issue since Θ is not our numerical solution; and aslong as G and M are computed using monotone fluxes, we won’t run into stability issues.Nevertheless, if one wishes to ensure monotone flux for Θ as well, all that needs to be doneis to compute the limiter φ Θ , nk − / by applying (4.35) to X = G + N , and include this φ Θ , nk − / inthe minimum of the right hand side of (4.38). Finally, we compute the u -fluxes for M by:(4.39) F M , nu , j = η j ( R n ) F upw u nj R n ; M nj − / , M nj + / for first-order accuracy or:(4.40) F M , nu , j = η j ( R n ) F muscl u nj R n ; M nj − / , M nj − / , M nj + / , M nj + / , φ nj − / , φ nj + / for second-order accuracy, where the limiters are the same ones computed by (4.38).4.5. Modified fluxes: Part II.
To construct a cubic-preserving flux F X , nR (cid:48) , j , however, wecannot follow the same strategy as in the previous section. Indeed, if this flux is defined as: F X , nR (cid:48) , j = η j ( R n ) F upw (cid:32) − η j R (cid:48) n R n , X nj − / , X nj + / (cid:33) , supposing R (cid:48) n ≥ X nk − / ≡ F , nR (cid:48) , j − = − η j − R (cid:48) n R n , F , nR (cid:48) , j = − η j R (cid:48) n R n ⇒ ∆ η (cid:16) F , nR (cid:48) , j − F , nR (cid:48) , j − (cid:17) = − (cid:32) η j − / + ∆ η (cid:33) R (cid:48) n R n , which is di ff erent from − η j − / R (cid:48) n R n , as required by the cubic-preserving property.To proceed, we recognize that a higher-order and nonlinearly stable flux can be obtainedby a polynomial reconstruction of the solutions on each interval such that the total variationdoes not increase, and apply the upwind flux to the two reconstructed values on both sidesof the interval face. In the MUSCL scheme, the reconstruction is achieved by limiting theslope of a linear function that preserves the interval average; in this section, we adopt theaverage-preserving and monotone cubic reconstruction of the Piecewise Parabolic Method(PPM) [7], but construct the flux di ff erently.Let X be a generic variable as before, the PPM reconstruction in normalized coordinates ξ = ( η − η j − ) / ∆ η on the interval [ η j − , η j ] reads: X ( ξ ) = X j − / , − + ξ (cid:16) X j − / , + − X j − / , − + X , j − / (1 − ξ ) (cid:17) , (4.41) X , j − / == X j − / − X j − / , − + X j − / , + ) . (4.42)Here X j − / , − and X j − / , + are the two end values that are defined as: X j − / , − = X j − / + φ Xj − / , − ( X j − − X j − / ) , (4.43) X j − / , + = X j − / + φ Xj − / , + ( X j + − X j − / ) , (4.44) where X k def ==
712 ( X k − / + X k + / ) −
112 ( X k − / + X k + / ) , k = j − , j . (4.45)The value X k , k = j − , j are third-order reconstructions of the face values when the datais smooth; and the two limiters φ Xj − / , ± are decides as follows:(a) If ( X j − X j − / )( X j − − X j − / ) ≥
0, we have a local extrema and set:(4.46) φ Xj − / , − = φ Xj − / , + = . VM OF INFILTRATION IN TUMOR GROWTH 19 (b) If (a) is not true, and if (cid:12)(cid:12)(cid:12) X j − X j − / (cid:12)(cid:12)(cid:12) > (cid:12)(cid:12)(cid:12) X j − − X j − / (cid:12)(cid:12)(cid:12) or (cid:12)(cid:12)(cid:12) X j − − X j − / (cid:12)(cid:12)(cid:12) > (cid:12)(cid:12)(cid:12) X j − X j − / (cid:12)(cid:12)(cid:12) ,the corresponding reconstructed profile is not monotone on the interval [ η j − , η j ] andwe compute:(4.47) φ Xj − / , + = − X j − − X j − / ) X j − X j − / in the former case, and(4.48) φ Xj − / , − = − X j − X j − / ) X j − − X j − / in the latter case.(c) Otherwise, set the remaining limiter, which is φ Xj , − , or φ Xj , + , or both, to one.Finally, denoting X nj − / = η j − / X nj − / we compute the flux ˆ F nj of (4.17) as:ˆ F nj (cid:16) { X nj − / , X nj − / , X nj − / , X nj + / , X nj + / , X nj + / } , { φ X , nj − / , + , φ X , nj + / , − } (cid:17) (4.49) def == F upw (cid:16) − R (cid:48) n R n ; X nj − / , + , X nj + / , − (cid:17) . Here X nk − / , ± are computed according to (4.43) and (4.44) with data X nk − / ; and the lim-iters φ X , nj − / , ± are computed as φ Xj − / , ± according to (a-c) given previously. Theorem 4.7.
The flux F X , nR (cid:48) , j = ˆ F nj , which is given by (4.49), satisfies (4.20), i.e.,: ∆ η (cid:16) F , nR (cid:48) , j − F , nR (cid:48) , j − (cid:17) = − η j − / R (cid:48) n R n for j ≥ .Proof. We suppress the superscript n for simplicity and let X k − / ≡ , ∀ k , then the recon-structed values X k , k ≥ X k = (cid:16) η k − / + η k + / (cid:17) − (cid:16) η k − / + η k + / (cid:17) = η k − η k ∆ η . Now we compute the limiters φ Xk − / , ± for k ≥
3. First of all: X k − X k − / = η k − η k ∆ η − η k − / = η k − / ∆ η + η k − / ∆ η > . X k − − X k − / = η k − − η k − ∆ η − η k − / = − η k − / ∆ η + η k − / ∆ η < , thus the condition in (a) does not hold. Continuing to check the conditions in (b), we have:2 (cid:12)(cid:12)(cid:12) X k − X k − / (cid:12)(cid:12)(cid:12) − (cid:12)(cid:12)(cid:12) X k − − X k − / (cid:12)(cid:12)(cid:12) = η k − / ∆ η + η k − / ∆ η > , (cid:12)(cid:12)(cid:12) X k − − X k − / (cid:12)(cid:12)(cid:12) − (cid:12)(cid:12)(cid:12) X k − X k − / (cid:12)(cid:12)(cid:12) = η k − / ∆ η − η k − / ∆ η > , hence neither condition in (b) is true. Thus we conclude that φ Xk − / , ± = k ≥ X k − / , − = X k − and X k − / , + = X k .Finally we can calculate the flux F , nR (cid:48) , j for all j ≥
3. Because X j + / , − = X j = X j − / , + ,regardless of the sign of R (cid:48) n , there is:(4.51) F , nR (cid:48) , j = − R (cid:48) n R n X j = − R (cid:48) n R n (cid:32) η j − η j ∆ η (cid:33) , thus for j , j − ≥ ∆ η (cid:16) F , nR (cid:48) , j − F , nR (cid:48) , j − (cid:17) = − R (cid:48) n R n ∆ η (cid:34)(cid:32) η j − η j ∆ η (cid:33) − (cid:32) η j − − η j − ∆ η (cid:33)(cid:35) = − R (cid:48) n R n ∆ η (cid:34)(cid:32) ( η j − / + ∆ η ) − ( η j − / − ∆ η ) (cid:33) −
14 ( η j − η j − ) ∆ η (cid:35) = − R (cid:48) n R n ∆ η (cid:32) η j − / ∆ η + ∆ η − ∆ η (cid:33) = − η j − / R (cid:48) n R n . This concludes that the constructed flux is cubic-preserving. (cid:3)
Theorem 4.7 addresses the cubic-preserving for intervals far away from the bound-aries; next we consider boundary intervals and focus on those near the origin first. Morespecifically, we need to consider the cubic-preserving on the first three intervals [0 , ∆ η ],[ ∆ η, ∆ η ], and [2 ∆ η, ∆ η ]. Following a similar procedure in the preceding proof, cubic-preserving on these three intervals amounts to: F , nR (cid:48) , = − R (cid:48) n R n (cid:32) η − η ∆ η (cid:33) , (4.52a) F , nR (cid:48) , = − R (cid:48) n R n (cid:32) η − η ∆ η (cid:33) , (4.52b) F , nR (cid:48) , = . (4.52c)Note that (4.52c) is enforced automatically as the boundary condition at the origin; hencewe focus on the first two, which are guaranteed if: X = ∆ η , X = ∆ η , and that φ X , n / , ± = φ X , n / , ± = X ’s are 1. The required X is precisely given by theformula in (4.43) or (4.43); thus we just need to look at X and the limiters. To this end,utilizing the knowledge that when X is defined as η X , we expect X = η = X by the modified formula:(4.53) X =
712 ( X / + X / ) −
112 ( X / − X / ) . According to this definition, φ X , n / , ± = X ≡ φ X , n / , − = φ X , n / , + = / , ∆ η ] utilizing the knowledge that X ∼ η near η = X ( ξ ) = ξ ( X / , + + X , / (1 − ξ )) , X , / == X / − X / , + , X / , + = X / + φ X / , + ( X − X / ) . As before X , / is defined such that the mean of X ( ξ ) is X / for any right end value X / , + .Because X / , − ≡ X ≡
0, there is no way to design the limiter such that X ( ξ ) is monotoneif X / (cid:44) X is too close to X / . To this end, we design φ X / , + as follows instead ofthe generic construction for other intervals:(a’) If X / X ≤ (cid:12)(cid:12)(cid:12) X (cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12) X / (cid:12)(cid:12)(cid:12) , we set φ X / , + = VM OF INFILTRATION IN TUMOR GROWTH 21 (b’) If (a’) is not true, and if (cid:12)(cid:12)(cid:12) X (cid:12)(cid:12)(cid:12) > (cid:12)(cid:12)(cid:12) X / (cid:12)(cid:12)(cid:12) , we define(4.54) φ X / , + = X / X / − . (c’) Otherwise, φ X / , + = X / = η / = ∆ η and X = ∆ η , the limiter isprecisely φ X , n / , + = F G , nR (cid:48) , j and F M , nR (cid:48) , j by constructing the limiter φ nj , ± for both species as follows. For j ≥ M j − M j − / )( M j − − M j − / ) ≥ G j − G j − / )( G j − − G j − / ) ≥
0, we set:(4.55) φ nj − / , − = φ nj − / , + = . (B) Otherwise, we compute: α = min (cid:12)(cid:12)(cid:12) M j − − M j − / (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) M j − M j − / (cid:12)(cid:12)(cid:12) , (cid:12)(cid:12)(cid:12) G j − − G j − / (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) G j − G j − / (cid:12)(cid:12)(cid:12) , and α = max (cid:12)(cid:12)(cid:12) M j − − M j − / (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) M j − M j − / (cid:12)(cid:12)(cid:12) , (cid:12)(cid:12)(cid:12) G j − − G j − / (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) G j − G j − / (cid:12)(cid:12)(cid:12) . (B1) If α > α , use (4.55).(B2) Otherwise if α <
1, set:(4.56) φ nj − / , − = , φ nj − / , + = α ;and if α >
1, set:(4.57) φ nj − / , − = α − , φ nj − / , + = . Similarly for the first interval, we have:(A’) If X / X ≤ (cid:12)(cid:12)(cid:12) X (cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12) X / (cid:12)(cid:12)(cid:12) is true for either X = G or X = M , we set φ n / , + = φ n / , + = min , G / G / − , M / M / − . Modified fluxes: Part III.
In the last part of the modified fluxes, we consider theintervals near the right boundary η =
1. Particularly, these are the intervals whose fluxcalculation requires solutions beyond the computational domain.At η N η , we notice that the velocity V / R − η R (cid:48) / R ≡ F X , nV , N η + F X , nR (cid:48) , N η = , where X = G or X = M . In fact, practically we set both fluxes F X , nV , N η and F X , nR (cid:48) , N η to zero forconvenience.The remaining flux F M , nu , N η is computed as:(4.60) F M , nu , N η = η N η ( R n ) F upw u nN η R n ; M nN η − / , M bc ( t n ) , no matter which flux function is used for interior nodes.At η N η − , calculating F X , nV , N η − using the upwind flux (4.33) does not require any intervalsbeyond η = X N η + / . To this end, we avoid the linear reconstruction on the lastinterval and modify the numerical flux as:(4.61) F X , nV , N η − = η N η − ( R n ) F upw V nN η − R n ; X nN η − / + φ nN η − / ∆ X nN η − , X nN η − / , where X = G or X = M and φ nN η − / is computed according to (4.38). Similarly, the flux F M , nu , N η − is given by:(4.62) F M , nu , N η − = η N η − ( R n ) F upw u nN η − R n ; M nN η − / + φ nN η − / ∆ M nN η − , M nN η − / , where φ nN η − / is the same limiter used in (4.61).In order to compute F X , nR (cid:48) , N η − using the modified PPM method in Section 4.5, we needto use a biased stencil to interpolate X to obtain:(4.63) X N η − =
112 (3 X N η − / + X N η − / − X N η − / + X N η − / ) . Using this definition, when X ≡ X N η − = η N η − − η N η − ∆ η /
4, c.f. (4.50). Fol-lowing the same procedure in the proof of Theorem 4.7, we have φ XN η − / , ± = η N η − , η N η − ].Finally, to ensure cubic-preserving on the next interval [ η N η − , η N η − ], all that needs tobe done is to design X N η properly such that the limiter φ XN η − / , − takes the value 1 when X ≡
1. This can be achieved by the extrapolating formula:(4.64) X N η =
112 (25 X N η − / − X N η − / + X N η − / − X N η − / ) . Note that X N η is only used for computing the limiter φ XN η − / , − .4.7. Higher-order accuracy in time.
The preceding sections fully specify the discretiza-tion with first-order and second-order accuracy in space, and first-order accuracy in time.Here the temporal integration is achieved by the forward Euler (FE) method; hence ex-tension to higher-order accuracy in time can be easily achieved by using Total VariationDiminishing (TVD) Runge-Kutta methods [24, 14]. In particular, we consider the second-order TVD Runge-Kutta, denoted by TVD-RK2 in the rest of the paper, to match thespatial order of accuracy when the MUSCL fluxes are used. For this purpose, we denotethe solution as S = { G , M , V , R } and let the method proposed before with FE integrator besummarized as:(4.65) S n + = M ∆ τ n ( S n ) , here the subscript ∆ τ n denotes the time-step size; then the method using TVD RK2 reads: S (1) = M ∆ τ n ( S n ) , S (2) = M ∆ τ n ( S (1) ) , (4.66) S n + = (cid:16) S n + S (2) (cid:17) . VM OF INFILTRATION IN TUMOR GROWTH 23
Here the average is defined in the natural way, for example, G n + j = ( G nj + G (2) j ) / R n + = ( R n + R (2) ) /
2, etc.Before concluding this section, we have three remarks.
Remark 1 . The computation of the time step size ∆ τ n is according to the classicalCourant condition for linear stability. However, the original formula needs to be adjustedsince segregated advection velocities are used in the enhanced methods. In the case ofthe upwind flux combining with first-order explicit time-integrator, the method remainsconditionally stable and the analysis as well as the formula for the corresponding Courantcondition are provided in Appendix A. Remark 2 . The methodology extends naturally to the original tumor growth prob-lem. In particular, the fluxes for G , N , M are segregated similarly; and in extension tothe MUSCL flux or PPM flux, the limiter synchronization needs to take into account allspecies. Remark 3 . When the tumor growth model (2.2) is considered, one typically requires animplicit time-integrator for updating the chemoattractant concentration A to avoid tiny timestep sizes. It is then very natural to ask how the enhancement can be applied with implicitsolvers. We briefly address this issue in Appendix B in the case of the backward-Eulermethod and the class of Diagonally Implicit Runge-Kutta (DIRK) methods, and providebrief numerical results for comparison.5. N umerical A ssessment We assess the performance of the numerical methods of Section 4 using various bench-mark tests. First, we consider the model problem (3.9) and verify the DTCL and DGCLproperties of the proposed method.5.1.
The model problem.
To assess the numerical performance, we consider a series oftests that are characterized by spatially constant solutions, “prescribed” growth, and non-monotonic radius change, respectively. For all the tests, we set the initial radius as R (0) = .
0. The purpose of these tests is to assess the ability of the enhanced methods to satisfy thetotality conservation law and geometric conservation law discretely. Hence for each testbelow, we compare the numerical results that are obtained by using four di ff erent methods: • The conventional finite volume method with upwind fluxes as described in Sec-tion 2.2. • The conventional finite volume method with MUSCL fluxes with the TVD-RK2time-integrator as described in Section 4.7. • The enhanced finite volume method with upwind V and u fluxes and cubic-preserving R (cid:48) fluxes, as required by Theorem 4.6. • The enhanced finite volume method with MUSCL V and u fluxes and cubic-preserving R (cid:48) fluxes, as required by Theorem 4.6, and TVD-RK2 time-integratoras described in Section 4.7.The first two methods are denoted by “Conv. Upwind” and “Conv. MUSCL”, respectively;and the two enhanced ones are denoted by “Enhc. Upwind” and “Enhc. MUSCL”, respec-tively, in the subsequent tests. For all the methods, the fixed Courant number α cfl = . Test 1: Spatially constant solutions – single species.
In the first two tests, the ve-locity field u is manufactured such that (3.1) allows solutions that are independent of the Time: t L - no r m o f G + M - Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL F igure d θ on a 50-interval grid.spatial coordinate. Indeed, assuming f ≡ h ≡ G ( r , t ) = G ( t ), (3.1a) indicates: G (cid:48) ( t ) + r ∂∂ r ( r V ( r , t )) G ( t ) = . Thus V has to vary linearly in r and with abuse of notation we write V ( r , t ) = rV ( t ), where V ( t ) satisfies: G (cid:48) ( t ) + V ( t ) G ( t ) = ⇒ G ( t ) = e − (cid:82) t V ( s ) ds G (0) . Similarly, supposing further M ( r , t ) = M ( t ), (3.1b) indicates u ( r , t ) also varies linearly in r and u ( r , t ) = ru ( t ); hence we have: M (cid:48) ( t ) + V ( t ) + u ( t )) M ( t ) = ⇒ M ( t ) = e − (cid:82) t ( V ( s ) + u ( s )) ds M (0) . Lastly, (3.1c) is satisfied if and only if V ( t ) + u ( t ) M ( t ) ≡
0. To this end, we let V ( t ) = V > G ( t ) = e − V t G (0) , M ( t ) = − e − V t G (0) , u ( t ) = − V M ( t ) = − V − e − V t G (0) . It is easy to check that (5.1) indeed solves the model problem, providing the boundary data:(5.2) M bc ( t ) = − e − V t G (0) . In this case, the radius growth is exponential:(5.3) R (cid:48) = RV ⇒ R ( t ) = e V t . In the first test, we set V = . G ( r , = . , M ( r , = . . It is easy to tell that when the initial data G ( r ,
0) is zero, so is G ( r , t ) for all t >
0. This isindeed satisfied by all our numerical solutions, whether using the conventional methods orthe enhanced ones; hence we will not plot G in the next results. Solving the problem on agrid of 50 uniform intervals until T = .
0, the histories of the incompressibility constraintviolation index (2.11) are plotted in Figure 5.1. We clearly observe that the incompressibil-ity constraint is satisfied by both enhanced methods, whereas both conventional methodslead to increasing violation of this constraint as t grows.In Figure 5.2, we plot the radius histories and the profile of M ( r , T ) in the left panel andthe right panel, respectively. Here in Figure 5.2a, the reference radius growth curve (5.3)is plotted against the numerical ones; and we can see that all numerical solutions are closeto the reference one, with the MUSCL fluxes provide slightly more accurate results thanthe upwind ones. Figure 5.2b show that conventional methods fail to preserve constantsolutions for a single species, indicating the violation of the geometrical conservation law;whereas both enhanced methods satisfy DGCL.The same tests are performed on finer grids, with 100, 200, and 400 cells, respectively;and we have very similar plots as before. In Figure 5.3, the final radius is plotted for each VM OF INFILTRATION IN TUMOR GROWTH 25
Time: t T u m o r r ad i u s : R ( t ) ExactConv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( a ) Radius growth history. Radial coordinate: r M Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( b ) Cell numbers for M at T = . F igure
50 100 200 400
Number of intervals T e r m i na l r ad i u s ExactConv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL F igure able
1. Numerical errors in radius of test 1 at T = . N η Conv. Upwind Conv. MUSCL Enhc. Upwind Enhc. MUSCLError Rate Error Rate Error Rate Error Rate50 -6.97e-3 -1.86e-2 -2.12e-2 -8.37e-3100 -4.01e-3 0.80 -8.46e-3 1.14 -1.07e-2 0.98 -4.27e-3 0.97200 -2.21e-3 0.86 -3.67e-3 1.20 -5.40e-3 0.99 -2.15e-3 0.99400 -1.18e-3 0.91 -1.55e-3 1.24 -2.71e-3 1.00 -1.08e-3 0.99T able L -errors in M ( η, T ) of test 1 at T = . N η Conv. Upwind Conv. MUSCL Enhc. Upwind Enhc. MUSCLError Rate Error Rate Error Rate Error Rate50 6.03e-2 1.38e-2 3.89e-16 8.62e-16100 3.46e-2 0.80 6.67e-3 1.05 8.18e-16 N / A 5.94e-16 N / A200 1.95e-2 0.83 3.24e-3 1.04 3.71e-15 N / A 7.43e-16 N / A200 1.09e-2 0.85 1.58e-3 1.04 9.80e-16 N / A 1.34e-16 N / Amethod on the sequence of four grids; and they’re compared to the exact value. In addition,quantitative comparison is provided in Table 1 and Table 2, which summarizes the errorsin the final radius and the numerical solutions in M , respectively. In order to evaluate theerrors in M , we consider the L -error in the normalized coordinate that is defined as: N η (cid:88) i = ∆ η (cid:12)(cid:12)(cid:12)(cid:12) M N τ k − / − M ∗ ( η k − / , T ) (cid:12)(cid:12)(cid:12)(cid:12) , where N τ is the time step at T and M ∗ is the exact solution given by (5.1); for this particularproblem, we have M ∗ ≡
1. From Table 2, we see that the numerical error in M by theenhanced methods is at the scale of the machine accuracy, which indicates that they satisfy Radial coordinate: r G Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( a ) Cell numbers for G at T = . Radial coordinate: r M Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( b ) Cell numbers for M at T = . Time: t T u m o r r ad i u s : R ( t ) ExactConv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( c ) Radius growth history. Time: t L - no r m o f G + M - Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( d ) Histories of d θ . F igure
50 100 200 400
Number of intervals T e r m i na l r ad i u s ExactConv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL F igure Test 2: Spatially constant solutions – two species.
In the second test, we set again V = . G ( r , = M ( r , = . , and modify the boundary condition accordingly, so that the exact solution is given by (5.1).On the coarsest grid with 50 uniform cells, the numerical solutions at T = . G and M , see Figures 5.4a and 5.4b; comparing the conventional and enhanced methods,however, we see clearly that the enhanced ones produce solutions with much less over-shoots or undershoots. In Figure 5.4d, once more we observe the satisfaction of DTCL bythe enhanced methods, as the incompressibility constraint is well preserved.Similar as the previous test, the final radii computed by all four methods on a sequenceof four meshes are plotted in Figure 5.5; and the numerical errors are summarized in Ta-ble 3–5 for quantitative comparison. For this test, the interaction between the two cellnumbers causes the L -errors in Table 4–5 to be much larger than the previous case; how-ever, the enhanced methods still produce much more accurate solutions than the conven-tional ones. In addition, it is no coincidence that for both enhanced methods, the L -errors VM OF INFILTRATION IN TUMOR GROWTH 27 T able
3. Numerical errors in radius of test 2 at T = . N η Conv. Upwind Conv. MUSCL Enhc. Upwind Enhc. MUSCLError Rate Error Rate Error Rate Error Rate50 -6.25e-3 -2.53e-2 -1.93e-2 -7.71e-3100 -3.79e-3 0.72 -1.27e-2 0.99 -9.75e-3 0.98 -3.94e-3 0.97200 -2.12e-3 0.84 -6.30e-3 1.01 -4.91e-3 0.99 -1.99e-3 0.98400 -1.13e-3 0.91 -3.11e-3 1.02 -2.46e-3 1.00 -1.00e-3 0.99T able L -errors in G ( η, T ) of test 2 at T = . N η Conv. Upwind Conv. MUSCL Enhc. Upwind Enhc. MUSCLError Rate Error Rate Error Rate Error Rate50 1.86e-3 1.57e-3 5.03e-4 1.09e-3100 1.13e-3 0.72 9.79e-4 0.68 2.67e-4 0.91 5.76e-4 0.93200 6.63e-4 0.77 5.83e-4 0.75 1.40e-4 0.94 2.98e-4 0.95400 3.81e-4 0.80 3.39e-4 0.78 7.18e-5 0.96 1.53e-4 0.97T able L -errors in M ( η, T ) of test 2 at T = . N η Conv. Upwind Conv. MUSCL Enhc. Upwind Enhc. MUSCLError Rate Error Rate Error Rate Error Rate50 6.18e-2 1.43e-2 5.03e-4 1.09e-3100 3.55e-2 0.80 6.99e-3 1.03 2.67e-4 0.91 5.76e-4 0.93200 2.01e-2 0.82 3.44e-3 1.02 1.40e-4 0.94 2.98e-4 0.95400 1.12e-2 0.84 1.70e-3 1.02 7.18e-5 0.96 1.53e-4 0.97in M are the same as those in G on the same grids; this is because when DTCL and DGCLare satisfied, the two variables sum up to a constant value, whose numerical error is on thescale of machine precision, as shown in the previous test.5.1.3. Test 3: Monotone growth with constant boundary condition.
In the view of (3.7),we can setup the velocity u and boundary condition M bc accordingly to obtain almost anydesired monotonic growth pattern. To be more specific, suppose a growth curve ˆ R ( t ) withˆ R (cid:48) > u ( ˆ R ( t ) , t ) < , u ( ˆ R ( t ) , t ) M bc ( t ) = − ˆ R (cid:48) ( t ) . Indeed for the simplified model (3.1), if u ( R ( t ) , t ) < t , then the growth of R ( t ) iscompletely determined by the boundary velocity and the boundary condition. In this test,we consider a constant boundary condition M bc ( t ) = .
5, and set up u such that R growslinearly as R ( t ) = R (0) + V t , where V = . u ( r , t ) = − V sin (cid:32) π r R (0) + V t ) (cid:33) . Here u ( r , t ) is nonlinear in space, c.f. the previous test; and we do not expect the solutionsto stay constant across the domain.In Figure 5.6, we plot the numerical solutions for G and M at T = . d θ . All methods predict well the linear growth ofthe radius. Comparing the conventional methods and the enhanced ones, when the formerare used, clear overshoots near the origin and spurious oscillations near the right boundaryin both G and M are observed; however, both enhanced methods seem to lead to smooth Radial coordinate: r G Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( a ) Cell numbers for G at T = . Radial coordinate: r M Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( b ) Cell numbers for M at T = . Time: t T u m o r r ad i u s : R ( t ) ExactConv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( c ) Radius growth history. Time: t L - no r m o f G + M - Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( d ) Histories of d θ . F igure
50 100 200 400
Number of intervals T e r m i na l r ad i u s ExactConv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL F igure d θ in this case is not at the machine precision level for the reason that the proposed methodsare DGCL and DTCL for the interior nodes, whereas our theory does not address whetherit is possible to satisfy these properties with arbitrary incoming data M bc . This is exactlywhat happened here – because of the jump in the boundary condition and the numericalsolution at the last interval, small incompressibility violation is created and propagatedtowards the origin of the domain. Nevertheless, the enhanced methods show significantimprovement over their conventional counterparts.In Figure 5.7 and Table 6 we provide the convergence of the final radii by all methodson the same sequence of grids as well as quantitative comparisons. Clearly, the enhancedmethods provide much more accurate results than the conventional ones.5.1.4. Test 4: A prediction problem with non-monotone radius change.
Finally, we con-sider a test whose radius change cannot be predicted, by considering the velocity:(5.7) u ( r , t ) = V sin( r (1 + t )) , where V = . VM OF INFILTRATION IN TUMOR GROWTH 29 T able
6. Numerical errors in radius of test 3 at T = . N η Conv. Upwind Conv. MUSCL Enhc. Upwind Enhc. MUSCLError Rate Error Rate Error Rate Error Rate50 3.46e-2 1.30e-2 -2.54e-3 -2.08e-3100 1.57e-2 1.14 5.37e-3 1.28 -1.30e-3 0.97 -1.06e-3 0.97200 7.60e-3 1.04 2.58e-3 1.06 -6.59e-4 0.98 -5.33e-4 0.99400 3.77e-3 1.01 1.27e-3 1.02 -3.31e-4 0.99 -2.67e-4 1.00
Radial coordinate: r G Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( a ) Cell numbers for G at T = . Radial coordinate: r M Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( b ) Cell numbers for M at T = . Time: t T u m o r r ad i u s : R ( t ) Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( c ) Radius growth history. Time: t L - no r m o f G + M - Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( d ) Histories of d θ . F igure Time: t L - no r m o f G + M - Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( a ) d θ histories on a 200-interval grid. Time: t L - no r m o f G + M - Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( b ) d θ histories on a 400-interval grid. F igure d θ histories of test 4 on two finer grids.as before, the enhanced methods produce much smoother solutions than the conventionalones, and the overshoots near the origin is in much smaller magnitudes. Figure 5.8d revealsthat the incompressibility constraint is much better preserved by the enhanced methods,whose small violation is due to the boundary conditions.For this test, the overshoot in G near the origin by the conventional methods increasesrapidly as the mesh is refined, and eventually kill the computations when 400 uniformintervals are used to discretize the domain. In Figure 5.9 we plot the d θ histories on a 100-interval grid and a 200-interval grid in the left panel and right panel, respectively. Table 7summarizes the final radii by all methods on the sequence of the four grids; note that we donot compute the numerical error as before since the exact value is unknown. In the table, T able
7. Numerical solutions of the final radius of test 3 at T = . N η Conv. Upwind Conv. MUSCL Enhc. Upwind Enhc. MUSCL50 1.0777 1.0757 1.0713 1.0711100 1.0763 1.0747 1.0735 1.0734200 1.0754 1.0745 1.0745 1.0745400 nan 1.0747 1.0749 1.0749
Radial coordinate: r -0.100.10.20.30.4 G Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( a ) Cell numbers for G at T = . Radial coordinate: r M Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( b ) Cell numbers for M at T = . Time: t T u m o r r ad i u s : R ( t ) Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( c ) Radius growth history. Time: t L - no r m o f G + M - Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( d ) Histories of d θ . F igure t = .
2, c.f. Figure 5.9b; in addition, the conventional MUSCL method seems toproduce non-monotone “convergence” pattern, whereas both enhanced methods seem toprovide monotonic and convergent solutions.5.2.
The tumor growth model.
Next we consider the tumor growth model (2.2) and ourfirst test revisits the case study in Section 2.3. Then, a set of parameters are chosen accord-ing to the study of [21] to assess the impact of using the enhanced methods in practicalpredictions. The PDE system (2.2) involves one more equation for the velocity field U ; forall the four methods including the enhanced ones, we use the same discretization methodas described in Section 2.2 to update A .5.2.1. The case study revisited.
Using the same initial and boundary conditions as in Sec-tion 2.3, the tumor growth model (2.2) is solved by the four methods until T = .
0. Thesample solutions on a grid of 50 uniform intervals are plotted in Figure 5.10. From Fig-ure 5.10d we clearly see that the enhanced methods lead to much smaller incompressibilityviolation than the conventional methods. The numerical solutions show similar pattern asthe simpler model in Section 5.1.3; and we observe alike oscillations in the solutions bythe conventional methods. Nevertheless, the radius growth histories seem to compare wellamong di ff erent methods; and this observation is made more precise by Table 8, whichsummarizes the terminal radii at T = . VM OF INFILTRATION IN TUMOR GROWTH 31 T able
8. Numerical final radius of the tumor case study at T = . N η Conv. Upwind Conv. MUSCL Enhc. Upwind Enhc. MUSCL50 2.1843 2.1748 2.1625 2.1670100 2.1809 2.1737 2.1666 2.1689200 2.1763 2.1718 2.1682 2.1693400 2.1730 2.1705 2.1688 2.1693T able
9. Parameters of the PDGF-driven glioma problem. For the bio-logical interpretation of each parameter, the readers are referred to [21].Parameter Value Dimension λ − µ − δ − ρ − ν · day − m + · ml − · day − β + · mm − γ + − α · ml · day − · pg − θ + · mm − The tumor problem of [21] . Finally, we consider the tumor model describing thePDGF-driven glioma cells in the previous work [21], where the set of parameters are cho-sen so that the survival length matches experimental data. Here the survival length isdefined as the time T term when the radius reaches 5 . mm . The parameters and their dimen-sions are summarized in Table 9.The initial tumor size is given by R (0) = . G ( r , = . θ, H ( r , = . θ, M ( r , = . θ, r ∈ [0 , R (0)) ; A ( r , = − r ) , r ∈ [0 , + ∞ ) . The boundary condition for M describes the environmental number density for the immunecells: M bc ( t ) = . θ . Numerical solutions to this problem on a uniform grid of 50intervals are plotted in the left panel of Figure 5.11, where the solutions of the radius, theglioma cells ( G ), and the total number of cells ( G + H + M ) at T term are plotted from top tobottom. For comparison, the solutions computed on a uniform grid with 200 intervals areprovided side-by-side in the right panel of the same figure. Although the exact solutions tothis problem is unknown, we have the following observations: • The conventional methods seem to underestimate the growth rate of the tumor; andthe modified methods provide much faster convergent results, c.f. Figure 5.11a. • The solutions to cell species are very di ff erent between the conventional methodsand the enhanced ones – particularly near the tumor boundary the conventionalFVMs produce overshoots that grows significantly on finer grids, whereas the en-hanced ones predict a flat plateau, that could possibly represent the “rim” that isreported in many existing studies [4]. See Figure 5.11b. • The incompressibility condition is severely violated by both conventional meth-ods; whereas the enhanced ones respect this constraint very nicely, see Figure 5.11c.
Time (Days) T u m o r r ad i u s ( mm ) Time (Days) T u m o r r ad i u s ( mm ) ( a ) Radius growth histories. Radial coordinate: r (mm) G , m u l DH ( c e ll ) Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL 0 1 2 3 4 5
Radial coordinate: r (mm) G , m u l DH ( c e ll ) Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( b ) Glioma cell distributions at T term . Radial coordinate: r (mm) G + H + M ( c e ll ) Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL 0 1 2 3 4 5
Radial coordinate: r (mm) G + H + M ( c e ll ) Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( c ) Total cell distributions at T term . F igure able
10. Survival length (day) of the PDGF-driven model. N η Conv. Upwind Conv. MUSCL Enhc. Upwind Enhc. MUSCL50 48.8644 53.0510 58.5284 58.5276100 50.4446 54.5993 58.4433 58.4408200 54.0926 57.1748 58.4223 58.4212400 57.6348 58.4293 58.4177 58.4174Quantitative comparisons are provided in Table 10, which summarizes the survival length T term computed by all four methods on a sequence of four grids. This table reveals thatalthough the conventional methods predict T term that is somewhat di ff erent from that bythe enhanced methods, the predictions converge nevertheless to the same value as the gridis refined. Hence we claim that the enhanced methods indeed improve the accuracy of thenumerical simulation.Finally, it is worth noting that although the conventional methods compute very di ff erentsolutions in the cell numbers, they seem to give reasonable predictions on the tumor growthcurves, which explains why the numerical results compare reasonably well to experimentaldata in the previous work [21]. To close this section, we provide an explanation to thephenomenon by utilizing a relation that is analogous to (3.7). Particularly, when the tumorgrows monotonically as in the present case, a simpler formula determining the growth VM OF INFILTRATION IN TUMOR GROWTH 33
Time (Days) - ∂ A / ∂ r a t R ( t ) Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( a ) 50-interval grid. Time (Days) - ∂ A / ∂ r a t R ( t ) Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( b ) 200-interval grid. F igure ∂ A ( R ( t ) , t ) ∂ r on two uniform grids.pattern is: R (cid:48) ( t ) = − α ∂ A ( R ( t ) , t ) ∂ r M bc ( t ) . Since M bc is a constant, the growth is determined by the chemoattractant concentrationgradient at the tumor boundary. Because A is governed by the di ff usion-reaction equa-tion (2.2j), its profile is less a ff ected by di ff erent cell number solutions in G . Particularly,we plot numerical solutions to ∂ A ( R ( t ) , t ) ∂ r by various methods in Figure 5.12, and see thatthe di ff erence between the conventional methods and enhanced methods is less significantthan the di ff erence in the cell number solutions.6. C onclusions We propose a finite volume framework with segregated fluxes for numerical computa-tion of free boundary problems that model infiltration dynamics in spherically symmetrictumor growth. Under this framework, su ffi cient conditions for ensuring the geometricconservation law on a moving grid and the incompressibility constraint are derived; andclassical first-order and second-order finite volume methods are enhanced following theserequirements. The numerical performance of the enhanced methods are assessed by sev-eral representative tests, either for a simplified model or a full PDGF-driven tumor growthmodel; and their solutions exhibit significant improvements over those by conventionalmethods. More importantly, the cell-incompressibility condition is well respected by theenhanced methods but not by the conventional one; and it is shown to be crucial to deliverconvergent and stable solutions on refined grids.It worth noting that, although the MUSCL-type methods generally produce more ac-curate solutions than the upwind ones, they do not deliver second-order convergence evenwhen the solutions are smooth. This is probably due to the integro-di ff erential nature ofthe governing equation; and how to improve the second-order methods will be addressedin future work. Nevertheless, the methodology to ensure the DGCL and DTCL propertiesfor MUSCL-based methods is expected to remain the same; hence it is addressed in thecurrent paper instead of in future publications.A cknowledgements X. Zeng would like to thank University of Texas at El Paso for the start up support.P. Tian would like to thank the National Science Foundation of US for the support inmathematical modeling under the grant number DMS-1446139. R eferences [1] Martin Burger, Marco Di Francesco, and Yasmin Dolak-Struss. The KellerSegel model for chemotaxis withprevention of overcrowding: Linear vs. nonlinear di ff usion. SIAM J. Math. Anal., 38(4):1288–1315, 2006.[2] J. C. Butcher. Numerical Methods for Ordinary Di ff erential Equations. John Wiley & Sons, 3rd edition,2016.[3] Vincent Calvez and Jos´e Carrillo. Volume e ff ects in the Keller-Segel model: energy estimates preventingblow-up. J. Math. Pure. Appl., 86(2):155–175, August 2006.[4] J. J. Casciari, S. V. Sotirchos, and R. M. Sutherland. Mathematical modelling of microenvironment andgrowth in EMT6 / Ro multicellular tumour spheroids. Cell Proliferat., 25(1):1–22, January 1992.[5] Xinfu Chen and Avner Friedman. A free boundary problem for an elliptic-hyperbolic system: An applicationto tumor growth. SIAM J. Math. Anal., 35(4):974–986, 2003.[6] Alina Chertock, Yekaterina Epshteyn, Hengrui Hu, and Alexander Kurganov. High-order positivity-preserving hybrid finite-volume-finite-di ff erence methods for chemotaxis systems. Adv. Comput. Math.,44(1):327–350, February 2018.[7] Phillip Colella and Paul R. Woodward. The piecewise parabolic method (ppm) for gas-dynamical simula-tions. J. Comput. Phys., 54(1):174–201, April 1984.[8] Shangbin Cui and Avner Friedman. A hyperbolic free boundary problem modeling tumor growth. InterfaceFree Bound., 5(2):159–182, 2003.[9] Elio Espejo, Karina Vilches, and Carlos Conca. Sharp condition for blow-up and global existence in a twospecies chemotactic Keller-Segel system in R . Eur. J. Appl. Math., 24(2):297–313, April 2013.[10] Charbel Farhat, Philippe Geuzaine, and C´eline Grandmont. The discrete geometric conservation law andthe nonlinear stability of ale schemes for the solution of flow problems on moving grids. J. Comput. Phys.,174(2):669–694, December 2001.[11] Francis Filbet. A finite volume scheme for the Patlak-Keller-Segel chemotaxis model. Numer. Math.,104(4):457–488, October 2006.[12] Avner Friedman, Wenrui Hao, and Bei Hu. A free boundary problem for steady small plaques in the arteryand their stability. J. Di ff er. Equations, 259(4):1227–1255, August 2015.[13] Avner Friedman, Bei Hu, and Chuan Xue. Analysis of a mathematical model of ischemic cutaneous wounds.SIAM J. Math. Anal., 42(5):2013–2040, 2010.[14] Sigal Gottlieb and Chi-Wang Shu. Total variational diminishing runge-kutta schemes. Math. Comput.,67(221):73–85, January 1998.[15] Wenrui Hao and Avner Friedman. The LDL–HDL profile determines the risk of atherosclerosis: A mathe-matical model. PLoS ONE, 9(3):e90497, March 2014.[16] Wenrui Hao, Larry S. Schlesinger, and Avner Friedman. Modeling granulomas in response to infection inthe lung. PLoS ONE, 11(3):e014738, March 2016.[17] Evelyn F. Keller and Lee A. Segel. Model for chemotaxis. J. Theor. Biol., 30(2):225–234, February 1971.[18] Inwon Kim and Yao Yao. The Patlak-Keller-Segel model and its variations: Properties of solutions viamaximum principle. SIAM J. Math. Anal., 44(2):568–602, 2012.[19] Randall LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge Texts in Applied Mathe-matics. Cambridge University Press, 2002.[20] K. W. Morton and P. K. Sweby. A comparison of flux limited di ff erence methods and characteristic galerkinmethods for shock modelling. J. Comput. Phys., 73(1):203–230, November 1987.[21] Ben Niu, Xianyi Zeng, Frank Szulzewsky, Sarah Holte, Philip K. Maini, Eric C. Holland, and Jianjun PaulTian. Mathematical modeling of PDGF-driven glioma reveals the infiltrating dynamics of immune cells intotumors. 2018. Submitted.[22] A. L´opez Ortega and G. Scovazzi. A geometrically-conservative, synchronized, flux-corrected remap forarbitrary lagrangian-eulerian computations with nodal finite elements. J. Comput. Phys., 230(17):6709–6741, July 2011.[23] Benoˆıt Perthame and Anne-Laure Dalibard. Existence of solutions of the hyperbolic Keller-Segel model. T.Am. Math. Soc., 361(5):2319–2335, May 2009.[24] Chi-Wang Shu. Total-variation-diminishing time discretization. SIAM J. Sci. Stat. Comp., 9(6):1073–1084,November 1988.[25] Bram van Leer. Towards the ultimate conservative di ff erence scheme V. A second-order sequel to Godunov’smethod. J. Comput. Phys., 32(1):101–136, July 1979.[26] Xianyi Zeng. A general approach to enhance slope limiters in MUSCL schemes on nonuniform rectilineargrids. SIAM J. Sci. Comput., 38(2):A789–A813, 2016. VM OF INFILTRATION IN TUMOR GROWTH 35 A ppendix A. S plitted velocities for advection equations
In Section 4, we split the advection velocity for species and compute the fluxes sep-arately. In this section, we briefly investigate the stability associated with the splittingstrategy by the classical von Neumann analysis. Since i is used to denote the imaginaryunit, we’ll use j to denote the grid index. Let us consider the following one-dimensionaladvection equation:(A.1) ∂ G ∂τ + V ∂ G ∂η = , where G is the advected quantity and V is the constant advection velocity. Following thesplitting strategy, we rewrite V as the sum of K constants:(A.2) V = V + V + · · · + V K and solve the corresponding equation by the first-order upwind fluxes and first-order for-ward Euler time-integrator:(A.3) G n + j − / − G nj − / ∆ τ + K (cid:88) k = F upw ( V k ; G nj − / , G nj + / ) − F upw ( V k ; G nj − / , G nj − / ) ∆ η = , where ∆ τ is the time step size and F upw is given by (2.8).Clearly, the fluxes collapse into two groups, namely those associated with positivevelocities and those associated with negative ones. Denoting V + = (cid:80) ≤ k ≤ K , V k > V k and V − = (cid:80) ≤ k ≤ K , V k < V k , (A.3) simplifies to:(A.4) G n + j − / − G nj − / ∆ τ + V + ( G nj − / − G nj − / ) ∆ η + V − ( G nj + / − G nj − / ) ∆ η = . Following the standard von Neumann analysis, we write:(A.5) G nj − / = a n e i j κ ∆ η , where a is the so called amplifier coe ffi cient and κ is the arbitrary wave number; the nu-merical method is stable if and only if there is a ∆ τ c >
0, such that for all 0 ≤ ∆ τ ≤ ∆ τ c wehave | a | ≤ κ ∈ R .Denoting θ = κ ∆ η for simplicity, plugging (A.5) into (A.4) we have: a n + e i j θ − a n e i j θ ∆ τ + V + ∆ η ( a n e i j θ − a n e i ( j − θ ) + V − ∆ η ( a n e i ( j + θ − a n e i j θ ) = a = − V + ∆ τ ∆ η (1 − e − i θ ) − V − ∆ τ ∆ η ( e i θ − . Let the Courant numbers corresponding to V + and V − be α + = V + ∆ τ/ ∆ η ≥ α − = − V − ∆ τ/ ∆ η ≥
0, respectively, it is easy to compute that: a = − ( α + + α − )(1 − cos θ ) + i ( α + − α − ) sin θ ⇒| a | = (1 − ( α + + α − )(1 − cos θ )) + (( α + − α − ) sin θ ) = (1 − ( α + + α − )(1 − cos θ )) + (( α + + α − ) sin θ ) − α + α − sin θ = − α + + α − )(1 − α + − α − )(1 − cos θ ) − α + α − sin θ = − [2( α + + α − )(1 − α + − α − ) + α + α − (1 + cos θ )] (1 − cos θ ) ;it follows that | a | ≤ θ ∈ R if and only if for these θ :2( α + + α − )(1 − α + − α − ) + α + α − (1 + cos θ ) ≥ , or equivalently, α + + α − ≤
1. Hence the explicit split method is conditionally stable and theCourant condition is:(A.6) ( | V + | + | V − | ) ∆ τ ≤ ∆ η or K (cid:88) k = | V k | ∆ τ ≤ ∆ η . Similarly, using the implicit first-order backward Euler time-integrator instead, onecomputes: a = + α + (1 − e − i θ ) − α − ( e i θ − = + ( α + + α − )(1 − cos θ ) + i ( α + − α − ) sin θ , and | a | ≤ + ( α + + α − )(1 − cos θ )) + ( α + − α − ) sin θ ≥ , which holds naturally for all θ . To conclude, the implicit split method is unconditionallystable. A ppendix B. I mplicit E nhanced F inite V olume M ethods In this appendix we extend the enhanced method of Section 4 to implicit time-integrators.First, let us consider the backward Euler time-integrator, which is unconditionally stablecombined with the segregated upwind flux, as shown in the previous appendix.To illustrate the idea, in order to update the solutions from τ n to τ n + , all spatial dis-cretizations happen at τ n + instead of τ n , c.f. the explicit methods. For example, thecounterpart of (4.9) reads: η j − / [( R n + ) G n + j − / − ( R n ) G nj − / ] ∆ τ n + ∆ η (cid:104) F G , n + j − F G , n + j − (cid:105) (B.1) = η j − / ( R n + ) f n + j − / − η j − / R (cid:48) n + R n + G n + j − / , where F G , n + j approximates:(B.2) F G , n + j ≈ (cid:32) VR − η R (cid:48) R (cid:33) η R G (cid:12)(cid:12)(cid:12)(cid:12) η = η j , τ = τ n + , and it is segrated into F G , n + j = F G , n + V , j + F G , n + R (cid:48) , j .In Section 4, (4.2) is obtained from the explicit formula (4.1); hence it needs to bemodified to:(B.3) R (cid:48) n + = ( R n + ) − ( R n ) ∆ τ n R n + . In analogous of (4.4), we compute R (cid:48) n + such that:(B.4) R (cid:48) n + = (cid:32) − ∆ η (cid:33) − V n + N η . Now we have similar to Theorem 4.6 the following result:
Theorem B.1.
The numerical method given by the implicit version of (4.9), (4.10) and(4.15), c.f., (B.1), satisfies both DGCL and DTCL if: (1) F n + j is additive and V-consistent,(2) ˆ F n + j is additive and cubic-preserving, (3) F M , n + u , j = F M , n + u , j , and (4) R (cid:48) n + equals theright hand side of (B.3). VM OF INFILTRATION IN TUMOR GROWTH 37
Radial coordinate: r M Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( a ) Cell numbers for M at T = . Time: t T u m o r r ad i u s : R ( t ) ExactConv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( b ) Radius growth history. Time: t L - no r m o f G + M - Conv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( c ) Histories of d θ .
50 100 200 400
Number of intervals T e r m i na l r ad i u s ExactConv. UpwindConv. MUSCLEnhc. UpwindEnhc. MUSCL ( d ) Final radii convergence. F igure B.1. Implicit solutions to the test in Section 5.1.1 on a 50-cell grid.The proof is completely analogous and omitted here.Extension to higher-order implicit time-integrators are straightforward by using the Di-agonally Implicit Runge-Kutta (DIRK) methods, which are well documented in many textson numericla methods for ordinary di ff erential equations, such as [2]. In essense, a DIRKmethod is a multi-stage method with higher time accuracy, where each stage is equivalentto a backward Euler step; hence the method described before extends naturally to thesetime-integrators. The particular one that we will use in combine with MUSCL in space isthe second-order DIRK method given in Section 361 of [2].To demonstrate the numerical performances, we repeat the test in Section 5.1.1 using amuch larger Courant number α cfl = .
0. The solutions on a grid of 50 uniform cells as wellas the convergence plots of the terminal radii are plotted in Figure B.1. Comparing thesefigures to Figures 5.1–5.3, we can draw very similar conclusions, indicating the successfulextension of the enhanced methods in Section 4 to implicit time-integrators. D epartment of M athematical S ciences ,, C omputational S cience P rogram , U niversity of T exas at E l P aso ,E l P aso , TX 79902, U nited S tates ., T el .: + E-mail address , Corresponding author, X. Zeng: [email protected] C omputational S cience P rogram , U niversity of T exas at E l P aso , E l P aso , TX 79902, U nited S tates . E-mail address , M. Saleh: [email protected] D epartment of M athematical S cience , N ew M exico S tate U niversity , L as C ruces , NM 88003, U nited S tates . E-mail address , J. Tian:, J. Tian: