Numerical resolution of a model of tumour growth.

Mathematical Medicine and Biology Advance Access published February 11, 2015 Mathematical Medicine and Biology (2015) Page 1 of 29 doi:10.1093/imammb/dqv004

Numerical resolution of a model of tumour growth

[Received on 3 December 2013; revised on 30 October 2014; accepted on 21 January 2015] We consider and solve numerically a mathematical model of tumour growth based on cancer stem cells (CSC) hypothesis with the aim of gaining some insight into the relation of different processes leading to exponential growth in solid tumours and into the evolution of different subpopulations of cells. The model consists of four hyperbolic equations of first order to describe the evolution of four subpopulations of cells. A fifth equation is introduced to model the evolution of the moving boundary. The coefficients of the model represent the rates at which reactions occur. In order to integrate numerically the four hyperbolic equations, a formulation in terms of the total derivatives is posed. A finite element discretization is applied to integrate the model equations in space. Our numerical results suggest the existence of a pseudo-equilibrium state reached at the early stage of the tumour, for which the fraction of CSC remains small. We include the study of the behaviour of the solutions for longer times and we obtain that the solutions to the system of partial differential equations stabilize to homogeneous steady states whose values depend only on the values of the parameters. We show that CSC may comprise different proportions of the tumour, becoming, in some cases, the predominant type of cells within the tumour. We also obtain that possible effective measure to detain tumour progression should combine the targeting of CSC with the targeting of progenitor cells. Keywords: tumour evolution; cancer stem cells; free boundary problems; numerical resolution; finite element methods.

1. Introduction Mathematical modelling might be highly useful for simulating the dynamics of cancer initiation and progression. During the last decades systems of partial differential equations (PDEs) to model tumour growth have been studied and experiments have confirmed the existence of subpopulations of cancer stem cells (CSC) inside the tumours of most cancer’s types. The model studied here was presented and mathematically studied in Tello (2013). The model is based on the CSC hypothesis which basically postulates that most tumours arise by consecutive genetics changes in a subpopulation of cells that have intrinsics characteristics similar to those of normal cells. Distinct properties of these so-called CSC are the longevity or even immortality, the self-renewal capacity, the unlimited proliferation and the ability to produce more such CSC as well as non-stem cancer cells (CC) (see, e.g. Hillen et al., 2013). Experimental studies evidence CSC as responsible for the long-term survival of some type of cancer after therapies, while other experiments are focused on the role of CSC in metastatic progression of cancer (see Collins et al., 2005; Dylla et al., 2008), nevertheless the knowledge about these cells is still limited. CSC’s mitosis may originate two CSC or two progenitor cells through symmetric division or one of each class through asymmetric division. Regulation of symmetric or asymmetric division is a complex process which depends on a range of conditions (see Morrison & Kimble, 2006; Boman c The authors 2015. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

Downloaded from http://imammb.oxfordjournals.org/ at Florida International University on November 14, 2015

Ana I. Muñoz Department of Applied Mathematics, Rey Juan Carlos University, Tulipán Street, Móstoles, Madrid, E 28933, Spain Corresponding author: [email protected]

2 of 29

A. I. MUÑOZ ET AL.


et al., 2007; Dick, 2008; Turner et al., 2009). The regulation process still possesses several steps not well understood. Reports have suggested that regulation of the ratio between symmetric and asymmetric division might be crucial for the development and progression of cancer (see Dingly & Michor, 2006, Caussinus & Hirth (2007); Giebel & Bruns (2008); Pacheco et al., 2009). An effective therapeutic strategy for fighting against cancer relates to targeting of CSC killing them (see Collins et al., 2005; Dingly & Michor, 2006; Ganguly & Puri, 2007; Maitland & Collins, 2008). The need to not only specifically target CSC, but also the progenitor subpopulation to effectively detain tumour progression has been suggested before in literature (see Kern & Shibata, 2007). Cancer therapy cannot be successful unless it eliminates the CSC. Thus, understanding the nature and the evolution of CSC is crucially important in the development of effective therapies. Characterization of tumour growth kinetics in terms of clinically relevant parameters is increasingly required for optimizing and personalizing treatments. The results of designing a numerical model that is robust in the parameter space of realistic physiological scenarios, and further validated with accurate experimental measurements, suggests its potential as a tool for analysing different therapy strategies against tumours. Modelling cancer tumour evolution is a complex issue that has been explored in recent years resulting in the production of numerous mathematical models with different approaches. Some tumour growth models represent the tumour as a mass of discrete inner layers separated by moving boundaries, and other models of tumour growth are formulated in terms of a continuum of proliferating, quiescent or necrotic cells (see Bellomo & Preziosi, 2000; Friedman, 2004, 2006, 2007; Byrne, 2009). Some numerical simulations of models of tumour growth can be found in Valenciano & Chaplain (2003), Zheng et al. (2005), Barrea & Turner (2005), Wang & Tian (2008), Chaplain (2008) and Hao et al. (2012). For example, Valenciano & Chaplain (2003) employed a spectral element method in order to obtain high accurate solutions of a tumour angiogenesis model, and Hao et al. (2012) considered a free boundary problem which arises in a model of tumour growth with a necrotic core. In this case, Hao et al. proposed a numerical algorithm based on recent developments in numerical algebraic geometry (see references therein). Next, we shall comment, as an example, some papers devoted to the numerical resolution of models quite similar to the one presented here and introduced in Tello (2013). In Barrea & Turner (2005), it is solved a model proposed in Cui & Friedman (2003) by using a spectral numerical method. A mathematical model proposed in Friedman et al. (2006) to study the efficacy of tumour virotherapy is numerically solved in Wang & Tian (2008) through a finite difference second-order accuracy both, space and time. It is the aim of this paper to study the behaviour of the solutions to the model proposed in Tello (2013), and as a consequence, to obtain some inferences regarding the evolution of the different subpopulations of cells within the tumour. In order to perform our numerical simulations, we consider data and parameters, representing the rates at which the reactions occur, previously considered by Molina-Peña & Alvarez (2012), which lead to feasible and consistent scenarios. In fact, in Molina-Peña & Alvarez (2012) it is proposed and studied a mathematical model consisting of a system of ordinary differential equations with an analytical solution which is capable of describing the basic kinetic features of tumour growth during the exponential phase and that is consistent with the CSC hypothesis. The model numerically solved in this paper consists of four hyperbolic equations of first order to describe the evolution of four different subpopulations of cells: CSC, differentiated cells, cells in an intermediate stage between CSC and differentiated cells called progenitor cells and finally dead cells. A fifth equation is introduced to model the evolution of the moving boundary. The system includes non-local terms of integral type in two of the parameters. In Tello (2013), under some restrictions in the parameters it is shown that there exists a unique homogeneous steady state which is stable. The model is considered for the early stage of

NUMERICAL RESOLUTION OF A MODEL OF TUMOUR GROWTH

3 of 29


the cancer when the tumour size is small and necrosis is not present. Experiments show that the growth of the tumour at this stage follows an exponential growth. The scheme of resolution we use in this paper to numerically solve the model presented in Tello (2013) is based on a combination of the characteristics method with a reformulation of the convective term in terms of the total derivative, and the use of a finite element method. Up to our knowledge, the technique employed here has not been used before to deal with the numerical resolution of models of tumour growth similar to the one solved here. As a result of our simulations, we can derive some conclusions regarding the evolution of different subpopulations of cells within the tumour and we can also suggest some therapeutic strategies against cancer. To precise, we obtain that for the early stage of the tumour, our results are in complete agreement with the ones presented in Molina-Peña & Alvarez (2012), that is to say that the CSC subpopulation remains almost invariant and small when one assumes that CSC represent initially less than 1% of the tumour cells. We obtain that the tumour follows an exponential growth as expected. From our numerical results for early times, we can derive that a possible effective measure to detain tumour progression should combine the targeting of CSC with the targeting of progenitor cells, as was also inferred in Molina-Peña & Alvarez (2012). In the present work, we include also the study of the behaviour of the solutions to our PDEs system for longer times in order to obtain some insight about a possible stabilization of the fractions of the different subpopulations of cells. The study of the stabilization of the CSC population was previously addressed in Beretta et al. (2012), where a system of ordinary differential equations or a system of delay differential equations, depending on the fact of not considering or considering an underlying field, were studied in order to describe the dynamics of CSC population in culture. They obtained that the percentage of CSC population is maintained at the same level during a certain time, later, after a rapid decrease, the percentage stabilizes to a small typical value. Our numerical results suggest the stabilization of the densities of the different subpopulations of cells to certain constant values which depend on the set of parameters. In fact, the results show the existence of a very rapid adjustment to a pseudo-equilibrium stage, where the proportion of CSC remains small, which is followed by a slower change to a final equilibrium stage. Maybe, the most remarkable fact that we observe in our results is that the fraction of CSC may vary significantly during the progression of the tumour, becoming on some occasions the predominant type of cells within the tumour. Note that according to recent studies (see, e.g. Johnston et al., 2010 and references therein), CSC subpopulation may comprise any possible proportion of the tumour. Maybe this is a difference between CSC and normal stem cells. Unlike normal stem cells that represent a very small population of cells in the tissue they reside in, CSC can make a rather significant portion of tumour cells (see Shipitsin & Polyak, 2008). The numerical results that we obtain evidence that for longer times, CSC can constitute, depending on the parameter values present in the model, the majority of the tumour cells. Such variations in the CSC fraction are not only consistent with the CSC hypothesis but lend support to it as one of the expected byproduct of the dynamical interactions that are predicted to take place among a relatively small CSC population, its non-stem CC counterpart and the host compartment over time (see Enderling et al., 2013). Recent research papers evidence the tumour growth paradox, which describes the effect that a larger death rate of CC can liberate dormant CSC who in turn produce a larger tumour (see Enderling et al., 2009; Hillen et al., 2013; Borsi et al., 2014). The article is organized as follows. In Section 2, we present the mathematical model which consists of a system of PDEs and a moving boundary. Section 3 is devoted to describe the numerical scheme of resolution. In Section 4, some numerical simulations are presented. Finally, in view of the numerical results some conclusions are drawn.

