Appl. Math. Inf. Sci. 7, No.

5, 1857-1863 (2013) 1857

Applied Mathematics & Information Sciences
An International Journal

http://dx.doi.org/10.12785/amis/070524

On a Model for the Growth of an Invasive Avascular

Tumor
J. A. J. Avila1,∗ and G. Lozada-Cruz 2
1 Department of Mathematics and Statistics, Universidade Federal de São João del Rei, 36307-904, São João del Rei, MG, Brazil
2 Department of Mathematics,IBILCE, Universidade Estadual Paulista, 15054-000, São José do Rio Preto, São Paulo, Brazil

Received: 17 Jan. 2013, Revised: 20 May. 2013, Accepted: 21 May. 2013

Published online: 1 Sep. 2013

Abstract: This paper is devoted to study the 1D model of invasive avascular tumor growth, which takes into account cell division,
death, and motility, proposed by Kolobov and collaborators in 2009. First, we examine the existence and uniqueness of the solution to
this model. Second, we studied qualitatively and numerically the traveling wave solutions. Finally, we show some numerical simulations
for the cell density and nutrient concentration.

Keywords: Invasive avascular tumor, traveling wave, existence of solution, numerical simulation.

1 Introduction representation of an avascular tumor with necrotic region

in the center. The mathematical tool that models the
Cancer is a leading cause of death worldwide and
accounted for 7.6 million deaths (around 13% of all
deaths) in 2008, [1]. On the other hand, in Brazil, the
National Cancer Institute (NCI) released estimates of
cancer incidence for 2012, which will be valid also for
2013. They point to the occurrence of approximately
518,510 new cases of cancer, [2]. The tumors principals
in males are due to non-melanoma skin cancer and
prostate cancer. For females stand out from the
non-melanoma skin cancer and breast cancer.
It is considered a tumor as a colony of cancer cells Fig. 1: Schematic illustration showing the 2D avascular tumor.
(live and dead) that grow uncontrollably, surrounded by
normal tissue. Among the various types of tumors, the growth of tumors are the partial differential equations of
most common are: solid tumors and invasive tumors, the parabolic type. In the case of invasive avascular tumors
difference is in the consistency of the expansive power. are considered cell division, death and motility, as
Solid tumors grow expanding as a compact mass with essential variables in the dynamics of cell density and
well-defined edge between cancer cells and normal tissue, nutrient concentration governed by reaction-diffusion
however, the invasive tumors grow rapidly expanding but equation.
the fraction of cancer cells in the tissue is low. The modeling of the avascular tumor growth is the
The tumors are generally vascular, i.e, there are a network first step towards building fully vascularized tumor
of blood vessels in it, with great possibilities to develop models. Some references regarding the qualitative
angiogenesis. Initially these tumors can be considered analysis of the dynamics of solid avascular tumor growth
without the presence of blood vessels, which are called are found in [3, 4]. Sherratt and Chaplain [5] formulate a
avascular tumors. Thus, the concentration of nutrients new mathematical model for avascular tumor growth.
(oxygen) diffuses to tumor cells from remote enough Ferreira et al. [6] use the equation of reaction-diffusion
blood vessels. In Figure 1 we showed a schematic model for the growth of a tumor avascular with numerical
∗ Corresponding author e-mail: avila jaj@ufsj.edu.br
c 2013 NSP
Natural Sciences Publishing Cor.
1858 J. A. J. Avila , G. Lozada-Cruz: On a Model for the Growth of...

