MATHEMATICAL MODEL AND ITS FAST NUMERICAL METHOD FOR THE TUMOR GROWTH

In this work, we reformulate the diffuse interface model of the tumor growth (S.M. Wise et al., Three-dimensional multispecies nonlinear tumor growth-I: model and numerical method, J. Theor. Biol. 253 (2008) 524–543). In the new proposed model, we use the conservative secondorder Allen–Cahn equation with a space–time dependent Lagrange multiplier instead of using the fourth-order Cahn–Hilliard equation in the original model. To numerically solve the new model, we apply a recently developed hybrid numerical method. We perform various numerical experiments. The computational results demonstrate that the new model is not only fast but also has a good feature such as distributing excess mass from the inside of tumor to its boundary regions.

1. Introduction.The morphological evolution of a growing tumor is the result of many factors, such as cell-cell and cell-matrix adhesion, mechanical stress, cell motility, and transport of oxygen, nutrients, and growth factors [3].Mathematical modeling of cancer gives unique and important insights into tumor progression, helps explain experimental and clinical observations, and helps provide optimal treatment strategies [24].In the past several years, a considerable amount of research on mathematical models of cancer has been conducted, and numerical simulations of tumor growth have been performed [6,10,15,25,26,29,31,34,36,38,40,41,47,49,52,54,55,56,59,63].A variety of modeling strategies are available to investigate one or more aspects of cancer.Discrete cell-based models (e.g., cellular automata [2,9,30,44,46,62] and agent-based models [7,43,50]), where individual cells are tracked and updated according to a specific set of biophysical rules, are particularly useful for studying carcinogenesis, natural selection, genetic instability, and interactions of individual cells with each other and the microenvironment.In larger-scale systems, continuum methods provide a good modeling alternative [14,16,17,18,19,27,37,39,48].The governing equations are typically of the reaction-diffusion type.Recently, nonlinear continuum modeling has been 1174 HYUN GEUN LEE, YANGJIN KIM AND JUNSEOK KIM performed to study the effects of shape instabilities on avascular, vascular, and angiogenic solid tumor growth.Cristini et al. [27] performed the first fully nonlinear simulations of a continuum model of avascular and vascularized tumor growth in two dimensions using a boundary integral method.Li et al. [48] extended this model to three dimensions using an adaptive boundary integral method.Zheng et al. [65] also extended this model to include a hybrid continuum discrete model of angiogenesis (based on earlier work of Anderson and Chaplain [5]) and investigated the nonlinear coupling between growth and angiogenesis in two dimensions using finite element/level-set method.Wise et al. simulated tumor growth [64] and angiogenesis [35] in three dimensions using a diffuse interface, multiphase mixture model.Chen et al. [24] extended the model of Wise et al. and incorporated the effect of a stiff membrane to model tumor growth in a confined microenvironment.
Here, we reformulate the diffuse interface continuum model of multispecies tumor growth of Wise et al. [64].The model consists of fourth-order nonlinear advection-reaction-diffusion equations of Cahn-Hilliard-type (CH) [21] for the cell species volume fractions coupled with reaction-diffusion equations for the substrate components.Because the original model involves fourth-order equations, it is challenging to develop accurate and efficient numerical methods.For example, explicit methods suffer from severe time step restrictions (∆t < C(∆x) 4 ) and thus are computationally expensive to handle large systems.Therefore, in the new proposed model, we use the conservative second-order Allen-Cahn (AC) equation with a space-time dependent Lagrange multiplier instead of using the fourth-order CH equation in the original model.The classical AC equation was originally introduced as a phenomenological model for antiphase domain coarsening in a binary alloy [4], but does not conserve mass of the mixture.Rubinstein and Sternberg [57] introduced a nonlocal AC equation with a time dependent Lagrange multiplier to enforce conservation of mass.However, with their model, it is difficult to keep small features since they dissolve into the bulk region.One of the reasons for this is that mass conservation is realized by a global correction using the time-dependent Lagrange multiplier.To resolve the problem, we use a space-time dependent Lagrange multiplier to preserve the mass of the mixture.And, to numerically solve the new model, we use a recently developed hybrid numerical method [45].
This paper is organized as follows.In Section 2, we reformulate the diffuse interface model of Wise et al. [64] by using the conservative second-order AC equation with a space-time dependent Lagrange multiplier.A numerical algorithm using an operator splitting method is described in Section 3. Various numerical experiments are presented in Section 4. Finally, conclusions are drawn in Section 5.
2. Mathematical model.In this section, we present a mathematical model of tumor growth.We begin by recalling the nondimensional tumor growth model from Wise et al. [64], where a thermodynamically consistent diffuse interface continuum model of multispecies tumor growth was developed, analyzed, and simulated.The authors take into account mechanical interactions, mainly focused on cell-cell adhesion between a tumor and host.The dimensionless dependent variables defined in a bounded tissue domain Ω ⊂ R d (d = 2 or 3) are as follows: φ, ψ, and ξ are the volume fractions of tumor cells, dead tumor cells, and host tissue, respectively (see Figure 1).Here, we note that φ − ψ is the volume fraction of viable tumor cells.Hence, φ is the sum of the volume fractions of viable and dead tumor cells.u, p, and n are the tumor velocity, pressure, and nutrient concentration, respectively.The original governing equations for the tumor growth in [64] are where M > 0 is a mobility, µ is the chemical potential, F (φ) = φ 2 (1 − φ) 2 /4 is a double well bulk energy, > 0 is a parameter related to the thickness of the diffuse interface that separates the tumor and host domains, and γ is a parameter related to the cell-cell adhesion force.S T and S D are the net sources of tumor and dead cells, respectively, and are defined as where λ M , λ L , λ A , and λ N are the rates of volume gain or loss due to cellular mitosis, lysing, apoptosis, and necrosis, respectively, H is a Heaviside step function, and n N is the nutrient limit for cell viability.The diffusion coefficient D(φ) and nutrient capillary source term T C (φ, n) are, respectively, where D H is the nondimensional nutrient diffusion coefficient in the host domain, ν H p and ν T p denote the nutrient transfer rates for preexisting vascularization in the host and tumor domains, respectively, and n c is the nutrient level in the capillaries.Q(φ) is an interpolation function and is defined as , and Q(0) = 0. ν U is the nutrient uptake rate by the viable tumor cells.Eqs. ( 1)-( 6) are valid on the whole domain Ω and not just on the tumor volume.There are no boundary conditions required for φ and ψ at the tumor boundary.At the outer boundary, we choose the following boundary conditions where n is the unit normal vector to ∂Ω.
The diffuse interface model of Wise et al. [64] involves the fourth-order CH Eq. ( 1) and ( 2) with a source term.In this paper, we propose an alternative model for Eqs. ( 1) and ( 2).The proposed model consists of the following two equations and Eqs. ( 3)-( 6): First, we update φ according to Eq. ( 7), and then we relax φ using Eq. ( 8).In Eq. ( 8), is the classical AC equation which was originally introduced as a phenomenological model for antiphase domain coarsening in a binary alloy [4].Since the classical AC equation does not have the mass conservation property, Brassel and Bretin [11] introduced a nonlocal AC equation with a space-time dependent Lagrange multiplier (β(t) F (φ)) to enforce conservation of mass.Here, β(t) satisfies β(t) = Ω F (φ) dx/ Ω F (φ) dx.The proposed model involves a second-order equation and we will apply the recently developed hybrid numerical method [45] to numerically solve it.
3. Numerical solution.In this section, we describe an operator splitting algorithm for solving Eqs. ( 3)- (8).For simplicity and clarity of exposition, we shall discretize Eqs. ( 3)-( 8) in two-dimensional space, i.e., Ω = (a, b) × (c, d).Threedimensional discretization is defined analogously.Let the computational domain be partitioned into a uniform mesh with mesh spacing h.The center of each cell, Ω ij , is located at (x i , y j ) = ((i − 0.5)h, (j − 0.5)h) for i = 1, . . ., N x and j = 1, . . ., N y .N x and N y denote the number of cells in the x-and y-directions, respectively.Cell vertices are located at (x i+ 1 2 , y j+ 1 2 ) = (ih, jh).In this paper, tumor and dead cells, pressures, and nutrients are stored at the cell centers and velocities at cell faces.Let φ k ij be the approximations of φ(x i , y j , k∆t), where ∆t = T /N t is the time step, T is the final time, and N t is the total number of time steps.The other terms are similarly defined.
In this paper, we use an operator splitting method, in which we numerically solve Eqs. ( 7) and ( 8) by solving successively a sequence of simpler problems: First, we solve Eq. ( 9) by applying the explicit Euler's method: where ∇ d • is the discrete divergence operator.Next, we solve Eq. ( 10) by applying the explicit Euler's method: where M k = M φ k ij and ∆ d is the discrete Laplacian operator.And Eq. ( 11) is solved analytically using the method of separation of variables [60] and the solution is given as .
Finally, we discretize Eq. ( 12) as By Eq. ( 13), we get ), then by the property of mass conservation Thus, Eqs. ( 3)-( 6) are discretized as where u and v are the horizontal and vertical components of u, respectively.The discrete differentiation operators are and ∇ d is the discrete gradient operator.Apply the divergence operator to Eqs. ( 14) and ( 15) and get a Poisson equation for p k+1 ij : The resulting linear systems of Eqs. ( 16) and ( 17) are solved by a fast solver, such as a linear multigrid method [12,61].
4.1.Time scaling between the Cahn-Hilliard and conservative Allen-Cahn models.In this paper, we use the conservative second-order Allen-Cahn (CAC) equation with a space-time dependent Lagrange multiplier instead of using the fourth-order Cahn-Hilliard (CH) equation in the original model.The CAC, constant mobility CH, and variable mobility CH equations provide an approximation to motion by the volume preserving mean curvature flow [8,11,13,42,58], the Mullins-Sekerka flow [1,22,23,28,33,53], and the surface diffusion flow [20,32,51], respectively.Thus, there is a need for a time scaling to consider a difference between the motion of the interface for the CAC, constant mobility CH, and variable mobility CH equations.To evaluate a time scaling between the CH (Eqs.( 1) and ( 2)) and CAC (Eqs.( 7) and ( 8)) models, we consider the following initial condition: 20], with h = 20/128, ∆t = 0.01, and = 0.1 √ 2. In this test, the effects of velocity u and net source of tumor cells S T are negligible and we consider a constant mobility case.The numerical solution is computed to time T = 200.
Figures 2 (a) and (b) show the time evolutions of y = 10 of 0.5-level of φ obtained by solving the CH and CAC equations without and with time scaling, respectively.Here, a time scaling factor is about 3.2, that is, M CAC = 3.2M CH , where M CH and M CAC are constant mobilities of the CH and CAC equations, respectively.Note that the time scaling factor depends on the initial morphology of the interface.In this simulation, we solve Eqs. ( 1)-( 6) (the CH model) and ( 3)-( 8) (the CAC model) with the following parameters: M = 10 for the CH model, M = 32 for the CAC model, γ = 0.0, ν U = 1.0, λ M = 8.0, λ L = 1.0, λ A = 0.0, λ N = 3.0, n N = 0.6, D H = 1.0 × 10 3 , ν H p = 0.0, ν T p = 0.0, and n c = 1.0.Note that we take the same parameter values as in [64] except for λ M .To investigate the difference in distributing excess mass, we choose 8 times larger than the value in [64].
Figures 3 and 4 show the time evolutions of tumor cells obtained by solving the CH and CAC models, respectively.As we can see in Figure 3, the CH model does not distribute well excess mass from the inside of tumor to its boundary regions and thus excess mass builds up inside and the volume fraction of tumor cells becomes much larger than one.On the other hand, in the CAC model, excess mass (obtained by solving Eq. ( 9)) diffuses according to Eq. ( 10) and then the volume fraction of tumor cells becomes close to one according to Eq. ( 11).Finally, mass (changed by the AC step ( 10) and ( 11)) is corrected in the interfacial region by the mass correction step (12).Therefore, excess mass is well distributed to boundary regions of tumor cells.In this section, we present two-dimensional simulations of tumor growth using the parameters given in Table 1.Because the diffusivity of the nutrient in the host medium is 1000 times larger than that in the tumor interstitium, we use D(φ) as in [64] To validate our new model and numerical algorithm, we take the same initial condition (18)   Table 1.Nondimensional parameters used in the two-dimensional numerical simulations.
Figure 5 shows the evolution of the tumor.Initially, there are no dead cells in the tumor.But, as the nutrient concentration falls below level needed for viability, dead cells quickly begin to accrue.At time t = 5, the tumor has a fully developed necrotic core.One can observe a slight bulge oriented along the x-direction.At later times, the instability becomes more pronounced and the tumor develops buds that elongate into protruding fingers.The instability enables the tumor to increase its exposure to nutrient as its surface area increases relative to its volume.This allows the tumor to overcome the diffusional limitations to growth.The tumor will grow indefinitely as the instability repeats itself on the buds and protruding fingers.This result is in good agreement with the result in [64].Furthermore, during the simulation, our new model ( 3)-( 8) and numerical method take 500.624sCPU time, whereas the CH-type model ( 1)-( 6) and numerical method in [24] take 3034.718sCPU time.Since our new model involves a second-order equation and we use the recently developed hybrid numerical method for solving the model, we save significant computational time.In the following subsections, we will investigate the effect of biophysical parameters given in Table 1.The viable tumor cells are primarily contained between the inner and outer contours.The biophysical parameters are given in Table 1.