4 of 29

A. I. MUÑOZ ET AL.

2. The model

CSC → 2CSC and symmetric division to give rise to progenitor transient cells P, CSC → 2P. We do not consider asymmetric division of CSC as in Molina-Peña & Alvarez (2012), where CSC are allowed to undergo asymmetric division, so that one CSC can give rise to another CSC and a P cell. Experimental determination of the relative probability of occurrence of symmetric and asymmetric CSC division is difficult. The regulation of the ratio between symmetric and asymmetric CSC has been currently treated in literature (see, e.g. Morrison & Kimble, 2006; Caussinus & Hirth, 2007; Dingly et al., 2007; Giebel & Bruns, 2008; Turner et al., 2009). It seems that symmetric division rates are the key in the evolution of the tumour and that therapies for effective treatment of cancers should have as a goal to reduce or eliminate symmetric CSC division so that to inhibit their proliferation. In the model studied here, it is considered that the rate of growth of CSC is a non-local function kss which depends on the concentration of CSC in the following way: 1 k1 (t, x, y)s(y, t) dy, kss (x, t) = k0 − |Ω(t)| Ω(t) where k1 is a positive function which measures the influence of CSC in the tumour. The consideration of a non-local integral term in the coefficients is in accordance with the fact that immediate surrounding of a cell plays a determinant role in the division or replication process. We shall assume that k1 s ks (x, t) = k0 − s(y, t) dy, (2.1) |Ω(t)| Ω(t) where k0 and k1 are positive constants. We consider this definition in order to follow exactly at this point the model proposed in Tello (2013), in which k1 is assumed to be constant for simplicity reasons. Other definitions of k1 , will be subject of study in future works. For example, k1 can be considered to be a function depending on the distance between cells, k1 (|x − y|), which would allow us to take into


In this section, we shall describe the mathematical model proposed in Tello (2013). We refer the reader to Tello (2013) and references therein for more details. The four cell subpopulations densities which are considered in the model are ‘s’ CSC density, ‘p’ transit progenitor cells density,‘m’ differentiated cells → density and ‘d’ dead cells density. The variable ‘− v ’ will stand for the velocity of tumour cells, ‘Ω(t)’ for the interior of the tumour and finally ‘∂Ω(t)’ for the boundary of the tumour. We shall assume that the different types of cells are physically identical with a continuous distribution into the tumour and that cells interact through the exchange of molecules and may evolve spontaneously from one state to another. The model is considered for the early stage of the cancer, so there are no nutrient limitations and hence, angiogenesis occurs at a rate that ensures the accessibility of nutrients to sustain constant growth, necrosis is not present and the tumour follows an exponential growth. The CSC division process is regulated by a chemical feedback with the cell’s neighbourhood which determines the type of division (symmetric or asymmetric). In our case, we only consider two process regarding the evolution of CSC, which are symmetric division to give rise to new CSC, that is to say


5 of 29

consideration the immediate surrounding of a cell. Anyway, note that the measure of the concentration of CSC in the domain is present when one considers definition (2.1), due to the presence of the non-integral term. CSC may produce progenitor cells at ratio ksp in a similar way that in the previous growth rate kss . We shall consider k3 p s(y, t) dy, (2.2) ks (x, t) = k2 + |Ω(t)| Ω(t)

∂s → + div(− v s) = kss s − ksp s, 0 < t < T, x ∈ Ω(t), ∂t ∂p → + div(− v p) = ksp s + kpp p − kpm p − kp p, 0 < t < T, x ∈ Ω(t), ∂t ∂m → + div(− v m) = kpm p − km m, 0 < t < T, x ∈ Ω(t), ∂t ∂d → + div(− v d) = kp p + km m − kd d, 0 < t < T, x ∈ Ω(t), ∂t

(2.3)

where the coefficients kss , ksp , kpm , kp , km and kd have been described before. The system (2.3) is completed with appropriate initial data s(0, x) = s0 (x),

p(0, x) = p0 (x),

m(0, x) = m0 (x),

d(0, x) = d0 (x), x ∈ Ω0 .

(2.4)

The conservation of the mass laws for the densities of the cells, assumed homogeneous tumour density, gives s + p + m + d = constant, (2.5) where the constant is assumed to be 1. Adding (2.3) and thanks to (2.5), we obtain that the balance of mass is then given by → (2.6) div(− v ) = kss s + kpp p − kd d. Note that being the tumour tissue treated as a porous medium, we can assume that the velocity of the fluid is described by Darcy’s law − → v = −∇σ , where σ is the pressure, and then obtain that, −Δσ = kss s + kpp p − kd d,


where k2 and k3 are assumed to be positive constants. On the other side, P-cells can either self-renew or they can differentiate into differentiated M -cells at constant rates kpp and kpm , respectively. We assume that M -cells have a neglected capacity to proliferate, and therefore the corresponding rate growth factor does not appear. Finally, we consider that subpopulation tumour cells P die at constant rates kp and D cells dead cells) decompose at rate kd . The death rate of CSC, ks is assumed null. In Section 4, we shall consider the case where ks is allowed to be positive by slightly modifying one of the equations of the model, so that we can analyse its role in the evolution of the subpopulation densities. Due to the proliferation and removal of cells, there is a continuous motion within the tumour. We → shall represent this movement by a velocity field − v (see Friedman, 2006). We can then write the mass conservation laws for the densities s, p, m and d, by means of the following system of first-order hyperbolic equations (see Tello, 2013):

6 of 29

A. I. MUÑOZ ET AL.

n

direction. We assume radially symmetric distribution of cells and spherical tumours, i.e. Ω(t) := {x ∈ R3 , such that |x| R(t)}, where R(t) denotes the radius of the tumour, and we consider therefore, the radial component of the velocity v. Note that for in vitro experiments, multicellular tumours have a nearly spherical shape. As previously mentioned, we assume that the velocity of the free boundary is equal to the velocity of the fluid flow at the boundary for continuity (see Friedman, 2006; Tello, 2013): R˙ =

dR = v(R(t), t) dt

for t > 0.

(2.7)

Now, let us introduce the spacial variable r ∈ I := (0, 1) such that r :=

|x| , R(t)

and define v(r, t) := v(x, t). Since → → → ∇ · (− v s) = s∇ · − v +− v · ∇s = s(kss s + kpp p − kd d) +

v ∂s R ∂r

the system (2.3–2.7) becomes (see Tello, 2013): R˙ v ∂s ∂s + −r + = kss s − ksp s + s(−kss s − kpp p + kd d), ∂t R R ∂r R˙ v ∂p ∂p + −r + = ksp s + p(kpp − kpm − kp − kss s − kpp p + kd d), ∂t R R ∂r ∂m R˙ v ∂m + −r + = kpm p − km m + m(−kss s − kpp p + kd d), ∂t R R ∂r R˙ v ∂d ∂d + −r + = kp p + km m − kd d + d(−kss s − kpp p + kd d), ∂t R R ∂r ∂ 2v r = r2 (kss s + kpp p − kd d), ∂r R

(2.8) (2.9) (2.10) (2.11) (2.12)

with boundary conditions v(1, t) =

dR dt

and

v(0, t) = 0

(2.13)


→ so, we could pose our problem in terms of the pressure σ instead of the velocity field − v , but we shall → follow the formulation considered in Tello (2013), which is written in terms of the velocity field − v. Regarding the boundary conditions, it is assumed that the pressure σ on the surface of the tumour is equal to the surface tension (see Friedman, 2004). The motion of the free boundary is given by the continuity equation ∂σ − → → = Vn , v ·− n =− ∂n → where − n is the outward normal and V is the velocity of the free boundary in the outward normal

7 of 29

NUMERICAL RESOLUTION OF A MODEL OF TUMOUR GROWTH −3

6

Results for p

Results for s

x 10

0.1812 time=60 time=20 time=10 time=5

5.5

time=60 time=20 time=10 time=5

0.1811 0.181

5 0.1809 0.1808

4

p

s

4.5

0.1807

0.1805 3 0.1804 2.5 2

0.1803 0

0.2

0.4

0.6

0.8

0.1802

1

0

0.2

0.4

0.6

0.8

1

r Results for d

Results for m 0.394

0.424 time=60 time=20 time=10 time=5

0.3935

time=60 time=20 time=10 time=5

0.4238 0.4236 0.4234

0.393 d

m

0.4232

0.3925

0.423 0.4228 0.4226

0.392

0.4224 0.4222

0.3915

0

0.2

0.4

0.6

0.8

1

0.422

0

0.2

0.4

0.6

0.8

1

r

r

Fig. 1. Results obtained with the set of parameter values suggested in Molina-Peña & Alvarez (2012), which are: k0 = 1, k1 = 0.01, p k2 = 0.01, k3 = 0.01, kp = 5.35, kpm = 4.28, kp = 0.1, km = 1, kd = 0, ks = 0. We present the values for s (top left), p (top right), m (bottom left) and d (bottom right), for different times.

and initial data s(r, 0) = s0 ,

p(r, 0) = p0 ,

m(r, 0) = m0 ,

d(r, 0) = d0

and

R(0) = R0 .

(2.14)

It is assumed that initial values are regular functions and that they satisfy 0 < s0 < 1,

0 < p0 < 1,

0 < m0 < 1

and

0 < d0 < 1.

