Mathematical Biosciences 255 (2014) 83–90

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Mathematical modeling of efficient protocols to control glioma growth J.R. Branco a, J.A. Ferreira b,⇑, Paula de Oliveira b a b

CMUC & Polytechnic Institute of Coimbra, ISEC, DFM, Coimbra, Portugal CMUC & Department of Mathematics, University of Coimbra, Coimbra, Portugal

a r t i c l e

i n f o

Article history: Available online 21 July 2014 Keywords: Glioma Viscoelastic behavior Chemotherapy Protocols Numerical simulation

a b s t r a c t In this paper we propose a mathematical model to describe the evolution of glioma cells taking into account the viscoelastic properties of brain tissue. The mathematical model is established considering that the glioma cells are of two phenotypes: migratory and proliferative. The evolution of the migratory cells is described by a diffusion–reaction equation of non Fickian type deduced considering a mass conservation law with a non Fickian migratory mass flux. The evolution of the proliferative cells is described by a reaction equation. A stability analysis that leads to the design of efficient protocols is presented. Numerical simulations that illustrate the behavior of the mathematical model are included. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Cancer is a complex disease which leads to the uncontrolled growth of abnormal cells, destruction of normal tissues and invasion of vital organs. There are different stages of tumor development with varying duration, starting from genetic changes at the cell level and finishing with detachment of metastasis and invasion. Tumor cell transport and proliferation are the main contributors to the malignant dissemination [44]. Extensive research has been done to model cancerous growth, specially on solid tumors, in which growth primarily comes from cellular proliferation. It is far beyond the aim of the present paper to list exhaustively the many significant contribution in the topic. Refs. [16,22,23,25,33,43,44] and the references therein represent some of these contributions. Gliomas are diffusive and highly invasive brain tumors accounting for about 50% of all primary brain tumors and, unfortunately, the prognosis for patients with gliomas is very poor. Median untreated survival time for high grade gliomas ranges from 6 months to 1 year and even lower grade gliomas can rarely be cured. Theorists and experimentalists believe that inefficiency of treatments results from the high mobility of glioma cells. Additionally gliomas can exhibit very high proliferation rates. The understanding of malignant glioma growth is still far from being complete, mostly because gliomas proliferate as solid tumors and invade the surrounding brain parenchyma actively.

Proliferation and specially migration of gliomas represent a very challenging problem from a mathematical viewpoint. Cancer research has been a fertile ground for mathematical modeling, beginning with the early concept of simple exponential growth of solid tumors doubling at a constant rate. The introduction of logistic or gompertzian growth (characterized by increased doubling time and decreased growth fraction as a function of time) allowed to slow the growth in the later stages. With the recognition that tumor cells might spread outside the grossly visible mass, invading locally and metastasizing distantly, and that some cells die during the development process, the mathematical concepts necessarily became more complex than those used in the original simple models for solid tumors. The initial answer to the question of how to measure the growth of an infiltrating glioma was provided by Murray in the early 90s [33]. He formulated the problem as a conservation law where the rate of change of tumor’s cell population results from mobility and net proliferation of cells. An equation of type

@c þ r  J F ¼ f ðcÞ in X  ð0; 1Þ @t

was used, where X  Rn ; n ¼ 1; 2; 3, is the glioma domain, cðx; tÞ denotes the tumor cell density at location x and time t; f ðcÞ denotes net proliferation of tumor cells, and r defines the spatial gradient operator. Under the assumption of the classical Fick’s law for the mass flux J F

e rc; JF ¼  D e is the diffusion tensor, Eq. (1) can be written as where D

⇑ Corresponding author. E-mail addresses: [email protected] (J.R. (J.A. Ferreira), [email protected] (P. de Oliveira). http://dx.doi.org/10.1016/j.mbs.2014.07.002 0025-5564/Ó 2014 Elsevier Inc. All rights reserved.

Branco),

[email protected]

ð1Þ

ð2Þ

84

@c e rcÞ þ f ðcÞ in X  ð0; 1Þ: ¼ r  ðD @t

J.R. Branco et al. / Mathematical Biosciences 255 (2014) 83–90

ð3Þ

The mathematical model is complemented by boundary conditions which impose no migration of cells beyond the brain boundary, that is,

J F  g ¼ 0; on the boundary, where g denotes the exterior unit normal to the brain region, and by initial conditions cðx; 0Þ ¼ c0 ðxÞ; x 2 X, where c0 defines the initial spatial distribution of malignant cells. Tumor growth is generally assumed to be exponential, so the cell growth term is given by f ðcÞ ¼ qc, where the net proliferation rate q is constant. Logistic and gompertzian growths have been also considered but found to be unnecessary in the time frames considered for gliomas development [25]. To apply the modeling approach to specific patients, a more realistic look at the brain geometry and structure was necessary. Swanson et al. introduced in [43] the complex geometry of the brain and allowed diffusion to be a function of the spatial variable to reflect the observation that glioma cells exhibit higher motility in the white matter than in the gray matter. The partial differential Eq. (3), of parabolic type, was established combining the mass conservation law (1) with Fick’s law (2) for mass flux. It is well known that, in this case, if a sudden change on the cell concentration takes place somewhere in the space, it will be felt instantaneously everywhere. This means that Fickian approach gives rise to an infinite speed of propagation which is not physically observable. To avoid the limitation of Fickian models an hyperbolic correction has been proposed in different contexts (see [9,28,14,15,24,34] and the references cited in these papers). It is accepted by the biomedical research community that biochemical and biophysical properties of the brain tissue, namely of the extracellular matrix (ECM), are key factors in the proliferation and migration of glioma cells. The aggressiveness of the gliomas is determined by its unique pattern of interaction with ECM. Experimental studies show that the mechanical properties of ECM are regularization factors in the evolution of several cell types in particular glioma cells [12,13,31,45,47]. In fact it was observed that the growth, differentiation and functionalities of glioma cells are determined by the stiffness of the ECM. These observations are explained by the fact that extracellular matrix stiffness induces complex biochemical phenomena that depend on the type of diffusive cells and microenvironment properties which are not yet clarified. The complete understanding of such complex biochemical effect can be used to develop tumor treatments based on the characteristics of the mechanical milieu where the cancer cells move. The aim of the present paper is to study the influence of these properties on glioma growth and treatment. To this end we establish a class of non Fickian models that take into account the viscoelastic behavior of the brain tissue. The mathematical model that we consider is defined in a simple geometry. To apply the modeling approach to specific patients, a more realistic look at the brain geometry and structure is necessary. In this case we can follow [43], where, for a Fickian model, a complex geometry of the brain was considered as well as a space dependent diffusion coefficient was taken into account to reflect the observation that glioma cells exhibit higher motility in the white matter than in gray matter [22]. Finally we observe that the most popular treatments used to combat gliomas are chemotherapy and radiotherapy. Some mathematical models that describe the effect of these treatments were proposed in the literature. Without being exhaustive we mention [29,38,46]. Chemotherapy involves the use of drugs to disrupt the cell cycle and to block proliferation. The success of chemotherapy agents varies widely, depending on cell type and the type of