4.3.1.
Effect of λ L .In the net sources S T and S D of tumor and dead cells, λ L is the rate of volume loss due to cellular lysing.Thus, more tumor and dead cells are lysed as λ L increases, and the more lysed tumor and dead cells translate to the host tissue.To investigate the effect of λ L , we take the same initial condition (18) and parameter values used to create Figure 5 except for λ L .We vary λ L as λ L = 0.4, 3.0 and λ L = 1.0 (see Figure 5).Figure 6 shows the time evolutions of the tumor for different λ L values.From the results, we observe that the growth of tumor is inhibited as λ L increases.
where n N is the nutrient limit for cell viability.When n falls below n N , tumor cells are dead proportional to λ N .Thus, more tumor cells translate to dead cells as n N decreases.To investigate the effect of n N , we take the same initial condition (18) and parameter values used to create Figure 5 except for n N .We vary n N as n N = 0.4, 0.8, 0.9 and n N = 0.6 (see Figure 5).Figure 7 shows the time evolutions of the tumor for different n N values.As we expected, we can see more and more dead cells as n N decreases.

Effect of ν H
p and ν T p .The nutrient capillary source term T C (φ, n) is defined as , where ν H p and ν T p denote the nutrient transfer rates for preexisting vascularization in the host and tumor domains, respectively.As ν H p and ν T p increase, more nutrients are supplied in the host and tumor domains.In this test, we take the same initial condition (18) and parameter values used to create Figure 5 except for ν H p and ν T p . Figure 8 shows the evolution of the tumor with ν H p = ν T p = 0.4.Since more nutrients are supplied in the host and tumor domains, the tumor grows bigger.We choose the parameter values used to create Figure 5. Figure 9 shows the time evolutions of the tumor for different r values.Depending on the initial condition, we can see various patterns of tumor growth.4.4.Three-dimensional tumor growth.In this paper, we also perform a threedimensional simulation of tumor growth using the parameters given in Table 1 Figure 10 shows the evolution of the tumor.This simulation demonstrates the capability of feasibly simulating complex tumor progression in three dimensions.5. Conclusions.In this paper, we reformulated the diffuse interface model of the tumor growth of Wise et al. [64].In the new proposed model, we used the conservative second-order Allen-Cahn equation with a space-time dependent Lagrange multiplier instead of using the fourth-order Cahn-Hilliard equation in the original model.To numerically solve the model, we applied a recently developed hybrid numerical method.Through numerical examples, we observed that the new model is not only fast but also has a good feature such as distributing excess mass from the inside of tumor to its boundary regions.We also performed various numerical experiments varying the biophysical parameters.

Figure 1 .
Figure 1.Schematic of a growing tumor.φ and ψ are the volume fractions of the tumor and dead tumor cells, respectively.ξ is the volume fraction of the host tissue.

Figure 2 .
Figure 2. Time evolutions of y = 10 of 0.5-level of φ obtained by solving the CH and CAC equations.

Figure 5 .
Figure 5. Evolution of the contours φ − ψ = 0.5 during growth.The viable tumor cells are primarily contained between the inner and outer contours.The biophysical parameters are given in Table1.

Figure 10 .
Figure 10.Evolution of the contours φ − ψ = 0.5 during growth.The viable cells are primarily contained between the inner and outer surfaces.