We shall see in Section 4 that in all the simulations presented in Figs 1–8, where homogeneous initial data satisfying s0 + p0 + m0 + d0 = 1 are considered, the densities s, p, m and d satisfy 0 < s < 1,

0 < p < 1,

0 < m < 1,

0 < d < 1 and

s+p+d +m=1

for t > 0. For initial data not satisfying this constraint, we have that s + p + m + d tends to 1 as t increases.


0.1806

3.5

8 of 29

A. I. MUÑOZ ET AL. 7

9

Results for v

x 10

Results for v 250 time=60 time=20 time=10 time=5

8 7

time=20 time=10 time=5 200

6 150 v

v

5 4

2

50

1 0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

r

0.6

0.8

1

r 8

Result for the radius 700

3

600

Result for the radius

x 10

2.5

500

radius

radius

2 400

1.5

300 1 200 0.5

100

0

0

5 −3

1

10 t

15

20

30 t

40

50

60

50

60

Result for the total density

1.8 1.6 1.4 1.2

total density

0.7 0.6 0.5 0.4

1 0.8

0.3

0.6

0.2

0.4

0.1

0.2 0

10

2 time=60 time=20 time=10 time=5

0.9

0

0

Results for u

x 10

0.8

u

0

20

0.2

0.4

0.6 r

0.8

1

0

0

10

20

30 t

40

Fig. 2. We present the results obtained for v at times t = 5, 10, 20, 60 (top left) and for t = 5, 10, 20 (top right). The results for the radius at t = 20 (middle left) and t = 60 (middle right) are also shown. Finally, we present the results for u at times t = 5, 10, 20, 60 (bottom left) and the evolution for the total density (bottom right). These results have been obtained with the same parameter values that in Fig. 1.


100 3

9 of 29

NUMERICAL RESOLUTION OF A MODEL OF TUMOUR GROWTH Results for s

Results for p

0.2

0.2

0.18 0.19

0.16 0.14

0.18

p

s

0.12 0.1

0.17

0.16

0.04

0.15

0.02 0

0

100

200

300 t

400

500

600

0

100

200

Results for m

300 t

400

500

600

400

500

600

Results for d

0.75

0.45

0.7

0.4

0.65

0.35

0.6 0.3 d

m

0.55 0.25

0.5 0.2 0.45 0.15

0.4

0.1

0.35

0

100

200

300 t

400

500

600

0.05

0

100

200

300 t

Fig. 3. We present the results obtained, with the same parameter values as in Fig. 1, for the variables s (top left), p (top right), m (bottom left) and d (bottom right) at the point r = 1 in order to show the behaviour of the densities for longer times (t ∈ [0, 600]). These results suggest the existence of different stages in the evolution of the densities. Note that at early times, there exists a very rapid adjustment to a pseudo-equilibrium state, which is followed by a slower change to a final equilibrium state.

3. Numerical resolution In this section, we shall describe the scheme of resolution we employ to solve numerically the model consisting of (2.8–2.14). The numerical resolution of the model and subsequent numerical simulations allow us to analyse the effect of variations in the relationships between the different cellular events on the fractions of the different subpopulations of cells and also to study the stabilization of the densities of these subpopulations to homogeneous steady states. Due to the hyperbolic nature of Equations (2.8)–(2.11), the numerical method we are going to used follows the one described in Bercovier et al. (1983) (see also the references cited there). This technique has been successfully used in other fields (see, e.g. Bayada et al., 1998; Calvo et al., 1999, 2002). The method is based on the method of characteristics and the reformulation of the convection terms in Equations (2.8–2.11) in terms of the total derivative. Finally, a finite element method is employed. We shall assume that the variables satisfy some regularity requirements that allow us to develop the


0.08 0.06

10 of 29

A. I. MUÑOZ ET AL. −3

6

Results for p with t=20

Results for s with t=20

x 10

0.2 type 1 type 2 type 3 type 4 type 5

0.195

4

0.19

3

0.185

p

s

5

0.18

1

0.175

0

0.2

0.4

0.6

0.8

0.17

1

0

0.2

0.4

Results for m with t=20

1


0.404 0.402

0.402

0.4

0.4

0.398

0.398

0.396

0.396

0.394

0.394

0.392

0.392

0.39

0.39

0.388

0.388 0

0.2

0.4

0.6

0.8

type 1 type 2 type 3 type 4 type 5

0.404

d

m

0.8

Results for d with t=20

0.406

0.386

0.6 r

r

1

0.386

0

0.2

0.4

r

0.6

0.8

1

r −3

Results for v with t=20 1

450 type 1 type 2 type 3 type 4 type 5

400 350

Results for u with t=20

x 10


0.9 0.8 0.7

300

0.6 v

u

250 0.5

200 0.4 150

0.3

100

0.2

50 0

0.1 0

0.2

0.4

0.6 r

0.8

1

0

0

0.2

0.4

0.6

0.8

1

r

Fig. 4. Results obtained with the same values for the parameters as the ones considered in Fig. 1, with the following exceptions depending on the type. In type 1, the parameters are the same than in Fig. 1. In type 2, we have considered kpm = 4.32, instead p p of 4.28. In type 3, we have considered kp = 5.38 instead of 5.35. In type 4, we choose kpm = 4.18 and in type 5, kp = 5.32. We present the results with t = 20 for the variables s (top left), p (top right), m (middle left), d (middle right), v (bottom left) and u (bottom right).


2

0



11 of 29

Results for the radius with t=20 1200 type 1 type 2 type 3 type 4 type 5

1000

600

400

200

0

0

5

10 t

15

20

Fig. 5. We present the results of the tumour radius R(t) corresponding to the parameters considered in Fig. 4 with t = 20.

numerical scheme. We would like to remark that, to the best of our knowledge, the technique we use has not been previously considered regarding the numerical resolution of tumour growth models of the same mathematical nature as the one solved here, being spectral techniques and finite differences schemes the most commonly used (see Barrea & Turner, 2005; Wang & Tian, 2008). Next we describe the steps we follow in order to solve numerically the problem. Note that despite the fact that we are dealing with a free boundary problem, the spatial domain of resolution is fixed due to the change of variable r := |x|/R(t), to precise Ω = (0, 1). An important remark is that no boundary values need to be prescribed for the hyperbolic equations as the characteristics remain inside the domain [0, 1] × (0, T] (see Friedman, 2006). Let us introduce some notation we shall follow through the exposition for the right-hand sides of the respective hyperbolic equations: fs = kss s − ksp s + sh, fp = ksp s + p(kpp − kpm − kp + h), fm = kpm p − km m + mh, fd = kp p + km m − kd d + dh, where h is defined by h(s, p, d) = −kss s − kpp p + kd d. We shall used a time marching scheme to solve (2.8–2.14). Let us consider a fixed time step of discretization Δt. Note that at t = 0, we have initial conditions for the variables s, p, m and d and for the radius R. First we have to compute the values for kss and ksp , where the integral non-local terms are computed using the trapezoidal rule of numerical integration. The values of these two parameters are updated at every step in time considered. Then, after initializing the velocity v using a finite different scheme, we can initiate our scheme of resolution.


radius

800

12 of 29


3


x 10

Results for p with t=20 0.1814 type 1 type 2 type 3 type 4 type 5

0.1813

2

0.1812

1.5

0.1811

p

s

2.5

0.181

0.5

0.1809

0

0.2

0.4

0.6

0.8

0.1808

1

0

0.2

0.4

r

0.6

0.8

1

r



0.394


0.3938


0.4248 0.4246

0.3936

0.3934

d

m

0.4244 0.4242 0.424 0.3932 0.4238 0.393 0.4236 0.3928

0

0.2

0.4

0.6

0.8

0.4234

1

0

0.2

0.4

r

0.6

0.8

1

r Results for the radius with t=20

600 type 1 type 2 type 3 type 4 type 5

500

radius

400

300

200

100

0

0

5

10 t

15

20

Fig. 6. Results obtained with the same values for the parameters as the ones considered in Fig. 1 (taking as a reference and considered in type 1), with the following exceptions depending on the type. In type 2, we consider k2 = 1 (instead of k2 = 0.01). In type 3, we take k0 = 0.01 (instead of k0 = 1). In type 4, we choose k3 = 1 (instead of k3 = 0.01). And finally in type 5, we consider k1 = 10 (instead of k1 = 0.01). We present the results with t = 20 for the variables s (top left), p (top right), m (middle left), d (middle right) and R (bottom centre).


1

0


13 of 29


4

Results for p with t=20


x 10

0.184 type 1 type 2 type 3 type 4

3.8 3.6

type 1 type 2 type 3 type 4

0.183 0.182

3.4 0.181 3.2 s

p

0.18

3

0.178

2.6

0.177

2.4 2.2

0

0.2

0.4

0.6

0.8

0.176

1

0

0.2

0.4

0.6

0.8

1

r

r Results for m with t=20


0.397

0.434 type 1 type 2 type 3 type 4

0.396 0.395


0.432 0.43

0.394 0.428 d

m

0.393 0.426

0.392 0.424 0.391 0.422

0.39

0.42

0.389 0.388

0

0.2

0.4

0.6

0.8

1

0.418

0

0.2

0.4

r

0.6

0.8

1

r

Results for v with t=20

Results for the radius with t=20

250

700 type 1 type 2 type 3 type 4

200


600

500

v

radius

150

400

300

100

200 50 100

0

0

0.2

0.4

0.6 r

0.8

1

0

0

5

10 t

15

20