results. Jiang et al. [7] provide a multiscale model (levels: This model was obtained with the following
cellular, subcellular and extracellular) for growth of assumptions:
avascular tumors with numerical results. Roose et al [8]
1.The surroundings of the tumor, corresponding to
show an review outline of a number of illustrative
normal tissue, do not prevent the movement of live
mathematical models describing the growth of avascular
tumor cells and proliferation.
tumors. Bresch et al. [9] reported a viscoelastic model for
the growth of a avascular tumor that describes the 2.The growth of the tumor, in normal tissue, is without
evolution of three components: sane tissue, cancer cells the development of capillary networks.
and extracellular medium. Kolobov et al. [10] study the 3.Although cell division requires a large variety of
autowaves in a model of invasive avascular tumor growth nutrients, this model assumes that the oxygen is the
with 1-d numerical results. Again, Bresch et al. [11] only nutrient and, the missed him, causes the death of
studied the growth an solid avascular tumor in two and live tumor cells.
three dimensions with a focus on numerical methods. 4.The diffusion of oxygen begins in blood vessels distant
This work continues the studies of Kolobov et al. [10] from the tumor, as shown in Figure 1.
adding the part of mathematical analysis: existence and
uniqueness of the solution, and continues dependence of
the initial data. Finally, we present other numerical 3 Mathematical analysis
results, also.
The organization of this work is as follows. In Section In this section we study the existence and uniqueness of
2 we present the mathematical model for invase avascular the solution, and continues dependence of the initial data.
tumor growth. Section 3 is devoted to show the problem Also, we analyze qualitatively the solutions traveling wave
(1) - (2) is well posed in a suitable function space and, in type.
particular we study the traveling waves. Finally, in Section
4, we present some numerical simulations for the problem
(1) - (2). 3.1 Well-posed problem
We consider the system (1) with initial conditions
2 The model a(x, 0) = a0 (x) > 0, s(x, 0) = s0 (x) > 0, ∀x ∈ R. (3)
According Kolobov et al. in [10] the system of nonlinear For 1 6 p < ∞ we define the Banach space Xp by
partial differential equations that models the growth of an
one-dimensional invasive avascular tumor is given by the Xp = {u : R → R; u is bounded, uniform continuous and
reaction-diffusion equation u ∈ L p (R)}
 endowed with the norm
at = Da axx − P(s)a + Ba, −∞ < x < +∞, t > 0

 Z +∞ 1/p
(1) kuk p = sup |u(τ )| + |u(τ )| p d τ .
τ ∈R −∞

s = D s − qa,
t s xx −∞ < x < +∞, t > 0
subject with the boundary conditions For p > q we have the continuous inclusion W 1,q (R) ⊂
(n)
(a, sx ) → (0, 0) when x → −∞ Xp . Now, let n ∈ N and denote by Xp the Banach space
(2) (n)
Xp = {u ∈ Xp : u′ , u′′ , · · · , u(n) ∈ L p (R)} with the norm
(a, s) → (0, 1) when x → +∞,
where n
∑ ku( j) k p , u ∈ Xp
(n) (n)
kuk p = .
a = a(x,t) : cell density (density of live tumor cells) j=0
s = s(x,t) : nutrient concentration The equation (1) joint with (3) can be rewrite in the
and the parameters are defined by space X = Xp × Xp as an abstract evolution equations
(
Da : diffusion coefficient of tumor cells ut = Au + F(u)
Ds : diffusion coefficient of oxygen (4)
u(0) = u0 ,
B : division rate of live cells
Pm : maximal death rate where F : X → X is a nonlinear function given by F(u) =
(−P(s)a, −qa), u = (a, s), u0 = (a0 , s0 ) and A : D(A) ⊂
scrit : critical nutrient concentration X → X is given by
ε : characteristic deviation of s from scrit (2)
q : nutrient consumption rate of live tumor cells D(A) = {(a, s) ∈ X : a, s ∈ Xp , (a, sx ) → (0, 0), x → −∞,
Pm s − scrit (a, s) → (0, 1), x → +∞}
P(s) : cell death rate. P(s) = [1 − tanh( )]. Au = (Da axx + Ba, Ds sxx ).
2 ε

c 2013 NSP
Natural Sciences Publishing Cor.
Appl. Math. Inf. Sci. 7, No. 5, 1857-1863 (2013) / www.naturalspublishing.com/Journals.asp 1859

Using variation of constants formula, the solution of For any (t0 , φ0 ) ∈ [0, ∞) × X we have
(4) is given by
kT (t)(φ ) − T (t0 )φ0 k 6 kT (t)φ − T (t)φ0 k +
Zt (8)
kT (t)φ0 − T (t0 )φ0 k
(t−τ )A
tA
u(t) = e u0 + e F(u(τ ))d τ , t > 0 (5)
0 From the inequality (8), follows that T (t)φ is continuous
at (t0 , φ0 ), and thus the solution u(x,t, φ ) of (1) depends
where continuously on respect to initial data.
  From the definition of {T (t) : t > 0} it follows
T1 (t) 0 immediately that
etA =
0 T2 (t) i) T (0)u0 = u(·, 0, u0 ) = u0 .
ii) T (t1 +t2 )u0 = T (t1 )T (t2 )u0 , ∀t1 ,t2 > 0, this property is
Z
(T1 (t)a0 )(x) = eBt Γ1 (t, x − y)a0 (y)dy a consequence of uniqueness of solution of (1).
Z
R iii) The map [0, ∞) × X ∋ (t, u0 ) 7→ T (t, u0 ) = u(t, u0 ) ∈ X
(T2 (t)s0 )(x) = Γ2 (t, x − y)s0 (y)dy. is continuous, this follows of inequality (8).
R