drug being used. The effectiveness of a particular drug depends on the concentration of drug reaching the tumor, the duration of exposure and the sensitivity of the tumor cells to the drug. Tracqui et al. in [46] incorporated chemotherapy by introducing cell death as a loss term. If GðtÞ defines the rate of cells death then, assuming a loss proportional to the tumor cells density, Eq. (3) is replaced by

@c e rcÞ þ f ðcÞ  GðtÞc in X  ð0; T; ¼ r  ðD @t

ð4Þ

where

GðtÞ ¼



k;

when chemotherapy is being administered

0; otherwise:

ð5Þ

Here k describes the rate of cell death due to exposure to the drug. If f ðcÞ ¼ qc, for a tumor to decrease in size during chemotherapy, k must be larger than the growth rate q of the cell population. The main question in this paper is to define k and the periods of chemotherapy applications that lead to control the glioma mass. In chemotherapy protocols a specific drug or a cocktail of drugs is injected in the circulatory system and are homogeneously distributed in the human body. The effect of chemotherapy in glioma cells is described here by the function GðtÞ defined by (5). However, the death rate of cells that are exposed to the action of a drug should depend on location, duration of exposure and drug concentration. To incorporate all these interveners in the death rate we need to define a function G depending on x; t; cðx; tÞ and cd ðx; tÞÞ. To define such death effect on the tumor cells, the equation for the dynamic of cancer cells needs to be coupled with a diffusion equation for the concentration in the line of the models studied for instance in [26,35,48], where Michaelis–Menten kinetic or its generalization have been considered in the definition of the corresponding G. As in the present paper we do not consider such coupling, our assumption on the death rate GðtÞ while simple, is reasonable and can be mathematically manipulated to allow the definition of chemotherapy protocols with a prescribed effect. Partial differential equations of non-Fickian type that describe the evolution of cells in a medium where they die as they move can be established using the continuous time random walks approach. Without being exhaustive we refer [1,18,42] where such approach was considered in different contexts. This approach was also considered, for instance in [16,17], to establish non-Fickian diffusion models to describe the proliferation and migration of glioma cells in the absence of the death effect. Our aim in this paper is the modeling and analysis of glioma growth under the effect of the rheological properties of the brain tissue. The paper is organized as follows. The brain tissue presents a viscoelastic behavior that can be described by a Voigt–Kelvin model (see for instance [21,27,32]). Using this approach we present in Section 2 a class of non Fickian models to describe the space and time evolution of glioma cancer cells constructed by combining the diffusion process with the viscoelastic properties of the brain tissue. In Section 3 we study the behavior of the glioma mass when chemotherapy is considered. Criteria to define efficient protocols that lead to the decreasing of the tumor mass are established in this section. In Section 4 we introduce a semi-discrete model that mimics the continuous model in the sense that it presents the same qualitative properties. Plots depicting the evolution of gliomas are included in Section 5. The numerical results illustrate the theoretical results obtained. Finally, in Section 6 we present some conclusions. The present paper aims to extend the results obtained by the authors in [4–6].

85

J.R. Branco et al. / Mathematical Biosciences 255 (2014) 83–90

2. A viscoelastic model The class of non Fickian models that we present in what follows is established by taking into account the viscoelastic nature of the brain tissue. Following [8,10,11,20,30,41], if a diffusion process occurs in a medium with a viscoelastic behavior, then a modified diffusion equation

@c e rcÞ þ r  ð D e v rrÞ þ f ðcÞ in X  ð0; 1Þ; ¼ r  ðD @t

ð6Þ

should be used, where r represents the normal stress which is assumed to be the most relevant component of the stress developed e and D ev in the brain tissue when the tumor cells move in it. In (6) D stand for the diffusion and stress driven diffusion tensors respectively. The proliferation and migration of several cancer cells depend on the rigidity of ECM, growing significantly better on stiff matrices than on soft tissues [45,47]. This phenomenon is usually called durotaxis or mechanotaxis and it was firstly defined in [31] when the migration of fibroblasts in vitro from soft to stiff regions of the e v in Eq. (6) is a diagonal ECM was observed. Based on these facts, D tensor with negative entries. We assume that the viscoelastic behavior of the brain tissue is described by the Kelvin-Voigt model

@r @ þ br ¼ a1  þ a2 ; @t @t

ð7Þ