Fig. 7. Results obtained with the same values for the parameters as the ones considered in Fig. 1 (taking as a reference and considered in type 1), with the following exceptions depending on the type. In type 2, we introduce ks in the model and consider ks = 0.01. In type 3, we consider kp = 0.12 (instead of the original value kp = 0.1) and in type 4, we choose kp = 0.09. We present the results with t = 20 for the variables s (top left), p (top right), m (middle left), d (middle right), v (bottom left) and R (bottom right).


0.179 2.8

14 of 29


2.7158


x 10

Results for p with t=20 0.1808 type 1 type 2 type 3

2.7158

type 1 type 2 type 3

0.1808 0.1808

2.7158

0.1808 0.1808

s

p

2.7158

2.7158

0.1808

0.1808 0.1808

2.7158 0.1808 2.7158

0

0.2

0.4

0.6

0.8

1

0.1808

0

0.2

0.4

r

0.6

0.8

1

r



0.8

0.7 type 1 type 2 type 3

0.7


0.6

0.5

0.5

0.4 d

m

0.6

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

r

0.6

0.8

1

r

Results for v with t=20

Results for the radius with t=20 600

250 type 1 type 2 type 3


500

200

400

v

radius

150 300

100 200 50

0

100

0

0.2

0.4

0.6 r

0.8

1

0

0

5

10 t

15

20

Fig. 8. Results obtained with the same values for the parameters as the ones considered in Fig. 1 (taking as a reference and considered in type 1), with the following exceptions depending on the type. In type 2, we consider the extreme value km = 0 (instead of km = 1). In type 3, we consider km = 3. We present the results with t = 20 for the variables s (top left), p (top right), m (middle left), d (middle right), v (bottom left) and R (bottom right).


0.1808

2.7158


15 of 29

We shall first describe the method of resolution employed for solving Equation (2.8) relative to the unknown s. Note that due to the structure of Equations (2.9–2.11), the numerical schemes of resolution for them are almost identical to the one for (2.8) except for the different right-hand sides of the equations. We shall follow the notation and ideas presented in Bercovier et al. (1983). As before said, the scheme used here is based on the characteristics method and an approximation of the total derivative that appears when the convection term is written in conservative form. Therefore, let us consider the total derivative Ds/Dt defined as follows:

where, in this case, u would be an artificial velocity field given by

v R˙ u = −r + R R

.

Therefore, Ds = fs + s∇ · u. Dt

(3.1)

Note that when u is actually the velocity of the fluid and the fluid is considered to be incompressible, then ∇ · u is null. In this problem, u is an artificial velocity and ∇ · u is not necessary null so we will have to take it into consideration and compute it. To do that, as we are dealing with a spherically symmetric problem, 1 ∂ 2 1 ∂ R˙ 1 ∂ r2 v 2 r u = 2 ∇ ·u= 2 r −r + 2 , r ∂r r ∂r R r ∂r R hence, using the fact that

∂ 2v r = −r2 h, ∂r R

(3.2)

we get the identity ∇ · u = −3

R˙ + h. R

In our particular problem, we shall consider the total derivative of Js, given by DJs ∂ (r, t) = [J(r, t; τ )s(X (r, t; τ ), τ )]τ =t , Dt ∂τ where J is the Jacobian associated with the change of coordinates from Eulerian to Lagrangian coordinates. The Jacobian J(r, t; τ ) of the change of coordinates is defined as follows: J(r, t; τ ) := 1 −

τ

t

∇ · u(X (r, t; τ ), τ ) dτ .

(3.3)

Note that the presence of J arise from the application of the characteristics method when one considers the conservative form for the convection term. We shall denote the characteristics X by X (r, t; τ ), as r


Ds ∂s ∂s = + u + s∇ · u, Dt ∂t ∂r

16 of 29

A. I. MUÑOZ ET AL.

and t are parameters for the solution X of ⎧ ⎨ dX (r, t; τ ) = u(X (r, t; τ ), τ ), dτ ⎩X (r, t; t) = r,

(3.4)

∂J ∂s ∂X DJs ∂s ∂s = s+J + + ∇(su), = Dt ∂τ ∂r ∂τ ∂τ ∂t τ =t and the problem we have to solve numerically for the variable s turns out to be: DJs = fs + s∇ · u, 0 < t < T, r ∈ [0, 1], Dt s(0, r) = s0 (r), r ∈ [0, 1].

(3.5) (3.6)

As before said, it is no necessary to prescribe boundary conditions because in our problem X (r) ∈ [0, 1] as was proved in Friedman (2006). The variational formulation corresponding to the problem (3.5–3.6) is the following: Ω

DJs − fs − s∇ · u w dr = 0 ∀w ∈ L2 (0, 1), ∀t ∈ (0, T). Dt

(3.7)

To discretize in time, we use the following formula: DJs (r, tn+1 ) [J(r, tn+1 ; tn+1 )s(X (r, tn+1 ; tn+1 ), tn+1 ) − J(r, tn+1 ; tn )s(X (r, tn+1 ; tn ), tn )]/Δt. Dt

(3.8)

Note that if we consider t = tn+1 in (3.3–3.4), we get that J n+1 = J(r, tn+1 ; tn+1 ) = 1,

X (r, tn+1 ; tn+1 ) = r,

(3.9)

therefore (3.8) becomes DJs (r, tn+1 ) [sn+1 − J n (r)sn (X n (r))]/Δt, Dt

(3.10)

where sn (r) = s(r, tn ), J n (r) = J(r, tn+1 ; tn ) and X n (r) = X (r, tn+1 ; tn ) is the solution of (3.4) that will be in r at time tn+1 . The terms fs and s∇ · u in (3.7) are evaluated in the previous time step, therefore they will be known, in other words, we shall consider an explicit scheme. Coming back to the variational formulation (3.7), we have that

s Ω

n+1

w dr =

J (r)s (X (r))w dr + Δt n

Ω

n

n

Ω

(fsn + sn (∇ · u)n )w dr

∀w ∈ L2 (0, 1)

(3.11)


which means that X (r, t; ·) should be the particle path that passes at r at time t. The system (3.4) is backward in time but when u is regular enough, say Lipschitz continuous, it has a unique solution. Then, we obtain that


17 of 29

being fsn (r) = kss sn (r) − ksp sn (r) + sn (r)hn (r) = kss sn (r) − ksp sn (r) + sn (r)(−kss sn (r) − kpp pn (r) + kd d n (r)), (∇ · u)n = −3

R˙ n + hn (r), Rn

∂ ∂r

n+1 (ri ) 2v ≈ ri Rn vn+1 (ri ) =

2 vn+1 (ri−1 ) ri2 vn+1 (ri ) ri−1 − Rn Rn