The functions Γ1 and Γ2 are the heat kernel associated 3.2 Traveling waves
with the parabolic equations at = Da axx and st = Ds sxx
respectively. A travelling wave of (1), (2) is a solution of the type
(a(x,t), s(x,t)) = (A(ξ ), S(ξ )) ∈ C2 (R) × C2 (R) where
Note that the equations at = Da axx + Ba and ξ = x − ct, the functions A and S are the profile of the
st = Ds sxx generate the linear semigroups T1 (t) and T2 (t), travelling wave and c ∈ R is the speed of propagation of
respectively. the wave. The existence of this kind of solutions can be
Thus, the equation (5), for t > 0, take the form found in [12–14] among others.
The functions A and S reduce system (1) in the
Zt following systems of ODEs of second order
a(·,t, u0 ) = T1 (t)a0 − T1 (t − τ )P(s(τ )) a(τ )d τ ,
0
Da A′′ (ξ )+c A′ (ξ )−P(S)A(ξ )+BA(ξ ) = 0

(6)
Zt (9)
Ds S′′ (ξ ) + c S′ (ξ ) − qA(ξ ) = 0
s(·,t, u0 ) = T2 (t)s0 − q T2 (t − τ )a(τ )d τ
0 where −∞ < ξ < +∞, with the boundary conditions
Definition 1.We say that the function u : [0, +∞) → X is a
(
(A, S)(ξ ) → (0, σ ) as ξ → −∞
mild solution of (4) if u ∈ C([0, ∞), X) and satisfies (5). (10)
(A, S)(ξ ) → (0, 1) as ξ → +∞,
We have the following result on existence and
uniqueness of (1)-(2). where σ is a constant corresponding to the limit oxygen
concentration at ξ → −∞.
Theorem 1.Assume that F is global Lipschitz continuous. Note that the points (0, σ ) and (0, 1) are stationary
For any u0 = (a0 , s0 ) ∈ X the system (1) with the points of (9). Now, linearizing the system (9) around the
boundary conditions (2) has a unique mild solution stationary point (0, σ ) we have
u(x,t, u0 ) = (a(x,t, u0 ), s(x,t, u0 )) for all t > 0 with
Da A′′− (ξ )+c A′− (ξ )+(B−P(σ ))A− (ξ ) = 0

u(·, 0, u0 ) = u0 . (11)
′′ (ξ ) + c S′ (ξ ) − qA (ξ ) = 0
Ds S− − −
Proof. From the definition of P and
F(u) = (−P(s)a, −qa), it is not difficult to see that there where −∞ < ξ < +∞. Also, linearizing the system (9)
is L > 0 such that kF(u) − F(v)k 6 Lku − vk for all around the stationary point (0, 1) we have
u, v ∈ X. Thus, as a consequence of the abstract results
Da A′′+ (ξ ) + c A′+ (ξ ) + BA+ (ξ ) = 0

of [12], there exists a unique mild solution for (1). 
′′ (ξ ) + c S′ (ξ ) − qA (ξ ) = 0 (12)
Ds S+ + +
Now, define the family of operators {T (t) : t > 0} on
X by where −∞ < ξ < +∞.
Introducing a vector with the components x1 = A− ,
(T (t)φ )(x) = etA φ (x) x2 = A′− , x3 = S− , x4 = S−
′ , the system (11) can be written

= u(x,t, φ ) (7) as a system of ODEs of first order

= (a(x,t, φ ), s(x,t, φ )), ∀x ∈ R, t > 0. ẋ = Λx (13)

c 2013 NSP
Natural Sciences Publishing Cor.
1860 J. A. J. Avila , G. Lozada-Cruz: On a Model for the Growth of...

where

0 1 0 0
 k+ −
1 = k1 , k+ −
2 = k2
 P(σ )−B − Dca 0 0  √
Λ= Da ,
c(Ds − 2Da ) + Ds ∆ +
 
 0 0 0 1  (14)
q −
− Dcs 2qDa
 
Ds 0 0  
 
x = (x1 , x2 , x3 , x4 )T .  c D − D
c + √∆ +
− 2BD D 
 
 s a a s 
Similarly, putting y1 = A+ , y2 = A′+ , y3 = S+ , y4 = S+
′ , the
k+ 2qD2a
 
3 = (19)
 
system (12) can be written as a system of ODEs of first

 
order
 2Da 
 − √ 
ẏ = Ξy (15) 
 c+ ∆ + 

 
where
  1
0 1 0 0 √
c(Ds − 2Da )− Ds ∆ +
 
− B − c 0 0 
Ξ= Da Da , −
 0 0 0 1  2qDa
 
q
(16)  
0 0 − Dcs  
Ds  c D − D
c − √∆ +
− 2BD D 
 
s a a s 
y = (y1 , y2 , y3 , y4 )T . 
k+ 2qD2a
 
4 =
 
The matrices Λ and Ξ have eigenvalues respectively 


√  2Da 
c −c ± ∆ − − √
 
µ1− = 0, µ2− = − , µ3,4 −
= c+ ∆+
 
 
Ds 2D√a (17)  
−c ± ∆+ 1
µ1 = µ1 , µ2 = µ2 , µ3,4 =
+ − + − +
2Da
Thus, for the system (13), we have the following linear
where independent solutions
∆ − = c2 + 4(P(σ ) − B)Da xi (ξ ) = k− µi ξ
−
, i = 1, · · · , 4
i e (20)
∆ + = c2 − 4BDa
− µ2 ξ −
From (20) we have lim x1 (ξ ) = k−
1 , x2 (ξ ) = k2 e
and the corresponding eigenvectors are given, respectively, ξ →−∞
by − µ4− ξ
and x4 (ξ ) = k4 e are unbounded solutions when
Ds T
k− T
1 = (0, 0, 1, 0) , k−2 = (0, 0, − , 1) ξ → −∞ if c > 0 and P(σ ) > B respectively. And,
c µ3− ξ
√ x3 (ξ ) = k−
3e is bounded when ξ → −∞.
c(Ds − 2Da ) + Ds ∆ −
 
− If µ3− is a real positive (negative) number the stationary

 2qDa 
 point (0, σ ) of (11) is saddle (stable) node and if µ3− is a
complex number the stationary point (0, σ ) of (11) is a
 h 
 D c c − √∆ −
+ 2 P(σ ) − B
D 
 i 
 a s  stable focus.
k−
  And for the system (15), we have the following linear
3 = 2qD2a
  (18)


 independent solutions
2Da
 
µi ξ
+
− √ yi (ξ ) = k+
 
i e , i = 1, · · · , 4. (21)
c+ ∆
 
 − 
From (21) we have lim y1 (ξ ) = k+
1 , lim y2 (ξ ) = 0 and
 
1 ξ →+∞ ξ →+∞
√ µi+ ξ

c(Ds − 2Da )− Ds ∆ −
 the other solutions yi (ξ ) = k+
i e√ , i = 2, 3 also remains
 −
2qDa  bounded when ξ → +∞ if c > 2 BDa . Thus, we have that



 the stationary point (0, 1) of (12) is either a stable node or
 √
h √ i  a stable focus.
− c+ ∆ − c Ds − 2Da + Ds ∆ 


 − −
 
−
 
k4 = 
 4qDa2 

  4 Numerical results
2Da
 
− √
 

 c+ ∆ −

 In this section we will find numerically by finite difference
  method, a solution of type traveling wave, and full solution
1 of our problem (1) - (2).

c 2013 NSP
Natural Sciences Publishing Cor.
Appl. Math. Inf. Sci. 7, No. 5, 1857-1863 (2013) / www.naturalspublishing.com/Journals.asp 1861

4.1 Dimensionless equations

The units of measurements for the cell density and
0.015
nutrients concentration is, respectively,
c [cell/cm3 ] and s [mol/cm3 ].
The values and units of measurement of the parameters are
0.01
taken from Kolobov et al. [10],
Da = 2.5 × 10−8 [cm2 /s], Ds = 2.5 × 10−5 [cm2 /s] A

B = 10−7 [1/s], Pm = 0.2 × 10−6 [1/s]

scrit = 0.3 × 10 [mol/cm ], ε = 10−9 [cm3 /mol]
−7 3
0.005

q = 10−20 [mol/s · cell].

The characteristic scales are:
Length scale : L0 = 5 × 10−2 [cm] 0
-1000 -500 0 500 1000
ξ
6
Time scale : T0 = 10 [s]
Density scale : amax = 107 [cell/cm3 ] Fig. 2: Cell density profile A(ξ ).
Concentration scale : smax = 10−7 [mol/cm3 ].
The independent and dependent dimensionless variables
are given by
x∗ = x/L0 , t ∗ = t/T0 , a∗ = a/amax , s∗ = s/smax 1

We have used the following dimensionless data:

D∗a = 10 , D∗s = 104 , B∗ = 10√
−1 , P∗ = 0.2 , s∗ = 0.3 ,
m crit
0.8

ε ∗ = 10−2 , q∗ = 1 , cmin = 2 2 , c = cmin , σ = 0.18

S
Henceforth, the equations are written in dimensionless
0.6
form.

4.2 Traveling waves 0.4

Rewriting the equation (9) in the following form

(
A′′ = p1 A′ + q1 (S)A 0.2
-1000 -500 0 500 1000
(22)
S′′ = p2 S′ + q2 A ξ

where Fig. 3: Nutrient concentration profile S(ξ ).

c 1

p1 = − , q1 (S) = P(S) − B ,
Da Da
c q
p2 = − , q2 = ,
Ds Ds 4.3 Numerical solution
For the numerical solution of problem (1) - (2) we used
joint with the boundary conditions the classic and unconditionally stable Crank-Nicolson
A(−L) = 0 = A(L) Method. The domain (−L, L) is discretized with uniform
(23) grid of N = 1200 elements and L = 2000. The sizes of
S(−L) = σ , S(L) = 1,
spatial and temporal steps are, respectively, h = 2L/N and
where ξ ∈ (−L, L) and L > 0 is large. k = v0 h/c, where v0 = 0.26. The initial condition is given
2
Taking L = 1000 and using finite difference method by gaussian function a(x, 0) = 0.1e−0.0025x and constant
we have a numerical solutions for cell density and function s(x, 0) = 1.
nutrient concentration which are showed in the Figure 2 The asymptotic behavior of live cell density, at each
and 3, respectively. point where it reaches its maximum value, is showed in

c 2013 NSP
Natural Sciences Publishing Cor.
1862 J. A. J. Avila , G. Lozada-Cruz: On a Model for the Growth of...

the Figure 4. Note that of stationary state is attained in

t ≈ 200, i.e, in approximately 6.5 years. In Figure 5 and 6

0.7
0.5

0.6 s

0.5

2000
200
1500
0.4
150 1000
a

500

100 0
t x
0.3 -500
50 -1000
-1500
0 -2000
0.2

0.1
Fig. 6: The profile of nutrient concentration density for long-time
0
0 50 100 150 200
t

Fig. 4: Asymptotic behavior of maximal live cell density a. 5 Conclusion

The problem of growth of an invasive avascular tumor,
one dimensional, is well-posed. Physically, we observe in
we showed, in 3D, the density profile of live tumor cells figures and in numerical simulations that the tumor grows
and nutrient concentration profile, respectively, the long- in the direction where there are nutrients (oxygen). It is
time. We note that at steady state, the density profile of observed that it is invasive and aggressive from the
live tumor cells approaches the solution of type traveling beginning to about t = 60 (2 years), immediately after a
wave, given, in the Figure 2. relaxation which is going to reach steady state, at
approximately t = 200 (6.5 years), where the density of
live tumoral cells is very low. The tumor grows because
the tumor cells die. Thus, the volume of necrotic region
increases.

References
[1] World Healt Organization: Media centre. Cancer, Fact sheet
0.6 No 297, Reviewed January, (2013).
[2] Estimativa 2012: Incidência de Câncer no Brasil. Instituto
0.4 Nacional de Câncer José Alencar Gomes da Silva,
a Coordenação Geral de Ações Estratégicas, Coordenação de
0.2
Prevenção e Vigilância. Rio de Janeiro: Inca, 118 (2011).
[3] Ward, J. P. King, J. R., Mathematical modelling of avascular-
0 2000
200 1500
tumour growth. IMA Journal of Mathematics Applied in
150
500
1000
Medicine & Biology, 14, 39–69 (1997).
t 100
0
x
[4] Byrne, Helen, M. A weakly nonlinear analysis of a model
-500

50 -1000 of avascular solid tumour growth. Journal Mathematical

-1500
0 -2000
Biology, 39, 59–89 (1999).
[5] Sherratt, Jonathan A.; Chaplain, Mark A. J., A new
mathematical model for avascular tumour growth.
Fig. 5: The profile of live cell density the long-time. Mathematical Biology, 43, 291–312 (2001).
[6] Ferreira Jr., S. C.; Martins, M. L.; Vilela, M. J., Reaction-
diffusion model for the growth of avascular tumor. Physical
Review E, 65, 021907, (2002).

c 2013 NSP
Natural Sciences Publishing Cor.
Appl. Math. Inf. Sci. 7, No. 5, 1857-1863 (2013) / www.naturalspublishing.com/Journals.asp 1863

[7] Jiang, Yi; Pjesivac-Grbovic, J.; Cantrell, Charles; Freyer, J. P., Jorge Andrés Julca
A Multiscale model for avascular tumor growth. Biophysical Avila, PhD in Mechanical
Journal, 89, 3884–3894 (2005). Engineering from USP
[8] Roose, T.; Chapman, S. J.; Maini, P. K., Mathematical - University of São Paulo.
Models of Avascular Tumor Growth. SIAM Review, 49, 179– Assistant Professor at the
208 (2007). Department of Mathematics
[9] Bresch, Didier; Colin, Thierry; Grenier, Emmanuel; Ribba, of UFSJ - Federal University
Benjamin; Saut, Oliver, A viscoelastic model for avascular of São João del-Rei.
tumor growth. Discrete and Continuous Dynamical Systems, Has experience in Applied
supplement, 101–108 (2009).
Mathematics emphasizing
[10] Kolobov, A. V., Gubernov, V. V. and Polezhaev, A. A.,
Numerical Analysis,
Autowaves in a model of invasive tumor growth. Biophysics,
54, 232–237 (2009).
Functional Analysis and PDEs. For the moment, the
[11] Bresch, Didier; Colin, Thierry; Grenier, Emmanuel; Ribba, major interests are in the Navier-Stokes Equations,
Benjamin; Saut, Oliver, Computational modeling of solid Non Newtonian Fluids and Finite Elements. Is also
tumor growth: The avascular stage. SIAM J. Sci. Comput, experienced in Mechanical Engineering with emphasis in
32, 2321–2344 (2010). Fluid Mechanics and Fluid Dynamics. Is member of the
[12] Henry, D., Geometric theory of semilinear parabolic Scientific Committee of UFSJ.
equations, Lect. Notes in Math. Springer, 840, (1981).
[13] Aronson, D. G.; Weinberger, H. F., Multidimensional Germán Jesús Lozada
nonlinear diffusion arising in population genetics. Adv. in Cruz, PhD in Mathematics
Math., 30, 33–76 (1978). (2001) from University
[14] Volpert, Aizik I.; Volpert, Vitaly A.; Volpert, Vladimir A., of São Paulo at São Carlos.
Traveling wave solutions of parabolic systems. Translations Assistant professor at the
of Mathematical Monographs. American Mathematical Department of Mathematics
Society, Providence, RI, 140, (1994). of State University of
São Paulo. Has experience in
Mathematics with emphasis
in Partial Differential
Equations.
The areas of interest are asymptotic behavior of reaction
diffusion equations in dumbbell domains, existence of
attractors and continuity of attractors with relation to
small parameter.

c 2013 NSP
Natural Sciences Publishing Cor.