where  stands for the normal strain. Eq. (7) is based on a mechanistic model which is represented by a spring (restorative force component) and a dashpot (damping component) in parallel, connected with a free spring. In (7) the viscoelastic characteristic time b is given by b ¼ E0lþE1 , a1 ¼ El0 E1 ; a2 ¼ E0 where E1 is the Young 1 1 modulus of the spring element, l1 represents the viscosity and E0 stands for the Young modulus of the free spring (see [21,27,32,39]). Eq. (7) leads to the following expression for r

rðtÞ ¼

Z

t



ebðtsÞ a1 ðsÞ þ a2

0



@ ðsÞ ds þ ebt rð0Þ: @t

ð8Þ

If we assume that the strain  satisfies  ¼ kc where k is a positive constant(see [8]), we obtain from (6) and (8) an integro-differential equation

@c ¼ r  ðDrcÞ þ @t  ð0; 1Þ;

Z 0

(proliferation phenotype) the cancer cells do not migrate and only proliferation takes place with rate q; a cell of type 1 remains in state 1 during a random time period and then switches to a cell of type 2; b1 is the switching rate from state 1 to 2; a cell of type 2 remains in state 2 during a random time period and then switches to a cell of type 1; b2 is the switching rate from state 2 to 1. Let uðx; tÞ and v ðx; tÞ represent the density of migratory and proliferation cells at x and t, respectively. The dynamics of glioma cells is then described by (

@u @t @v @t

Rt ¼ r  ðDruÞ þ 0 ker ðt  sÞr  ðDv ruðsÞÞ ds  b1 u þ b2 v in X  ð0; T; ¼ qv þ b1 u  b2 v in X  ð0; T;

ð10Þ where T > 0 is fixed, D and Dv denote square matrices of order n. System (10) is complemented with initial conditions

ker ðt  sÞr  ðDv rcðsÞÞ ds þ f ðcÞ in X ð9Þ

e þ ka2 D e v ; Dv ¼ kða1  ba2 Þ D e v and ker ðsÞ ¼ ebs . where D ¼ D In this paper we consider that the viscoelastic behavior of the brain tissue is described by the Voigt–Kelvin model equation (7) and the mass flux of migration cells J is driven by the gradient of the concentration and by the gradient of the forces exerted by the e rc  D e v rr. The brain tissue into the glioma cells, that is J ¼  D stress r is given by (8) where the strain  is identified with the results of the action of the glioma cells into the brain tissue. Here, to simplify, we assume that such results depend linearly on the glioma cell concentrations. We do not take into account the microenvironment where glioma cells migrate and proliferate, their constituents and their interactions. Mathematical models based on mixture theory and interaction forces between intervenients in the cancer growth have been studied, for instance, in [2,3,37,36,40]. To establish a mathematical model to describe the space–time evolution of gliomas some medical information is needed to establish. According to [16,17] the following assumptions are considered in our model: the glioma cells are of two phenotypes – proliferation (state 1) and migratory (state 2); in state 2 (migratory phenotype) the cells randomly move but there is no cell fission; in state 1

in X;

and boundary conditions

J  g ¼ 0 on @ X;

ð11Þ

where @ X denotes the boundary of X; g represents the exterior unit normal and the non Fickian flux J is given by Rt JðtÞ ¼ DruðtÞ  0 ebðtsÞ Dv ruðsÞ ds. Condition (11) means that the glioma is located inside of the brain and the cancer cells do not cross the pia mater. In what follows we assume that D ¼ ½dij  and Dv ¼ ½dv ;ij  are diagonal matrices with diagonal entries di and dv ;i such that

0 < a0 6 di ; dv ;i in X;

i ¼ 1; . . . ; n:

ð12Þ

Let M1 ðtÞ be the mass of glioma cells in X, R M1 ðtÞ ¼ X ðuðtÞ þ v ðtÞÞ dX. We study in what follows the behavior of M1 ðtÞ. We start by remarking that

M01 ðtÞ ¼

 Z  @u @v ðtÞ þ ðtÞ dX: @t @t X

As u and

M01 ðtÞ

¼

ð13Þ

v are defined by the system (10), from (13) we obtain

Z

ðr  JðtÞ þ qv ðtÞÞ dX;

X

that leads to

M01 ðtÞ ¼ 

t

v ð0Þ ¼ v 0

uð0Þ ¼ u0 ;

Z

JðtÞ  g d@ X þ q

@X

Z

v ðtÞ dX:

X

R From (11) we conclude that M01 ðtÞ ¼ q X v ðtÞ dX; which means that the instantaneous time variation of the cancer mass depends only, as expected, on the mass of the proliferation cells and on the proliferation rate q. Assuming the positivity of u, we finally obtain the upper bound M1 ðtÞ 6 eqt M1 ð0Þ. To avoid the positivity assumption on u, we establish in what follows an upper bound for the mass related functional M 2 ðtÞ ¼ kuðtÞk2 þ kv ðtÞk2 , where k:k denotes the usual L2 norm pffiffiffiffiffiffiffi induced by the usual L2 inner product ð:; :Þ. As M 1 ðtÞ 6 jXj ðkuðtÞk þ kv ðtÞkÞ where jXj represents the measure of omega if   v ðtÞk pffiffiffiffiffi ; kp ffiffiffiffiffi P 1, we conclude that an upper we assume that min kuðtÞk jXj

jXj

bound for M 1 ðtÞ can be deduced from an estimate of M 2 ðtÞ. This assumption is biologically sound because it states that the tumor density is larger than 1.    As 12 M 02 ðtÞ ¼ @u ðtÞ; uðtÞ þ @@tv ðtÞ; v ðtÞ , we obtain from (10) @t Z pffiffiffiffi 1 0 M 2 ðtÞ ¼ JðtÞ  guðtÞ d@ X  k DruðtÞk2 2 @X Z t   ker ðt  sÞDv ruðsÞ ds; ruðtÞ 0

 b1 kuðtÞk2 þ ðb2 þ qÞkv ðtÞk2 þ ðb1 þ b2 ÞðuðtÞ; v ðtÞÞ;

ð14Þ

86

J.R. Branco et al. / Mathematical Biosciences 255 (2014) 83–90

where the inner product in L2 ðXÞ  L2 ðXÞ is denoted by ðð:; :ÞÞ and k:k represents the induced norm. Considering the boundary condition (11), the Cauchy–Schwarz inequality and the following equality

Z t 2 pffiffiffiffiffiffi d ker ðt  sÞ Dv ruðsÞ ds dt 0 Z t  ¼2 ker ðt  sÞDv ruðsÞ ds; ruðtÞ ð15Þ

we deduce from (14) that

ð16Þ

R 2 pffiffiffiffiffiffi t where EðtÞ ¼ M2 ðtÞ þ 0 ker ðt  sÞ Dv ruðsÞ ds . Inequality (16) leads to

M 2 ðtÞ 6 e2 max

b

2 b1 ;b1 b2 þq;b 2 2



t

M2 ð0Þ:

ð17Þ

We observe that if maxfb2  b1 ; b1  b2 þ 2q; 2bg ¼ 2b then b1 þ 2q þ 2b < b2 < b1  2b which is not possible. This means that we can drop 2b from the max expression. In the case b2  b1 > b1  b2 þ 2q we have b2  b1 > q; for b2 < b1 þ q the maximum is 2q  b2 þ b1 > q. In both cases the second member of (17) is an increasing function of t. As expected, under these assumptions, we can not select parameters b2 ; b1 ; q such that M2 ðtÞ is bounded in time. We remark that inequality (17) allow us to conclude the stability of the proposed mathematical model with respect to perturbations of the initial conditions in ½0; T, for fixed T > 0. 3. Chemotherapy: control of the glioma growth To take into account the chemotherapy effect, the viscoelastic model for glioma growth (10) is modified as follows (

@u ¼ @t @v ¼ @t

ð18Þ

where GðtÞ is defined by (5). From (18), following the proof of the upper bound (16), it can be shown that

 b  b1 b  b2 E0 ðtÞ 6 2 max 2  GðtÞ; 1 þ q  GðtÞ; b EðtÞ: 2 2



ð19Þ

In what follows we establish conditions on the parameters that lead to a decreasing of M 2 ðtÞ: 1. If the net proliferation rate is greater than the switching proliferation rate

q > b2  b1 ;

ð20Þ

and the difference between the rate of cells death and the switching proliferation rate is bounded by the viscoelastic characteristic time

GðtÞ 

b2  b1 < b; 2

ð21Þ

then Eq. (19) leads to



M2 ðtÞ 6 e

2

b1 b2 þq 2

 Rt t

0

t

GðsÞ ds
0;

 Z t b1  b2 þq t < GðsÞ ds; 2 0

Z

0

Z t 2 pffiffiffiffiffiffi  2b ker ðt  sÞ Dv ruðsÞ ds ;



 Z t b2  b1 t< GðsÞ ds; 2 0

ð29Þ

we conclude that M 2 ðtÞ decreases. Again we observe that the parameter b has influence on the definition of an admissible threshold for the chemotherapy treatment. We note that as from (27)



 Z t b1  b2 þqþb t > GðsÞ ds; 2 0

ð30Þ

we must impose that the difference between the net and switching proliferation rates is greater than the viscoelastic characteristic time

q  ðb2  b1 Þ > b

ð31Þ

in order to have compatibility between (29) and (30). When chemotherapy is applied, conditions (21) and (23) or conditions (27) and (29) can be used to determine an effective dosage that induces a rate k of cell death, due to the exposure to the drug, that allows to control the total tumor mass. Obviously the value of k depends of the protocol of chemotherapy. The typical bang–bang protocol corresponds to a treatment which alternate maximum doses of chemotherapy with rest periods when no drug is administered, as defined by (5) and illustrated in Fig. 1.

GðsÞ ds

M 2 ð0Þ:

ð22Þ

To conclude that M2 ðtÞ decreases we need to combine condition (20) and (21) with

4. A semi-discrete model To compute the artificial mass M 2 ðtÞ we use a numerical method which is obtained discretizing the spatial derivatives of

87

J.R. Branco et al. / Mathematical Biosciences 255 (2014) 83–90

k

ðwH ; qH ÞXH ¼

Nh Nh 1 X 2 X

xi;j wi;j qi;j ; wH ; qH 2 W H ðXH Þ;

ð37Þ

i¼0 j¼0

t Fig. 1. Chemotherapy protocol.

(18) using centered difference operators. In what follows we show that this discretization preserves the qualitative behavior of the initial boundary value problem studied in the last section. More precisely we establish discrete versions of the inequalities (22) and (28) under onditions (20), (21), (23), (25) or (26), (27), (29), (31), respectively. We assume that n ¼ 2; X is the square ½0; L  ½0; L and H ¼ ðh1 ; h2 Þ with hi > 0; i ¼ 1; 2. In X we introduce the spatial grid XH ¼ fðx1;i ; x2;j Þ; i ¼ 0; . . . ; Nh1 ; j ¼ 0; . . . ; Nh2 g, where x‘;i ¼ x‘;i1 þ h‘ ; i ¼ 1; . . . ; N h‘ ; x‘;0 ¼ 0; x‘;Nh‘ ¼ L, for ‘ ¼ 1; 2. By @ XH we represent the set of boundary points. We introduce the following auxiliary points x‘;1 ¼ x‘;0  h‘ ; x‘;Nh‘ þ1 ¼ x‘;Nh‘ þ h‘ ; ‘ ¼ 1; 2. To

simplify

the

presentation we use the notation

@u , where a is a scalar funcwi;j ¼ wH ðx1;i ; x2;j Þ. We discretize @x@1 a @x 1 tion using the usual second order finite difference discretization

r

 ^ h1 ðaH

 1 rh1 uH Þðx1;i ; x2;j Þ ¼ aiþ1=2;j Dx1 uiþ1;j  ai1=2;j Dx1 ui;j ; h1 ð32Þ

h2

1;i

2;j

@x2

@x2

defined analogously. The semi-discrete approximations for u and v in XH at time t; uh ðtÞ and v H ðtÞ, are defined by the following system of ordinary differential equations

8 X  X Rt 0 > rhi ðdi rhi uH ðtÞÞ þ 0 ker ðt  sÞ rhi ðdv ;i rhi uH ðsÞÞ ds > > uH ðtÞ ¼ > > i¼1;2 i¼1;2 < > > > > > :

ðb1 þ GðtÞÞuH ðtÞ þ b2 v H ðtÞ in XH ;

v 0H ðtÞ ¼ ðq  b2  GðtÞÞv H ðtÞ þ b1 uH ðtÞ in XH ; ð33Þ

complemented with the initial conditions

uH ð0Þ ¼ RH u0 ;

v H ð0Þ ¼ RH v 0

in XH ;

ðwH ; qH Þ@XH ;x1 ¼

Dgx ui;j ðtÞ ¼ 0;

i ¼ 0 . . . ; N h1 ;

ð34Þ

0

1

h2

ðwH ; qH Þh1

0 1 Nh Nh 1 1 2 X X X h2 ¼ h1 @ h2 wi;j qi;j þ wi;j qi;j A; 2 j¼1 i¼1 j¼0;N h2

being ðwH ; qH Þh2 defined analogously, kwH k2H ¼ ðwH ; wH ÞXH ; P 2 kwH k1;s ¼ i¼1;2 kDxi wH k2hi . The following identity has a central role in what follows

rh‘ ða‘ rh‘ wH Þ; wH

XH

^‘;H Dx‘ wH ;Dx‘ wH Þh þ ðDa‘ ;gx wH gx ; wH Þ ¼ ða ‘ ‘

j ¼ 0; Nh2 ;

;

In what follows we establish an upper bound for the semi-discrete version of M 2 ðtÞ

M2;H ðtÞ ¼ kuH ðtÞk2H þ kv H ðtÞk2H ; where uH ðtÞ and v H ðtÞ are defined by (33)–(35). Let EH ðtÞ be the semi-discrete version of EðtÞ

EH ðtÞ ¼ M2;H ðtÞ þ

2 Z qffiffiffiffiffiffiffiffiffiffi X t ^ ker ðt  sÞ d : u ðsÞ ds D v ;i;H xi H 0

i¼1;2

hi

Multiplying both equations of (33) by uH ðtÞ and v H ðtÞ, respectively, with respect to the inner product ð:; :ÞXH and taking into account (38) we deduce

X d ^ Dx uH ðtÞ; Dx uH ðtÞÞ  ðb þ GðtÞÞkuH ðtÞk2 kuH ðtÞk2H ¼  ðd i;H 1 H i i hi dt i¼1;2 Z t þ b2 ðv H ðtÞ; uH ðtÞÞXH  ker ðt  sÞ 0 X ^ ðd  v ;i;H Dxi uH ðsÞ; Dxi uH ðtÞÞh ds X

i

ðDgx uH ðtÞgxi ; uH ðtÞÞ

ð39Þ

d kv H ðtÞk2H ¼ ðq  b2  GðtÞÞkv H ðtÞk2H þ b1 ðuH ðtÞ; v H ðtÞÞXH : dt

ð40Þ

i

@ X;xi

ð35Þ and

t

ker ðt  sÞDv ;gxm ui;j ðsÞ ds;

@ XH ;x‘



ð38Þ

i¼1;2

ð36Þ

In (34) RH denotes the restriction operator and in (36) Da;gx is 1 defined by

Da;gx wi;j ¼

j¼1

1 X 1 h2 wi;j qi;j þ h2 wi;j qi;j A; 2 j¼0;N

for grid functions defined on @ XH , being ðwH ; qH Þ@ XH ;x2 defined analogously, and for wH ; qH 2 W H ðXH Þ

j ¼ 0; . . . ; N h2 ;

where, for m ¼ 1; 2,

Dgxm ui;j ðtÞ ¼ Dd;gxm ui;j ðtÞ þ

@ 1

þ

i ¼ 0; Nh1 ;

Z

2 X

i¼1;2

Dgx ui;j ðtÞ ¼ 0; 2

N h 1

i¼0;N h

and the boundary conditions 1

0

X

‘ ¼ 1;2:



where ai1=2;j ¼ a x1;i  h21 ; x2;j and Dx1 denotes the usual backward finite difference operator in the x1 direction. The second order

^H r uH Þðx ; x Þ of @ b @u is finite difference discretization r ðb h2

where xi;j ¼ h1 h2 in XH , xi;j ¼ 12 h1 h2 on @ XH  CH , xi;j ¼ 14 h1 h2 on CH , CH denotes the set of corner points of X and @ XH ¼ XH \ @ X. The norm induced by the inner product (37) is denoted by k:kH . To simplify the presentation we use the following notations:

 1 aiþ1=2;j Dx1 wiþ1;j þ ai1=2;j Dx1 wi;j ; 2

being Db;gx wi;j defined analogously. 2 To prove a discrete version of the upper bounds (22), (28) we follow [19] to introduce a convenient discrete functional context. By W H ðXH Þ we denote the space of grid functions defined in XH . In W H ðXH Þ we introduce the inner product

Considering the boundary conditions (35), the discrete version of (15)

Z t 2 qffiffiffiffiffiffiffiffiffiffi d ^ ker ðt  sÞ d u ðsÞ ds D x H v ;i;H i dt 0 hi Z t  ^ ¼2 ker ðt  sÞd v ;i;H Dxi uH ðsÞ ds; Dxi uH ðtÞ 0

Z t 2 qffiffiffiffiffiffiffiffiffiffi ^  2b ker ðt  sÞ dv ;i;H Dxi uH ðsÞ ds ;

hi

0

and the Cauchy–Schwarz inequality, we get from (39), (40)

88

J.R. Branco et al. / Mathematical Biosciences 255 (2014) 83–90

5. Results

Migration cells u, day 6

15 k=0

k=0.45

k=0.225

k=0.065

In this section we present some numerical results illustrating the behavior of the glioma cells concentrations defined by (33)–(35). The numerical results were obtained integrating in time the ordinary differential problem using the implicit Euler method and discretizing the integral term in (33) with a left rectangular rule. We consider a homogeneous square domain X ¼ ½0; 15 cm ½0; 15 cm, a growth rate q ¼ 0:012=day and switching parameters

10

1

3

1 0

b1 ¼ 106 =day, and b2 ¼ 0:036=day. These parameters were obtained from [39] and have biological meaning. According to

2

0

[32], the initial condition is defined by 105 cells=cm2 of proliferation tumor cells at middle point of the domain, E0 ¼ 3156 Pa, E1 ¼ 6E0

5

and

l ¼ 8:9  104 Pa s. We also consider k ¼ 1 cm2 , isotropic

~11 ¼ d ~22 ¼ 0:004 cm2 =day2 . To guarantee the positivbehavior with d ~ ~v ;22 ¼ 1014 cell=Pa day [7,30] which leads ity of Dv we take dv ;11 ¼ d 0 0

5

10

2

15

Fig. 2. Numerical results at day 6 (logarithmic scale).

Table 1 kmin defined by (21) and (23), for a protocol of 24 consecutive hours of chemotherapy. Protocol

kmin (./day)

Each 2 days Each 7 days Each 14 days

0.064 0.224 0.448

E0H ðtÞ 6 2 max

  b2  b1 b  b2  GðtÞ; 1 þ q  GðtÞ; b EH ðtÞ; 2 2

ð41Þ

If (20) and (21) hold then 2

 b

M 2;H ðtÞ 6 e

 Rt

1 b2 þq 2

t

0

GðsÞ ds

M 2;H ð0Þ

ð42Þ

Otherwise, if (26) and (27) hold then 2



M 2;H ðtÞ 6 e

R

t b2 b1 t 2 0

GðsÞ ds

M2;H ð0Þ;

ð43Þ

Finally, if (20), (21), (23) and (25) hold or (26), (27), (29) and (31), hold then the discrete artificial mass M 2;H ðtÞ is bounded by M 2;H ð0Þ. 

 Rt b1 b2 2 þq t GðsÞ ds 2 0 and In the first case the decreasing factor is e

Rt b b 2 2 2 1 t GðsÞ ds 0 e in the second case.

Migration cells u, day 104

15 k=0

k=0.45

0

k=0.225

10

10

k=0.065

10

0

1

k=0.225

0

2

3

0

k=0.45

4

0

k=0

1 0

2

Proliferation cells v, day 104

15 k=0.065

0

12

to d11 ¼ d22 0:004 cm2 =day and dv ;11 ¼ dv ;22 ¼ 0:001 cm2 =day . In Fig. 2 we plot the numerical solution at day 6. Solution is presented in logarithmic scale, which means that contour plots represent the power of 10 of the density of tumor cells. In this case we observe a very intense spreading of proliferation cells and we can conclude that migration cells are already quite far from the core. Let us consider now that the chemotherapy treatment defined by (5) is applied with a protocol as illustrated in Fig. 1. Conditions (21) and (23) are used to compute a profile for GðtÞ that lead to control the total tumor mass. We consider a 24 h dosage and different rest periods. In Table 1 we show the minimum value of k ((5)) defined by conditions (21) and (23), for a virtual patient as defined in the beginning of this section. In Fig. 3 we plot the cell distribution at day 104 for an untreated patient and three patients with different chemotherapy protocols. All three protocols start at day 7, and follow different rest periods represented in Fig. 1 (k ¼ 0:065=day is associated to a 24 h dosage at days 7, 9, 11, 13, etc.; k ¼ 0:225=day is associated to a 24 h dosage at days 7, 14, 21, 28, etc.; k ¼ 0:45=day is associated to a 24 h dosage at days 7, 21, 35, 49, etc.). We observe that glioma mass density is significantly reduced when chemotherapy is used, although we do not observe a significant reduction of the tumor’s area. We remind that all plots are presented using logarithmic scale of the density of cells. Finally, in Fig. 4 we compare glioma masses of the virtual patient when no chemotherapy is administered and the results of the adminstration of the above three chemotherapy protocols. We observe a significant reduction of glioma masses when

3

2 1

0

5

0

12

5 1

0

0

0

5

10

15

0

0

5

Fig. 3. Numerical results at day 104 (logarithmic scale).

10

15

J.R. Branco et al. / Mathematical Biosciences 255 (2014) 83–90

Regional Development Fund through the program COMPETE and by the Portuguese Government through the FCT – Fundação para a Ciência e Tecnologia under the Project PEst-C/MAT/UI0324/2013.

6 k=0

k=0.45

k=0.225

89

k=0.065

tumor mass [logaritmic scale]

5

References 4

3

2

1

0

0

50

100

150

200

250

300

350

time [days] Fig. 4. Glioma masses M 1 ðtÞ for 1 year (logarithmic scale).

compared to glioma’s untreated patient. The results presented in this figure shows the effectiveness of the our approach to define chemotherapy protocols. 6. Conclusions In this paper we studied a mathematical model to describe the evolution of glioma cells with and without chemotherapy. The model was established combining a mass conservation law with a non Fickian mass flux that takes into account the viscoelastic behavior of the brain tissue described by the Voigt–Kelvin model. Using the energy method we deduced estimates for a functional related glioma mass M 2 ðtÞ defined using L2 norm. These estimates allowed us to define sufficient conditions on the parameters that lead to the control of M 2 ðtÞ. Semi-discrete approximations for the proliferation and migratory cancer cells presenting the same qualitative behavior of the continuous counterparts were introduced. Sufficient conditions that relate the chemotherapy effect with the growing rates of the semi-discrete migratory and proliferations glioma cells that lead to M 2;H ðtÞ < M 2;H ð0Þ were also established. These conditions allow us to define efficient protocols that lead to a decreasing of the cancer mass. Numerical experiments illustrating the behavior of the glioma mass under the conditions deduced for the chemotherapy protocols are also included. The results obtained suggest that our approach is a promising one. Models that will take into account the space effect of chemotherapy will be addressed in a future work. In this case we need to incorporate in the diffusion equation for the drug concentration the viscoelastic effect of the brain tissue on the diffusion drug similar to the one considered here in the migration of glioma cells. As the complete model is composed by integro-differential quasilinear equations of diffusion–reaction type, its mathematical and numerical study is a challenging problem. Future work will also address mathematical models with a real geometry and that take into account the white and gray matter of the brain. Comparison of the model with existing medical protocols will be also considered. Acknowledgements This work was partially supported by the Centro de Matemática da Universidade de Coimbra (CMUC), funded by the European

[1] E. Abad, S. Yuste, K. Lindenber, Reaction-subdiffusion and reactionsuperdiffusion equations for evanescent particles performing continuous time random walks, Phys. Rev. 81 (031115) (2010) 1–7. [2] S. Astanin, L. Preziozi, Multiphase models of tumor cells, in: N. Bellomo, M. Chaplain, E. De Angelis (Eds.), Selected Topics in Cancer Modeling: Genesis, Evolution, Immune Competition, and Therapy, Birkauser, 2008, pp. 223–253. [3] H. Byrne, J.R. King, D. McElwain, A two-phase model of solid tumor growth, Appl. Math. Lett. 16 (2003) 567–573. [4] J.R. Branco, J.A. Ferreira, P. de Oliveira, Memory in mathematical modelling of highly diffusive tumors, in: J. Vigo Aguiar (Ed.), Proceedings of the 2011 International Conference on Computational and Mathematical Methods in Science and Engineering, Benidorm, Spain, June 26–30, 2011, pp. 242–253. [5] J.R. Branco, J.A. Ferreira, P. de Oliveira, Efficient protocols to control glioma growth, in: I.P. Hamilton, J. Vigo-Aguiar (Eds.), Proceedings of the 2013 International Conference on Computational and Mathematical Methods in Science and Engineering, Almeria, Spain, June 24–27, 2013, pp. 293–303. [6] J.R. Branco, J.A. Ferreira, P. Oliveira, A viscoelastic model for glioma growth, in: Tuncer Ãren, Janusz Kacprzyk, Leifur Leifsson, Mohammad S. Obaidat, Slawomir Koziel (Eds.), Proceedings of SIMULTECH 2013 – 3rd International Conference on Simulation and Modeling, Technologies and Applications, Methodologies, Reykjavik, Iceland, 29–31 July, 2013, pp. 689–695. [7] D. De Kee, Q. Liu, J. Hinestroza, Viscoelastic (non-Fickian) diffusion, Can. J. Chem. Eng. 83 (2005) 913–929. [8] D.A. Edwards, D.S. Cohen, An unusual moving boundary condition arising in anomalous diffusion problems, SIAM J. Appl. Math. 55 (1995) 662–676. [9] D.A. Edwards, D.S. Cohen, A mathematical model for a dissolving polymer, AIChE J. 41 (1995) 2345–2355. [10] D.A. Edwards, Non-Fickian diffusion in thin polymer films, J. Polym. Sci., Part B: Polym. Phys. Ed. 34 (1996) 981–997. [11] D.A. Edwards, A spatially nonlocal model for polymer–penetrant diffusion, J. Appl. Math. Phys. 52 (2001) 254–288. [12] A. Engler, S. Sen, H. Sweeney, D. Discher, Matrix elasticity directs stem cell lineage specification, Cell 126 (2006) 677–689. [13] A. Engler, H. Sweeney, D. Discher, J. Schwarzbauer, Extracellular matrix elasticity directs stem cell differentiation, J. Musculoskelet Neuronal Interact. 7 (2007) 335. [14] S. Fedotov, Traveling waves in a reaction–diffusion system: diffusion with finite velocity and Kolmogorov–Petrovskii–Piskunov kinetics, Phys. Rev. E 58 (4) (1998) 5143–5145. [15] S. Fedotov, Nonuniform reaction rate distribution for the generalized Fisher equation: ignition ahead of the reaction front, Phys. Rev. E 60 (4) (1999) 4958– 4961. [16] S. Fedotov, A. Iomin, Migration and proliferation dichotomy in tumor-cell invasion, Phys. Rev. Lett. 98 (2007) 118110(1)–118110(4). [17] S. Fedotov, A. Iomin, Probabilistic approach to a proliferation and migration dichotomy in tumor cell invasion, Phys. Rev. E 77 (2008) 1031911(1)– 1031911(10). [18] S. Fedotov, Non-Markovian random walks and nonlinear reactions: subdiffusion and propagation fronts, Phys. Rev. E 81 (011117) (2010) 1–7. [19] J.A. Fereira, P. de Oliveira, P. Silva, Reaction–diffusion in viscoelastic materials, J. Comput. Appl. Math. 236 (2012) 3783–3795. [20] J.A. Ferreira, M. Grassi, E. Gudiño, P. de Oliveira, A 3D model for mechanistic control of drug release, SIAM J. Appl. Math. 74 (2014) 620–633. [21] G. Franceschini, The mechanics of human brain tissue, PhD Thesis, University of Trento, 2006. [22] A. Giese, L. Kluwe, B. Laube, H. Meissner, M. Berens, M. Westphal, Migration of human glioma cells on myelin, Neurosurgery 38 (1996) 755–764. [23] S. Habib, C. Molina-París, T. Deisboeck, Complex dynamics of tumors: modeling an emerging brain tumor system with coupled reaction–diffusion equations, Physica A 327 (2003) 501–524. [24] S. Hassanizadeh, On the transient non-Fickian dispersion theory, Transp. Porous Media 23 (1996) 107–124. [25] H.L. Harpold, E.C. Alvord Jr, K.R. Swanson, The evolution of mathematical modeling of glioma proliferation and invasion, J. Neuropathol. Exp. Neurol. 66 (2007) 1–9. [26] P. Hinow, P. Gerlee, L. McCawley, V. Quaranta, M. Cionanu, S. Wang, J. Graham, B. Ayati, J. Claridge, K. Swanson, M. Loveless, A. Anderson, A spatial model of tumor host interaction: application of chemotherapy, Math. Biosci. Eng. 6 (2009) 521–545. [27] J. Humphrey, Continuum biomechanics of soft biological tissues, Proc. R. Soc. Lond. 459 (2003) 3–46. [28] D. Joseph, L. Preziosi, Heat waves, Rev. Mod. Phys. 61 (1989) 47–71. [29] N. Laperrierea, L. Zurawb, G. Cairncross, The Cancer Care Ontario Practice Guidelines Initiative Neuro-Oncology Disease Site Group, Radiotherapy for newly diagnosed malignant glioma in adults: a systematic review, Radiother. Oncol. 64 (2002) 259–273. [30] Q. Liu, X. Wang, D. De Kee, Mass transport through swelling membranes, Int. J. Eng. Sci. 43 (2005) 1464–1470.

90

J.R. Branco et al. / Mathematical Biosciences 255 (2014) 83–90

[31] C. Lo, H. Wang, M. Dembo, Y. Wang, Cell movement is guided by the rigidity of the substrate, Biophys. J. 79 (2000) 144–152. [32] A. Mehrabian, Y. Abousleiman, General solutions to poroviscoelastic model of hydrocephalic human brain tissue, J. Theor. Biol. 291 (2011) 105–118. [33] J.D. Murray, Mathematical Biology, Springer, 2002. [34] S.P. Neuman, D.M. Tartakovsky, Perspective on theories of anomalous transport in heterogeneous media, Adv. Water Resour. 32 (2009) 670–680. [35] E. Norris, J. King, M. Byrne, Modelling the response of spatial structured tumours to chemotherapy: drug kinetics, Math. Comput. Modell. 43 (2006) 820–837. [36] L. Preziozi, A. Tosin, Multiphase modelling of tumor growth and extracellular matrix interaction: mathematical tools and applications, J. Math. Biol. 58 (2009) 625–656. [37] L. Preziozi, G. Vitale, A multiphase model of tumour and tissue growth including cell adhesion and plastic re-organisation, Math. Models Methods Appl. Sci. 21 (2011) 1901–1932. [38] R. Rockne, E.C. Alvord Jr, J.K. Rockhill, K.R. Swanson, A mathematical model for brain tumor response to radiation therapy, J. Math. Biol. 58 (2009) 561–578. [39] A. Roniotis, K. Marias, V. Sakkalis, M. Zervakis, Diffusive modelling of glioma evolution: a review, J. Biomed. Sci. Eng. 3 (2010) 501–508. [40] G. Sciumè, S. Shelton, W. Gray, C. Miller, F. Hussain, M. Ferrari, P. Decuzzi, B. Schrefler, A multiphase model for three-dimensional tumor growth, New J. Phys. 15 (2013) 015005.

[41] S. Shaw, J.R. Whiteman, Some partial differential Volterra equation problems arising in viscoelasticity, in: Proceedings of the Conference on Differential Equations and their Applications, Brno, 1997, pp. 183–200. [42] I. Sokolov, M. Schmidt, F. Sagués, On reaction-subdiffusion equations, Phys. Rev. E 73 (2006) 031102. [43] K.R. Swanson, E.C. Alvord Jr, J.D. Murray, A quantitative model for differential motility of gliomas in grey and white matter, Cell Prolif. 33 (2000) 317–329. [44] K.R. Swanson, C. Bridge, J.D. Murray, E.C. Alvord Jr, Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion, J. Neurol. Sci. 216 (2003) 1–10. [45] R. Tilghman, C. Cowan, J. Mih, Y. Koryakina, D. Gioli, J. Slack-Davis, B. Blackman, D. Tschumperlin, J. Parsons, Matrix rigidity regulates cancer cell growth and cellular phenotype, PloS One 5 (2010) e12905. [46] P. Tracqui, G.C. Cruywagen, D.E. Woodward, G.T. Bartoo, J.D. Murray, E.C. Alvord Jr, A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth, Cell Prolif. 28 (1995) 17–31. [47] T.A. Ulrich, E.M. de Juan Pardo, S. Kumar, The mechanical rigidity of the extracellular matrix regulates the structure, motility, and proliferation of glioma cells, Cancer Res. 69 (2009) 4167–4174. [48] J. Ward, J. King, Mathematical modelling of drug transport in tumour multicell spheroids and monolayer cultures, Math. Biosci. 181 (2003) 177–207.

Mathematical modeling of efficient protocols to control glioma growth.

In this paper we propose a mathematical model to describe the evolution of glioma cells taking into account the viscoelastic properties of brain tissu...
487KB Sizes 3 Downloads 5 Views