1 2 hn+1 (ri−1 )), Δr = − (ri2 hn+1 (ri ) + ri−1 2

2 ri−1 1 1 r2 vn+1 (ri−1 ) − ΔrRn hn+1 (ri ) − i−1 ΔrRn hn+1 (ri−1 ) 2 2 2 ri2 ri

and hence,

n+1

v

Δr (ri ) = 1 − ri

2 n+1

v

1 1 (ri−1 ) − ΔrRn hn+1 (ri ) − 2 2

Δr 1− ri

2 ΔrRn hn+1 (ri−1 ),

(3.12)

which is the scheme we shall consider in order to update the values of the velocity for the simulations presented in Figs 1–9. Instead of considering the scheme (3.12) to solve the problem for the velocity v, we could consider the following implicit midpoint approximation: ∂ ∂r

ri2

vn+1 (ri ) Rn

2 hn+1 (ri−1/2 ), ≈ −ri−1/2

which can be performed in the following way:

n+1

v

Δr (ri ) = 1 − ri

2 n+1

v

Δr (ri−1 ) − 1 − 2ri

2

1 ΔrRn (hn+1 (ri−1 ) + hn+1 (ri )). 2

(3.13)


where Rn and R˙ n make reference to the values of the radius and its derivative in the instant tn = nΔt, and finally, sn (r) = s(r, tn ), pn (r) = p(r, tn ) and d n (r) = d(r, tn ). In order to discretize the problem (3.7) with respect to the coordinate r, we use Lagrange linear finite elements. After applying some quadrature formula for integration and some algebraic operations, = sn+1 (ri ), and i = 1, 2, . . . N, being {ri } we obtain a linear system of equations for the unknowns sn+1 i the partition of the domain [0, 1] considering Δr as the step for the space discretization. The matrix of coefficients associated to the system is a band tri-diagonal one. Then system is solved by using the Thomas algorithm. For the rest of the variables p, m and d, the scheme of resolution is the same as for s, with the exception of considering the corresponding r.h.s., fp , fm and fd . After solving the problems for the variables s, p, m and d, we pass to compute the values for v. The equation for v is solved using a trapezoidal approximation, which is unconditionally stable and of high accuracy. Note that v(0) = 0 for all times. The rest of the nodal values in an arbitrary (n + 1)th step in time are computed as follows. Taking into account Equation (3.2), we consider the following approximation:

18 of 29

A. I. MUÑOZ ET AL.

−3

3

Results for p

Results for s

x 10

0.3

2.8 2.6

0.25

2.4 0.2

2

p

s

2.2

0.15

1.6 1.4

1

0.1

time=0 time=5.000000e−002 time=1 time=20

1.2 0

0.2

0.4

0.6

0.8

0.05

1

time=0 time=5.000000e−002 time=1 time=20 0

0.2

0.4

r

0.6

0.8

1

r

Result for the radius

Mean total densities

800

1.12

700 1.1 600 1.08 mean density

radius

500 400 300

1.06

1.04 200 1.02 100 0

0

5

10 t

15

1

20

0

5

10 t

15

20

Results for u

Result for v with t=0.05 0.005

0.6

0

0.5

−0.005

0.4

−0.01

v

u

0.7

0.3

−0.015

0.2

−0.02

0.1

−0.025

0

0

0.2

0.4

0.6 r

0.8

1

−0.03

time=5.000000e−002 time=1 time=20 0

0.2

0.4

0.6

0.8

1

r

Fig. 9. We present the values obtained considering the same set of parameters as in Fig. 1, but with the following initial conditions: s0 (ri ) = 0.002 + 0.001 cos(π i/10), p0 (ri ) = 0.198 + 0.1 cos(π i/10), m0 = 0.8 and d0 = 0 for i = 1, . . . , N, at different times.


1.8

19 of 29


4.08

Steady states for s

Results for s

x 10

0.1968

Trapezoidal approx., dr=0.01, t=40 Implicit midpoint, dr=0.001, t=40 Implicit midpoint, dr=0.01, t=40

4.07

type1 type2 type3

0.1968 0.1968

4.06 0.1968

4.05 s

s

0.1968 0.1968 0.1968

4.03 0.1968

4.02

4.01

0.1968 0.1968

0

0.2

0.4

0.6

0.8

r

1

0

0.2

0.4

0.6

0.8

1

r

Fig. 10. On the left, we present the values obtained for s at t = 40 considering the same set of parameters and initial data as in Fig. 9, with the trapezoidal approximation for Δr = 10−2 , and with the implicit midpoint method for Δr = 10−3 and Δr = 10−2 . On the right, we show the steady states for the three following different cases. The steady state s = 0.1968492188072 in type 1 is obtained considering the trapezoidal approximation with Δr = 10−2 and initial data as in Fig. 9. In type 2, the steady s = 0.1968492188685 is obtained with the implicit midpoint approximation, Δr = 10−3 and initial data as in Fig. 1. And in type 3, the steady state, which takes values between sl = 0.1968492188665 and su = 0.1968492189600, is obtained with the implicit midpoint approximation, Δr = 10−2 and initial data as in Fig. 9.

The scheme (3.13), which in principle is more accurate for approximating the r.h.s. term in the equation for the velocity, −r2 h(r), proves to be sensitive to oscillations. Hence, in order to avoid undesirable oscillations in the solutions, we have to consider small enough steps Δr. Due to this fact, the scheme (3.13), in this case, results to be more expensive than the one given in (3.12), being this cost, from our personal point of view, not justified as no significant differences in the results are obtained. In Section 4, we shall illustrate and complete the above comments (see Table 2 and Fig. 10). Finally, we can find that R˙ n+1 = vn+1 (rN = 1), and the radius can be updated in the following way: n+1 R = Rn + vn+1 (rN )Δt. The values of the parameters kss and ksp are updated too, using the updated values for s. 4. Numerical simulations and discussion This section is devoted to describe the different numerical results regarding the evolution of the four subpopulations of cells that we have obtained for different sets of parameters, both for the early stage and for later stages of tumour progression. Experimentally studying the balance between the different cellular events involved in the process of tumour growth is not a trivial matter. In the model studied here, the rates at which the reactions, such us cell division, differentiation or death, occur, are represented by the parameters. In order to obtain feasible solutions, that is to say, solutions rendering dynamical behaviours consistent with experimental observations of the evolution of the tumour, Molina-Peña & Alvarez (2012) considered the following constraints. They consider that the constraint ds/dt ≈ 0 at early stages is of particular importance (see Collins et al., 2005; Maitland & Collins, 2008), as it is clear evidence of the crucial role of the CSC reservoir in the tumour. Other constraints at this stage consist in assuming that s < 0.01, p ≈ 0.2, and that differentiated M cells should constitute the majority of tumour cells, i.e. m ≈ 0.8. It is also assumed, based on biological knowledge, that the intrinsic apoptosis rate increases as cell becomes


4.04

20 of 29

A. I. MUÑOZ ET AL.

more differentiated. Therefore, ks < kp < km .

s0 = 0.002,

p0 = 0.198,

m = 0.8 and

d = 0.

In Fig. 1, we present the results we have obtained for the densities s, p, m and d for different times with the following set of parameters: k0 = 1, k1 = 0.01, k2 = 0.01, k3 = 0.01, kpp = 5.35, kpm = 4.28, kp = 0.1, km = 1, kd = 0 and ks = 0. These values are in correspondence with the ones considered in Molina-Peña & Alvarez (2012) with the following exceptions: First of all, note that the kj notation in our model does not follow exactly the same as in Molina-Peña & Alvarez (2012). In the model presented here ks , the rate for CSC death, does not appear so it is assumed to be zero, while in Molina-Peña & Alvarez (2012) it is assumed to be small but positive. Anyway, we shall introduce a slight modification in our model which will allow us to consider ks and analyse the results for small positive values as in Molina-Peña & Alvarez (2012). Finally, in Molina-Peña & Alvarez (2012) the constants related to symmetric division of CSC are k1 and k3 , while in the model studied here are kss and ksp . The coefficients kss and ksp in our model (see (2.1) and (2.2)) depend on the values of s through the integral no local term (|Ω(t)|)−1 Ω(t) s(y, t) dy pondered by the constants k1 and k3 , respectively (which are not the same k1 and k3 than in MolinaPeña & Alvarez, 2012). This integral term is numerically approximated with the trapezoidal rule. We do not consider the asymmetric division of CSC, so the constant k2 in Molina-Peña & Alvarez (2012) has no a equivalent rate in our model. The steps in space and time used in all the simulations presented in Figs 1–9 are Δr = 10−2 and Δt = 10−4 , and as we mentioned before, the scheme used to solve the equation for the velocity is the trapezoidal approximation given in (3.12). It can be observed in Fig. 1 that for early times, we obtain feasible solutions according to MolinaPeña & Alvarez (2012), as s < 0.01 and ds/dt ≈ 0. and the results are in correspondence with the ones


In our simulations, we shall consider as a reference the values of the parameters and data used in Molina-Peña & Alvarez (2012). Based on an exploration of the parameter space, we shall derive some conclusions. Regarding the evolution of the different subpopulations of cells at early times, we obtain that our results are in agreement with the ones presented in Molina-Peña & Alvarez (2012), as we shall see later. We also include the analysis of the behaviour of solutions for longer times. We shall show some simulations where it can be observed that after some time the constraint ds/dt ≈ 0 is not longer satisfied as it varies significantly, and that for longer times, the fractions of the different subpopulations of cells stabilize to certain values (see Figs 3 and 10 (right), and Tables 5 and 6). It is important to remark that CSC can comprise different proportions of the tumour, becoming on some occasions the predominant type of cells within the tumour (see Table 5). These results may be admissible according to recent research works (see Shipitsin & Polyak, 2008; Johnston et al., 2010) and may be not only consistent with the CSC hypothesis but lend support to it (see Enderling et al., 2013). According to Molina-Peña & Alvarez (2012), we shall consider as initial conditions for our simulations the following initial constant values for s, p, m and d:


21 of 29

presented in Molina-Peña & Alvarez (2012). The solutions are constants for each t fixed, in other words, for t fixed, they do not vary with respect to r, and this fact has been obtained in all the simulations with uniform initial data. In Fig. 2, we present the results obtained for the velocity v, note that it presents a linear profile for each t which is expected due to the equation for v. We justify this behaviour as follows. Note that the equation for v is the one given in (3.2), i.e.

and that for a fixed time, h is constant (in fact a negative constant) due to the fact that all the variables do not depend on r. R is also constant, so we find that a particular solution to the linear first order ODE is the following: −Rh r, (4.1) v= 3 satisfying that v(0) = 0 for all times. In our simulations, the results obtained for the slopes of the velocity using the numerical values for the velocity at the nodes are in agreement with the value obtained for the slope, −Rh/3, at times for which solutions are homogeneous, in other words, they do not depend on r. The results presented in Table 2, illustrate, in some way, this fact for the case considered in Fig. 1, and show the well performance of the scheme (3.12). In Fig. 2, the results obtained for R(t) are shown for t = 20 and t = 60. Note that it follows an exponential growth as expected. Regarding the values for the artificial velocity u, we obtained that the profile and values were essentially the same in all the simulations presented in Figs 1–8. Note that its values are of the order of 10−3 . With regards to the conservation of the densities (2.5), we see that in all the numerical simulations performed with initial conditions not depending on r and satisfying that s0 + p0 + m0 + d0 = 1 then it is satisfied that s + p + m + d = 1, for all times. In Fig. 2, we present the results obtained for the evolution of the total densities where it can be checked that it is satisfied that s + p + m + d = 1. In Fig. 9, we shall present the results obtained with the same set of parameters than in Fig. 1 but with initial data | 1. These results would illustrate the fact that when initial depending on r, and with s0 + p0 + m0 + d0 = data do not satisfy the constraint s0 + p0 + m0 + d0 = 1, then s + p + m + d tends to 1 as t increases. More comments about this point can be found at the end of this section. As mentioned before, the steps in space and time used in all the simulations presented in Figs 1–9, are Δr = 10−2 and Δt = 10−4 . In Table 1, we present the values for the densities s and d, and for the radius R at t = 20, obtained with the parameter values used in Fig. 1 for different steps of discretization both in space and time, in order to illustrate the almost null dependence of the results on the steps of discretization. This fact is due to the small values of the artificial velocity u which are of the order of 10−3 . Note that it is the error in time discretization that dominates the error in space discretization. In Table 2, we present the results obtained at t = 20 for the different variables, with the parameters and initial data considered in Fig. 1, using the different schemes (3.12) and (3.13) for solving the equation for the velocity. We show the values obtained with the implicit midpoint method (3.13) for different steps in space and also the ones obtained with the trapezoidal approximation (3.12). It can be observed that there are no significant differences in the results when one considers homogeneous initial data. Regarding the monotonicity of the densities with respect to time, we present, in Fig. 3, the results obtained for longer times considering the same the set of parameters as in Fig. 1, in order to illustrate the


∂ 2v r = −r2 h, ∂r R

22 of 29

A. I. MUÑOZ ET AL.

Table 1 Values for the densities s and p, and for the radius R at t = 20 with the parameter values considered in Fig. 1 for different steps of discretization both in space and time Δt

Δr

s

p

R

10−2 10−2 5 × 10−3 10−3 10−2

2.715795115974555E−03 2.715790004047166E−03 2.715808470467057E−03 2.715823199684347E−03 2.715738898661894E−03

0.180838535276993 0.180838536160834 0.180838532506344 0.180838529591444 0.180838544996741

663.769678859659 663.138984447575 662.977470233778 662.925792530568 656.878245830995

Table 2 Results obtained at t = 20 for Δt = 10−4 , with the trapezoidal approximation described in (3.12) and with the implicit midpoint scheme (3.13) for different steps in space.

s p m d Radius R v(r = 0.5) v(r = 1) Slope = −Rh/3

Trapezoidal approx., Δr = 10−2 0.00271579511597 0.18083853527699 0.39285805069938 0.42358761890738 663.769678859659 107.349759930725 214.667321373168 214.663530898392

Impl. mid., Δr = 10−2 0.00271579511598 0.18083853527702 0.39285805069948 0.42358761890762 663.446294162300 107.265277215653 214.546645832029 214.558948084275

Impl. mid., Δr = 10−3 0.00271582831245 0.18083852870769 0.39285803715382 0.42358760582698 663.552988180224 107.293149229959 214.586459399802 214.593452786106

evolution of the fractions of the different subpopulations of cells. Note that for a certain time, solutions are no longer feasible, according to Molina-Peña & Alvarez (2012), as ds/dt ≈ 0 is not satisfied and s > 0.01. The solutions stabilize to certain values and the densities add up to one for all the times. At the end of the section, we shall present more results regarding the long time evolution of solutions (see Tables 4 and 5, and Fig. 10 (right)). As before mentioned, this results might be admissible according to recent research regarding the proportion of the CSC population in solid tumours (see, e.g. Shipitsin & Polyak, 2008; Johnston et al., 2010; Enderling et al., 2013 and references therein). In order to illustrate the effect of varying the parameters kpm and kpp , we present some numerical results for t = 20 in Figs 4 and 5. The type 1 corresponds to the values obtained with the same set of parameters considered in Fig. 1, that we take as a reference. In type 2, we have considered kpm = 4.32, instead of the original value kpm = 4.28. In type 3, we have considered kpp = 5.38 instead of the original value kpp = 5.35. Finally, in type 4, we choose kpm = 4.18 and in type 5, kpp = 5.32. The numerical results show that when we increase the value of kpm , s and d increase, and p and m (in a smaller rate than p) decrease. The effect of these fact is a reduction in the radius of the tumour. This same effect can be also obtained reducing the value of kpp . However, when one decreases the value of kpm or increases the value of kpp , we get smaller values for s and d, bigger values for p and m, and the radius R, is considerably bigger than the one obtained in the reference case (type 1). These results are in agreement with the suggestion that the progenitor population is also of importance in order to detain the tumour progression (see Kern & Shibata, 2007), and that an effective therapy against cancer should target progenitor cells as well as CSC. Experimental studies evidence CSC as responsible for the long-term survival of some type of cancer after therapies, while other experiments are focused on the role of CSC in metastatic progression of


10−4 10−3 10−3 10−3 10−2


23 of 29

Table 3 Values for the densities s and p, and for the radius R at t = 20 obtained for different values of the parameters k0 and k2

s p R

k2 = 1 6.800496391089605E − 12 0.181308411383174 664.938743180467

k0 = 0.01 6.870601868509005E − 12 0.181308411273530 662.669559037713

k0 = 1, k2 = 0.01 2.715795115974555E − 03 0.180838535276993 663.769678859659

R˙ v ∂s ∂s + −r + = kss s − ksp s + s(−kss s − kpp p + kd d) − ks s, ∂t R R ∂r we are able to consider different values of ks . The numerical scheme of resolution does not change much as the only change in the overall scheme is to consider the new right-hand side in the equation for s. Again, in type 1 are considered the same values than the ones considered in Fig. 1 as a reference. In type 2, we take ks = 0.01. Note that according to Molina-Peña & Alvarez (2012), it should be satisfied ks < kp < km . Variations in the parameter ks have the expected effect on the densities. Note that basically, the consideration of ks > 0 only affects to s. Regarding the parameter kp , we observe that when kp is increased, then, as expected, p decreases and this means a decrease in the radius size as well. On the contrary, if kp is decreased, then p increases and a bigger value for the radius R is observed. In Fig. 8, we study the role of km in the model regarding the behaviour of solutions. Note that in type 2, we consider the extreme value km = 0, instead of km = 1 as in Fig. 1. This choice is not a feasible


cancer (see Collins et al., 2005). Therapies for effective treatment of cancer should act to control or eliminate symmetrical CSC division, and target the CSC subpopulation (Collins et al., 2005; Dingly & Michor, 2006; Ganguly & Puri, 2007). In Fig. 6, we present some numerical results for t = 20 that illustrate how changing the values of kss and ksp , in terms of k0 , k1 , and k2 , k3 , affects to the solutions. That is to say, we study the effect of symmetric division of CSC. The variables s, p, m and d, present in types 2, 3, 4 and 5, the expected behaviour. Note that type 2, where we consider k2 = 1 instead of k2 = 0.01, and type 3, where we choose k0 = 0.01 instead of k0 = 1, show overlapping curves. In Table 3, we present the results obtained for the densities s and p, and for the radius R at t = 20, with k2 = 1 and k0 = 0.01, in order to appreciate the similar effect that causes on one hand reducing kss (k0 = 0.01), and on the other hand increasing ksp (k2 = 1). We include results for the reference case with k0 = 1 and k2 = 0.01, to make easier the comparison of the results. Note that when k0 = 0.01 and k2 = 1, the values for s are almost null, the values for p do not change significantly and the same happens to the radius R. The numerical results show that the parameters relative to the symmetric division of CSC, kss and p ks play an important role in the reduction of s, the density of CSC. The regulation of symmetric or asymmetric division is a complex process which depends on a range of conditions (see Morrison & Kimble, 2006; Boman et al., 2007; Dick, 2008; Turner et al., 2009) and the ratio between symmetric and asymmetric CSC division might be possibly crucial for the development and progression of cancer (see Morrison & Kimble, 2006; Pacheco et al., 2009). In the model studied here, there is no parameter which measures the rate of asymmetric division of CSC, so may be it would be interesting incorporating it in future analysis. In Fig. 7, we study the role of the parameters ks and kp in the model. The results correspond to t = 20. Note that ks , the rate at which CSC die is not considered in the model introduced in Tello (2013) and studied here, but changing the equation for s in the following way:

24 of 29

A. I. MUÑOZ ET AL.

Table 4 Values for the densities s, p, m and d, and for the radius R at t = 20 obtained with different values of the parameter km km = 1 2.715795115974555E − 03 0.180838535276993 0.392858050699382 0.423587618907385 663.769678859659

km = 0 2.715795115974555E − 03 0.180838535276993 0.797805357577724 1.864031202904602E − 002 663.769678859659

km = 3 2.715795115974555E − 03 0.180838535276993 0.194952070235556 0.621493599371281 663.769678859659

one as we consider kp = 0.1, and therefore it is not satisfied that kp < km . Anyway, we observe that this fact does not imply significant variations in the densities s and p. With regards to the values for the radius R, we obtain the same effect, insignificant changes. In type 3, we consider km = 3. We see that there exists remarkable variations only in the densities for d and m as expected. In Table 4, we present the results obtained for s, p and the radius R at t = 20 for the different values of the parameter km considered in Fig. 8. In Fig. 9, we consider the following initial conditions: s0 (ri ) = 0.002 + 0.001 cos(π i/10),

p0 (ri ) = 0.198 + 0.1 cos(π i/10),

m0 = 0.8, d0 = 0, i = 1, . . . , N. Note that s0 and p0 have not an homogeneous profile, i.e. their values depend on r, and also that it is not satisfied that for all ri , s0 (ri ) + p0 (ri ) + m0 (ri ) + d0 (ri ) = 1, not even the mean values of s0 , p0 , m0 and d0 add up to 1. The set of parameters considered in this case is the same than the one considered in Fig. 1. The results for t = 0, t = 0.05 t = 1 and t = 20 are shown. We observe that oscillations in s and p disappear for t > T1 , for a certain T1 < 1 and the property of flat profile (no spatial dependence) is then satisfied for all t > T1 . The same effect is observed in the results for the artificial velocity u. Note that there exist T2 < 5, such that the sum of the densities reaches and maintain the value 1 for t > T2 . In Fig. 10 on the left, we present the results obtained for s at t = 40, considering the same set of parameters and initial data as in Fig. 9, with the trapezoidal approximation (3.12) for Δr = 10−2 , and with the implicit midpoint approximation (3.13) for Δr = 10−3 and Δr = 10−2 . We can observe that when one considers Δr = 10−2 with the implicit midpoint scheme, then small oscillations remain in the solution, but when one considers a smaller enough step of discretization, in this case, Δr = 10−3 , then the oscillations disappear. Regarding the long-time evolution of solutions, according to the numerical results obtained in the numerous simulations we have performed, not only the ones presented here and for initial data and values of the parameters under biological restrictions, but for a wide variety of them, we could conjecture that solutions to the system (2.8–2.13) stabilize to homogeneous steady states whose values depend on the values of the parameters and not on the initial data. This kind of conclusion is in accordance with the result presented in the paper by Tello (2013), where it is analytically proved, for a simplified version of (2.8–2.13) and under certain constraints on the parameters, that there exists of a unique homogeneous steady state which is stable and depends on the values of the parameters. We have also obtained that for homogeneous initial data satisfying s0 + p0 + m0 + d0 = 1, solutions satisfy s + p + m + d = 1 for t > 0. If this is not the case, then s + p + m + d tends to 1 as t increases. In Table 5, we present the values to which the solutions stabilize (homogeneous steady states) for the reference case (parameters


s p m d R


25 of 29

Table 5 Homogeneous steady states to which densities tend towards for parameters and initial data as in Fig. 1, for the case considered in Fig. 9 and for the case considered in Fig. 4, type 2 (kpm = 4.32) Densities

Figure 9

Figure 4, type 2 (kpm = 4.32)

0.196849218807238 0.147579041712783 0.318044864758442 0.337526874725320

0.196849218807238 0.147579041712783 0.318044864758441 0.337526874725319

0.401134935825285 0.108851389354070 0.237271461306435 0.252742213515063

Table 6 Homogeneous steady state, steady state and solution to the PDEs original system, to the ODEs system and to the system of algebraic equations, respectively, for the parameter values considered in Fig. 1 Densities s p m d

PDEs system 0.196849218807238 0.147579041712783 0.318044864758442 0.337526874725320

ODEs system 0.197588113446935 0.147448698113865 0.317756847964989 0.337206340474773

System of algebraic equations 0.197588113459191 0.147448698111669 0.317756847960044 0.337206340469096

and initial data as in Fig. 1), for the case considered in Fig. 9, and the case considered in Fig. 4, type 2 (kpm = 4.32). Note that the values for the fraction of the CSC subpopulation are of the same order than the ones of the other non-stem like subpopulations of cells, becoming the CSC population in the case for which kpm = 4.32, the predominant type of cells within the tumour. These results could evidence that a wide variety of fractions for the CSC populations could be expected (see Johnston et al., 2010 and references therein). Unlike normal stem cells that represent a very small population of cells in the tissue they reside in, putative CSC can make a rather significant portion of the tumour cells. For example, in breast cancer, CSC with CD24−/low /CD44+ phenotype constituted 12–60% of the tumour cells. In fact, based on recent mathematical calculation determining the probability of tumour growth after treatment, the assumption that CSC represent a small fraction of all CC in an advanced stage could not appear to be correct (see Shipitsin & Polyak (2008)). In Fig. 10, on the right, we show the steady states to which solution stabilizes for different scenarios as well. To precise, the steady state s = 0.1968492188072 in type 1 is obtained considering the trapezoidal approximation (3.12) with Δr = 10−2 and initial data as in Fig. 9. In type 2, the steady state s = 0.1968492188685 is obtained with the implicit midpoint scheme (3.13) and Δr = 10−2 and initial data as in Fig. 1. And in type 3, the steady state, which takes values between sl = 0.1968492188665 and su = 0.1968492189600, is obtained with the implicit midpoint method and Δr = 10−2 and initial data as in Fig. 9. We observe that with the midpoint scheme, if the step in space is not small enough, undesirable oscillations appear when the initial data are not homogeneous even in the steady state. The values to which solutions to (2.8–2.13) stabilize are in accordance with the steady states that are obtained when one considers the system of ODEs which results when no spatial dependence of the variables is assumed, and also with the solutions to the system of algebraic equations consisting of the r.h.s. of (2.8–2.11). So, we could conjecture that the long-time evolution of solutions to the PDEs system (2.8–2.13) could be predicted by solving the corresponding systems of ODEs or algebraic equations. In order to illustrate the conclusion mentioned above, we show, in Table 6, the homogeneous steady state, the steady state and solution for the PDEs system, the ODEs system and algebraic system,


s p m d

Figure 1

26 of 29

A. I. MUÑOZ ET AL.

respectively, obtained for the set parameters considered in Fig. 1. The ODEs system was solved using the forward Euler method and the system of algebraic equations was solved with a fixed point algorithm. 5. Concluding remarks


The CSC hypothesis postulates that there is a fraction of cells in a tumour, the so-called CSC that is uniquely able to initiate and re-initiate heterogeneous tumours at a secondary site or after cytotoxic treatment. Cancer therapy cannot be successful unless it eliminates the CSC subpopulation, thus understanding the nature and evolution of CSC in a cancer is crucially important for the development of effective therapies. Mathematical modelling is a well-established method to render understanding about complex biological systems. In this paper, we have solved numerically the model of tumour growth presented and mathematically analysed in Tello (2013). The purpose of its numerical resolution is to gain some insight into the relation of different processes leading to exponential growth in solid tumours and into the evolution of different subpopulations of cells within the tumour. The method of resolution we use to solve numerically the model proposed in Tello (2013) is based on the scheme described in the paper Bercovier et al. (1983) which consists in a combination of the characteristics method and the reformulation on the convection term in terms of the total derivative, and a finite element method. To the best of our knowledge, this numerical scheme has not been used previously in the literature regarding the numerical resolution of models of tumour growth of the same nature as the one solved here (see Friedman, 2004), being spectral methods and finite difference schemes the most commonly used (see, e.g. Barrea & Turner, 2005; Wang & Tian, 2008). Experimentally studying the balance between different cellular events involved in the process of tumour growth is not easy. Due to the fact that parameters represent the rates at which reactions, such us cell division, differentiation or death, occur, their tuning in order to obtain numerical simulations with solutions resembling experimental observations, offers a measurement of the role and importance of different process. From the simulations presented here for the early stage in the evolution of the subpopulations for different sets of parameters, one can derive that the results are in accordance with the role of the parameters in the different equations and therefore with the way in which they might affect the concentration of the different subpopulations of cells. In order to bring some light to the fundamental biology of the tumours growth when CSC hypothesis is assumed, we performed an exploration of the parameters involved in the model presented in Tello (2013). We take as a reference the parameter values considered in Molina-Peña & Alvarez (2012), and then, we consider different sets of parameters consisting on changing one of the parameters and remaining the others with their original values. We also consider initial densities according to the ones considered in Molina-Peña & Alvarez (2012). In doing this, we are able to study the role of every parameter in the progression of the tumour. It is important the fact that the parameter relative to the asymmetric division of the CSC into CSC and progenitor cells rate which is taken into consideration in Molina-Peña & Alvarez (2012), is not included in the model studied here, so we are not able to study the effect of encouraging CSC asymmetric division which is considered to play a very important role for fighting against tumours (see Dingly & Michor (2006), Caussinus & Hirth (2007), Giebel & Bruns (2008)). Anyway, we obtained numerical results that are in correspondence with the ones presented in Molina-Peña & Alvarez (2012), and we are able to derive important and interesting conclusions by analysing our numerical results which we shall comment next. Therapies for effective treatment of cancer should act to control or eliminate symmetrical CSC division, and target the CSC subpopulation (Collins et al., 2005; Dingly & Michor, 2006; Ganguly & Puri, 2007; Dylla et al., 2008; Maitland & Collins, 2008). The numerical results show that the role of the


27 of 29

Funding The author is partially supported by Research Projects MTM2011-26016 (MICINN) and TEC201239095-C03-02. References Barrea, A. & Turner, C. (2005) A numerical analysis of a model of growth tumor. Appl. Math. Comput., 167, 345–354.


parameters relative to the symmetric division of CSC, kss and ksp play an important role in the reduction of the density of CSC (see Fig. 6). When increasing the values of ksp or decreasing the values of kss , we observe that s decreases, however this fact does not imply a reduction of the radius of the tumour as the density of transit progenitor cells, p, increases. We have also considered the role of the CSC death rate, which originally was not considered in the model by Tello (2013), but that we have included here. It is observed that this rate, as expected, originates a decrease of s, but again, has no significant effect on the size of the tumour. From the numerical simulations, we would like to highlight the fact that targeting the progenitor subpopulation could be an effective measure to detain tumour progression as we can observe in the results shown in Figs 4, 5 and 7. Note that the fact of increasing the values of the rate kpm (rate at which P cells differentiate to M cells) and the rate kp (rate at which P cells die) has the effect of reducing the size of the tumour. Regarding the values of s, if we increase the rate kpp (rate at which P cells selfrenew with a decreased capacity compared with CSC) then s decreases. These conclusions have been suggested before in the literature (see Kern & Shibata, 2007). Therefore, regarding the evolution of the different subpopulations of cells at the early stage of the tumour, the results that we obtain are in accordance with the ones presented in Molina-Peña & Alvarez (2012). So, according to our numerical results, we obtain the same kind of conclusion that is commented in Molina-Peña & Alvarez (2012), which is that an effective therapy should combine the targeting of CSC population with the targeting of progenitor cells in order to fight against tumours. We go further in our analysis in order to obtain information about the behaviour of the solutions for longer times (see Figs 3, 10 (right) and Tables 5, 6). The numerical results suggest the stabilization of the densities to certain constant values which depend on the set of parameters. The results show the existence of a very rapid adjustment to a pseudo-equilibrium stage, where the proportion of CSC remains small, which is followed by a slower change to a final equilibrium stage. The most remarkable conclusion that we could derive from our results is that the fraction of CSC may vary significantly during the progression of cancer, becoming on some occasions the predominant type of cells within the tumour. These results that we obtain could be admissible and compatible with the CSC hypothesis according to recent research (see, e.g. Shipitsin & Polyak, 2008; Johnston et al., 2010; Enderling et al., 2013; Hillen et al., 2013 and references therein). Finally, we comment some other interesting results that we have obtained from the mathematical point of view. Note that the densities remain constant for a fixed time, in other words there is no spatial dependence of the variables for each time, when one considers homogeneous initial data. The radius of the tumour follows an exponential growth as expected. According to our simulations, we can conjecture that solutions to the system (2.8–2.13) stabilize to homogeneous steady states whose values depend only on the values of the parameters. This fact is in agreement with the result presented and proved in Tello (2013). We have also obtained that the conservation property is satisfied for all the times when homogeneous initial data satisfying s0 + p0 + m0 + d0 = 1 are considered. When this is not the case, we observe that s + p + m + d tends to 1 as t increases.

28 of 29

A. I. MUÑOZ ET AL.


Bayada, G., Chambat, M. & Vázquez, C. (1998) Characteristics method for the formulation and computation of free boundary cavitation problem. J. Comput. Appl. Math., 98, 191–212. Bellomo, N. & Preziosi, L. (2000) Modelling and mathematical problem related to tumor evolution and its interaction with the immune system. Math. Comput. Model., 32, 413–452. Bercovier, M., Pironneau, O. & Sastri, V. (1983) Finite elements and characteristics for some parabolichyperbolic problems. Appl. Math. Model., 7, 89–96. Beretta, E., Capasso, V. & Morozova, N. (2012) Mathematical modelling of cancer stem cells population behavior. Math. Model. Nat. Phenom., 7, 279–305. Boman, B. M., Wicha, M. S., Fields, J. Z. & Runquist, O. A. (2007) Symmetric division of cancer stem cells-a key mechanism in tumor growth that should be targeted in future therapeutic approaches. Clin. Pharmacol. Ther., 81, 893–898. Borsi, I., Fasano, A., Primicerio, M. & Hillen, T. (2014) Mathematical properties of a non-local integro-PDE model for cancer stem cells, Submitted. Byrne, H. M. (2009) Mathematical modeling of solid tumor growth: from avascular to vascular, via angiogenesis. Mat. Biol., 14, 219–287. Calvo, N., Díaz, J. I., Durany, J., Schiavi, E. & Vázquez, C. (2002) On a doubly nonlinear parabolic obstacle problem modelling ice sheet dynamics. SIAM J. Appl. Math., 63, 683–707. Calvo, N., Durany, J. & Vázquez, C. (1999) Numerical approach of temperature distribution in a free boundary model for polythermalice sheets. Numer. Math., 83, 557–580. Caussinus, E. & Hirth, F. (2007) Asymmetric stem cell division in development and cancer. Prog. Mol. Subcell Biol., 45, 205–225. Chaplain, M. A. J. (2008) Modeling aspects of cancer growth: insight from the mathematical and numerical analysis and computational simulation. Multiscale Problems in the Life Sciences. Lecture Notes in Mathematics, vol. 1940. Berlin: Springer, pp. 147–200. Collins, A. T., Berry, P. A., Hyde, C., Stower, M. J. & Maitland, N. J. (2005) Prospective identification of tumorigenic prostate cancer stem cells. Cancer Res., 65, 10946–10951. Cui, S. & Friedman, A. (2003) A hyperbolic free boundary problem modeling tumor growth. Inter. Free Bound., 5, 159–181. Dick, J. E. (2008) Stem cell concepts renew cancer research. Blood, 112, 4793–4807. Dingly, D. & Michor, F. (2006) Successful therapy must eradicate cancer stem cells. Stem Cells, 24, 2603–2610. Dingly, D., Traulsen, A. & Michor, F. (2007) Asymmetric stem cell replication and cancer. PLoS Comput. Biol., 3, e53. Dylla, S. J., Beviglia, L., Park, I. K., Chartier, C., Raval, J. et al. (2008) Colorectal cancer stem cells are enriched in xenogenetic tumors following chemotherapy. PLoS ONE, 3, e2428. Enderling, H., Anderson, A. R. A., Chaplain, M. A. J., Beheshti, A., Hlatky, L. & Hahnfeldt, P. (2009) Paradoxical dependencies of tumor dormancy and progression on basic cell kinetics. Cancer Res., 69, 8814–8821. Enderling, H., Hlatky, L. & Hahnfeldt, P. (2013) Cancer stem cell: a minor cancer subpopulation that redefines global cancer features. Front. Oncol., 3, 76.10.3389/fonc.2013.00076. Friedman, A. (2004) A hierarchy of cancer models and their mathematical challenges. Discrete Contin. Dyn. Syst. Ser. B, 4, 147–159. Friedman, A. (2006) Cancer models and their mathematical analysis. Tutorials in Mathematical Biosciences, III. Lecture Notes in Mathematics, vol. 1872. Berlin: Springer, pp. 223–246. Friedman, A. (2007) Mathematical analysis and challenges arising from models of tumor growth. Math. Model. Methods Appl. Sci., 17, 1757–1772. Friedman, A., Tian, J. P., Fulci, G., Chiocca, E. A. & Wang, J. (2006) Glioma virotherapy: the effects of innate immune suppression and increased viral replication capacity. Cancer Res., 66, 2314–2319. Ganguly, R. & Puri, I. K. (2007) Mathematical model for chemotherapeutic drug efficacy in arresting tumor growth based on stem cells hypothesis. Cell Prolif., 40, 338–354.


29 of 29


Giebel, B. & Bruns, I. (2008) Self-renewal versus differentiation in hematopoietic stem and progenitor cells: a focus on asymmetric cells division. Curr. Stem Cell Res. Ther., 3, 9–16. Hao, W., Hauenstein, J. D., Hu, B., Liu, Y., Sommese, A. J. & Zhang, Y.-T. (2012) Continuation along bifurcation branches for a tumor model with a necrotic core. J. Sci. Comput., 53, 395–413. Hillen, T., Enderling, H. & Hahnfeldt, P. (2013) The tumor growth paradox and immune system-mediated selection for cancer stem cells. Bull. Math. Biol., 75, 161–184. Johnston, M. D., Maini, P. K., Chapman, S. J., Edwards, C. M. & Bodmer, W. F. (2010) On the proportion of cancer stem cells in a tumor. J. Theor. Biol., 266, 708–711. Kern, S. E. & Shibata, D. (2007) The fuzzy math of solid tumor stem cells: a perspective. Cancer Res., 67, 8985–8988. Maitland, N. J. & Collins, A. T. (2008) Prostate cancer stem cells: a new target for therapy. J. Clin. Oncol., 26, 2862–2870. Molina-Peña, R. & Álvarez, M. M. (2012) A simple mathematical model based on the cancer stem cell hypothesis suggests kinetic commonalities in solid tumor growth. PLoS ONE, 7, e26233. Morrison, S. J. & Kimble, J. (2006) Asymmetric and symmetric stem-cell divisions in development of cancer. Nature, 441, 1068–1074. Pacheco, J. M., Traulsen, A. & Dingli, D. (2009) The allometry of chronic myeloid leukemia. J. Theor. Biol., 259, 635–640. Shipitsin, M. & Polyak, K. (2008) The cancer stem cell hypothesis: in search of definitions, markers, and relevance. Pathobiol. Focus, Lab. Invest., 88, 459–463. Tello, J. I. (2013) On a mathematical model of tumor growth based on cancer stem cells. Math. Biosci. Eng., 10, 263–278. Turner, C., Stinchcombe, A. R., Kohandel, M., Singh, S. & Sivaloganathan, S. (2009) Characterization of brain cancer stem cells: a mathematical approach. Cell Prolif., 42, 529–540. Valenciano, J. & Chaplain, M. A. J. (2003) Computing highly accurate solutions of a tumour angiogenesis model. Models Methods Appl. Sci., 196, 747–766. Wang, J. & Tian, J. (2008) Numerical study for a model of tumor virotherapy. Appl. Math. Comput., 196, 448–457. Zheng, X., Wise, S. M. & Cristini, V. (2005) Nonlinear simulation of tumor necrosis, neovascularization and tissue invasion via an adaptative finite-element/level-set method. Bull. Math. Biol., 67, 211–259.

Emergent properties of a computational model of tumour growth.

Oxygen-Driven Tumour Growth Model: A Pathology-Relevant Mathematical Approach.

Numerical solution of a microbial growth model applied to dynamic environments.

Numerical Simulation of a Tumor Growth Dynamics Model Using Particle Swarm Optimization.

A study of tumour growth based on stoichiometric principles: a continuous model and its discrete analogue.

Posttraumatic orbital emphysema: a numerical model.

Warthin's tumour - Resolution following FNA.

A delta-rule model of numerical and non-numerical order processing.

Oxygen transport in a three-dimensional microvascular network incorporated with early tumour growth and preexisting vessel cooption: numerical simulation study.

A nonlinear numerical model of percutaneous drug absorption.

Numerical geometry of map and model assessment.

The immunoreceptor NKG2D promotes tumour growth in a model of hepatocellular carcinoma.

A dynamical model of tumour immunotherapy.

Keratinocyte growth factor promotes epithelial survival and resolution in a human model of lung injury.

A proposed fractional-order Gompertz model and its application to tumour growth data.

Numerical analysis of resolution enhancement in laser scanning microscopy using a radially polarized beam.

A Validated Multiscale In-Silico Model for Mechano-sensitive Tumour Angiogenesis and Growth.

A Numerical Method for Estimating the Variance of Age at Maximum Growth Rate in Growth Models.

A commercially available preparation of Staphylococcus aureus bio-products potently inhibits tumour growth in a murine model of mesothelioma.

FEM numerical model analysis of magnetic nanoparticle tumor heating experiments.

Modified Numerical Simulation Model of Blood Flow in Bend.

Numerical model of dog mast cell tumor treated by electrochemotherapy.

Exogenous administration of protease-resistant, non-matrix-binding IGFBP-2 inhibits tumour growth in a murine model of breast cancer.

Numerical simulation of alga growth and control in Dalian Bay.