Mathematical Biosciences 248 (2014) 88–96
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
On assessing quality of therapy in non-linear distributed mathematical models for brain tumor growth dynamics A.S. Bratus b, E. Fimmel a,⇑, S.Yu. Kovalenko c a
Mannheim University of Applied Sciences, Paul-Wittsack-Str. 10, 68163 Mannheim, Germany Lomonosov Moscow State University, Faculty of Computational Mathematics and Cybernetics, GSP-1, 1/52, Leninskie Gory, 119991 Moscow, Russia c Federal Science and Clinical Center of the Federal Medical and Biological Agency, 28 Orehovuy boulevard, 115682 Moscow, Russia b
a r t i c l e
i n f o
Article history: Received 15 June 2012 Received in revised form 12 December 2013 Accepted 16 December 2013 Available online 31 December 2013 Keywords: Optimal therapy control Chemotherapy Cancer model
a b s t r a c t In this paper a mathematical model for glioma therapy based on the Gompertzian law of cell growth is presented. In the common case the model is considered with non-linear spatially varying diffusion depending on a parameter. The case of the linear spatially-varying diffusion arose as a special case for a particular value of the parameter. Effectiveness of the medicine is described in terms of a therapy function. At any given moment the amount of the applied chemotherapeutic agent is regulated by a control function with a bounded maximum. Additionally, the total quantity of chemotherapeutic agent which can be used during the treatment process is bounded. The main goal of the work is to compare the quality of the optimal strategy of treatment with the quality of another one, proposed by the authors and called the alternative strategy. As the criterion of the quality of the treatment, the amount of the cancer cells at the end of the therapy is chosen. The authors concentrate their efforts on finding a good estimate for the lower bound of the cost-function. Thus it becomes possible to compare the quality of the optimal treatment strategy with the quality of the alternative treatment strategy without explicitly finding the optimal control function. Ó 2013 Elsevier Inc. All rights reserved.
Introduction Mathematical models of brain tumor growth under consideration of diffusion have been developed and studied in many papers (see for example [31–33,16,34,3,25,27,35,17,20,21,24]). An optimal-control approach to cancer treatment, especially to brain cancer treatment, is well-established (see [8,9,11,18]). In particular, a comparison of optimal and suboptimal treatment protocols was already considered for some models of tumor anti-angiogenesis in [22]. The basic model we are developing further was suggested in [31,33] and is based on the model of Murray [26]. In these models the authors assume a linear spatially-varying diffusion and the exponential growth law for the tumor cells:
@cðx; tÞ ¼ rðDðxÞrcðx; tÞÞ þ qcðx; tÞ; @t
ðsee ½29; from p:536Þ
@cðx; tÞ ¼ rðDðxÞrcðx; tÞÞ þ qcðx; tÞ GðtÞcðx; tÞ; @t
ðsee ½36; 34Þ
⇑ Corresponding author. Tel.: +49 6212926243; fax: +49 6212926237. E-mail address: e.fi[email protected] (E. Fimmel). 0025-5564/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mbs.2013.12.007
where the function c denotes the density of tumor cells, DðxÞ has only two values DðxÞ ¼ DG , for x in gray matter and DðxÞ ¼ DW for x in white matter with DW > DG ; q 2 Rþ is a constant and GðtÞ ¼ k 2 R when the chemotherapy is being administered, and GðtÞ ¼ 0 otherwise. The model suggested in Section 1 is distinguished from the Murray et al. model mainly by the assumption of the Gompertzian law of cell growth and the non-linear spatially-varying diffusion depending on a parameter. The dynamics of the amount of the chemotherapeutic agent is in contrast to the model of Murray et al. a part of the model too. The class of the therapy functions considered is wider than in [31]. In the model suggested in this article we will adopt the Gompertzian growth law for the tumor cells since the Gompertzian models ensure that the growth rate of the tumor cells reduces as the maximal density of the tumor cells is approached while the exponential growth assumes unlimited growth of the density of tumor cells. Another significant difference between the model suggested and the model from [31,33] is a non-linear spatially-varying diffusion as in [1,29]. The case of the linear spatially-varying diffusion arose as a special case for a particular parameter value. The description of the diffusion of the malignant glioma cells is still a subject of discussion (see [1,2,13,7,14,30]). In [30] it is argued that
89
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
there is no significant difference in motility between cultures of a gross tumor and a histologically normal brain. On the other hand many authors stress an extensive non-linear character of the tumor cells distribution [12,1,13]. In [1] a non-linear diffusion in the form rðC k rCÞ depending on a parameter k is considered. For k > 0 the tissue is modeled rather as a porous medium and for k ¼ 0 as a special case with a linear diffusion studied in detail in [2]. Analogous law is usually used for the description of the diffusion in porous medium (see, for instance, [4,29]). We assume this description of the diffusion for the model suggested too. The nature of the interaction between the tumor size and the prescribed treatment is still up for debate. In the present work a wide class of increasing concave bounded therapy functions is examined. The main goal of the work is to compare the quality of the optimal strategy of treatment with the quality of another one, proposed by the authors in Section 4 and called ‘the alternative’. As the criterion of the quality of the treatment, the amount of the cancer cells at the end of the therapy is chosen. At any given moment the amount of the applied chemotherapeutic agent is regulated by a control function with a bounded maximum. Additionally, the total quantity of chemotherapeutic agent which can be used during the treatment process is bounded. Despite the essential progress achieved in the solving of optimization problems for distributed systems (see [23]) the search for the optimal control for problems described in terms of non-linear equations of parabolic type is still mathematically very difficult. Therefore the authors concentrate their efforts in Section 3 on finding a good estimate for the lower boundary of the cost-function. Thus it becomes possible to compare the quality of the optimal treatment strategy with the quality of the alternative treatment strategy without explicitly finding the optimal control function. The corridor between the analytical lower estimate and the value of the cost-function on the alternative control enables qualitative assessments and forecasts for the results of the treatment to be made.
notes the growth rate of tumor cells, clim the maximum density of the tumor cells, dh is the diffusion coefficient of the medicine, ch the dissipation rate of the medicine, GðhÞ describes the influence of the chemotherapeutic (called therapy function). We will consider the class of increasing concave bounded functions with Gð0Þ ¼ 0 e.g.,
1. Model and optimization problem formulation
(see Section 3 of the present article).
Let X R3 be a bounded domain of volume S with a smooth boundary C; n be the outer normal unit vector to C. For f 2 L1 ðXÞ we denote onward
The actual optimal control problem
f :¼
Z
GðhÞ ¼
kh ; 1þh
k 2 Rþ :
The control function uðx; tÞ 2 L1 ðX ½0; TÞ denotes the quantity of the chemotherapeutic agent applied to a patient at the moment t in x 2 X. For the cumulative amount of the chemotherapeutic agent the following inequality
Z
T
0
Z
hðx; tÞdxdt 6 Q 0
ð1:2Þ
X
with Q 0 > 0 holds. Note 1. For a ¼ 0 (respectively ah ¼ 0) the diffusions become linear
A0 cðx; tÞ ¼ rðDðxÞrcðx; tÞÞ;
A0 hðx; tÞ ¼ dh Dhðx; tÞ:
The model description is completed with the following initial conditions
cðx; 0Þ ¼ c0 ðxÞ > 0;
hðx; 0Þ ¼ 0
ð1:3Þ
and some boundary conditions which prohibit the migration of the tumor cells outside the brain boundary. It makes sense to consider two kinds of boundary conditions: either if a P 0
@cðx; tÞ j ¼ 0; @n Cð0;T
@hðx; tÞ j ¼ 0; @n Cð0;T
ð1:4Þ
or if a > 12
cðx; tÞjCð0;T ¼ 0;
@hðx; tÞ j ¼0 @n Cð0;T
ð1:5Þ
Let K # L1 ðX ½0; TÞ. In what follows we will consider three possibilities to define K:
f ðxÞdx:
X 1
Let cðx; tÞ 2 C ðX ½0; TÞ be the density of tumor cells and hðx; tÞ the amount of the chemotherapeutic in x 2 X at the moment t. We will consider the following model describing the growth of the tumor cells under the influence of chemotherapy:
– K 1 : fuðx; tÞ 2 L1 ðX ½0; TÞ j uðx; tÞ ¼ dðx x ÞuðtÞ; x 2 Xg where d be the Dirac function, – K 2 : fuðx; tÞ 2 L1 ðX ½0; TÞ j uðx; tÞ ¼ uðtÞg, – K 3 : fuðx; tÞ 2 L1 ðX ½0; TÞ j uðx; tÞ ¼ vC ðxÞuðtÞg, where
vC ðxÞ :¼
@cðx; tÞ ¼ qcðx; tÞð1 b ln cðx; tÞÞ þ Aa cðx; tÞ cðx; tÞGðhÞ @t @hðx; tÞ ¼ ch hðx; tÞ þ Aah hðx; tÞ þ uðx; tÞ; 0 6 t 6 T; @t
ð1:1Þ
where q; b; dh ; ch ; clim ; a; ah 2 Rþ are constants and b ¼ ln 1c ,
1;
x 2 C;
0;
x R C
is the characteristic function of a compact set C X of area SC . In each of the three cases the inequality
0 6 uðtÞ 6 q
ð1:6Þ
lim
Aa cðx; tÞ :¼
3 X @ @cðx; tÞ ; DðxÞ cðx; tÞ2a @xi @xi i¼1
Aa h hðx; tÞ :¼ dh
3 X @ @hðx; tÞ hðx; tÞ2ah @xi @xi i¼1
with DðxÞ ¼ Dw for white matter and DðxÞ ¼ Dg (see [26]) for gray matter. The parameters a; ah characterize the degree of the non-linearity of the respectively diffusions. The parameter q (see [26]) de-
must hold with q > 0. Consider the following optimal control problem. We need to find the time T 2 ð0; 1Þ and the control function uðx; tÞ 2 K such that uðx; tÞ; T provide the infimum to the cost-function:
Uðu; TÞ ¼
Z
ln cðx; TÞdx ¼: ln cðTÞ
X
Let
inf Uðu; TÞ ¼ Uðu ; TÞ: u2K
ð1:7Þ
90
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
In what follows (see Section 3) we will find a lower estimate U for Uðu ; TÞ. It is evident that we get an upper estimate U for Uðu ; TÞ if we substitute an arbitrary non-optimal control function u in Uðu; TÞ. Thus the inequality holds.
U 6 Uðu ; TÞ 6 U: In what follows (see Sections 4 and 5) we will suggest a suboptimal control function for each functional space K i , i ¼ 1; 2; 3, called the alternative control function so that the corresponding value of the cost-function is in some sense sufficiently close to the obtained lower estimate of the cost-function. If this holds, the results achieved can be useful in the medical prognosis.
The following Theorem helps us to understand the relationship between the important constant k (see Theorem 2) and the restrictive constant q (see (1.6)) for the dose of the drug for each one of three functional spaces considered K i ; i ¼ 1; 2; 3 : Theorem 1. Let X; h; C X be as above, G : Rþ ! R be an increasing concave bounded function. Then
GðhÞ 6 S k and (1) k :¼ GScq for uðx; tÞ 2 K 1 , h (2) k :¼ Gcq for uðx; tÞ 2 K 2 , h SC q (3) k :¼ G Sc for uðx; tÞ 2 K 3 .
2. Relations between parameters
h
Lemma 1. The following inequalities take place: hðtÞ 6 cq
(1) For uðx; tÞ 2 K 1 we have
Z
T
Z
0
hðx; tÞdxdt ¼
0
X
h
T Z
Z
Proof. For a concave integrable function f ðxÞ the Jensen inequality [15] has the form
and
t ech ðtsÞ uðsÞds dt 6 Q 0
f
0
Z
Z
T
0
Z
hðx; tÞdxdt ¼ S
hðtÞ 6 cSq
0
X
h
t
and
Z
T
0
Z
h G S
ech ðtsÞ uðsÞds dt 6 Q 0
0
hðx; tÞdxdt ¼ SC
Z
h
t
ch ðtsÞ
e
uðsÞds dt 6 Q 0
! ¼G
Z Z 1 1 1 hðxÞdx P GðhðxÞÞdx ¼ GðhÞ; S X S S X
3. A lower estimate for the cost-function
0
SC qT 2 6 : 2
Our aim in this section is to obtain a lower estimate for the objective function considered which does not depend on the choice of a control function.
Proof. Integrating the second equation from (1.1) over X and using Stokes’s formula we get
Theorem 2. Let Gðhðx; tÞÞ; c1 ðx; tÞ 2 L1 ðXÞ,
dhðtÞ ¼ ch hðtÞ þ uðtÞ dt
ln c 2 L2 ðXÞ for a ¼ 0: Also we assume for a > 0 that
with hð0Þ ¼ 0. Hence
hðtÞ ¼
Z
t
kðxÞdx X
S is the volume of the domain X. Since G is increasing and the estimate of h from Lemma 1 is valid, we obtain the upper estimates for GðhÞ. h
hðtÞ 6 ScC q and
T Z
0
X
X
Z
In our case
SqT 2 6 ; 2 (3) For uðx; tÞ 2 K 3 we have
kðxÞf ðhðxÞÞdx if kðxÞ P 0 and
¼ 1:
T Z
Z
Z P
X
qT 2 6 : 2 (2) For uðx; tÞ 2 K 2 we have
kðxÞhðxÞdx
ca ðx; tÞ 2 W 21 ðXÞ and
ech ðtsÞ uðsÞds 6
0
M
ch
ð1 ech t Þ 6
M
ch
if
uðtÞ 6 M ðÞ
2 Rþ :
Z
ln cðx; tÞdx < 1: X
Let k 2 Rþ 1 with GðhÞðtÞ 6 kS for all t 2 ½0; T and S being the volume of the domain X. Then the following estimates are valid:
Consequently we obtain: for uðx; tÞ 2 K 1 for uðx; tÞ 2 K 2 for uðx; tÞ 2 K 3
that that that
(1) for the boundary condition (1.4) with arbitrary a P 0
uðtÞ ¼ uðtÞ 6 q, uðtÞ ¼ S uðtÞ 6 Sq, uðtÞ ¼ SC uðtÞ 6 SC q.
dln c P Sðq kÞ qbln c; dt (2) for the boundary condition (1.5) and a > 12
This proves the first inequality in each case. Substituting in each RT R case uðtÞ in () we get the formulas to calculate 0 X hðx; tÞdxdt.
2 dln c Dg P Sðq kÞ qbln c þ ðln cÞ dt c1 Sa2
We get upper estimates for Q 0 in each case using the inequality
Z 0
T
hðtÞdt 6
Z 0
T
M
ch
ð1 ech t Þdt ¼
M
c
2 h
ð1 þ ech T þ T ch Þ 6
MT 2 2
if ch is sufficiently small and, in addition, Q 0 should be chosen to be as small as possible. h
where c1 is the constant from Friedrichs inequality which depends only on X.2 1 2
The number k can be viewed as a characteristic of treatment intensity. 2 2 For a rectangle X with sides a and b, for instance, c1 ¼ 24a2 b 2 [28, p. 193]. p ða þb Þ
91
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
Proof. Dividing the first equation of (1.1) by cðx; tÞ and integrating over X we get
Z
dln c ¼ qðS bln cÞ þ dt
X
Aa cðx; tÞ dx GðhÞ: cðx; tÞ
Aa cðx; tÞ 1 @ @cðx; tÞ dx dx ¼ DðxÞc2a ðx; tÞ cðx; tÞ cðx; tÞ @xi @xi i¼1 X 3 Z X @ @cðx; tÞ dx ¼ DðxÞc2a1 ðx; tÞ @xi @xi i¼1 X 3 Z X @ 1 @cðx; tÞ DðxÞc2a ðx; tÞ dx: @xi cðx; tÞ @xi i¼1 X
Z
3 Z X
X
In view of the Stokes’s theorem
@ @cðx; tÞ dx DðxÞc2a1 ðx; tÞ @xi @xi i¼1 X 3 Z X @cðx; tÞ DðxÞc2a1 ðx; tÞ cosðn; xi Þds: ¼ @xi i¼1 C
3 Z X
In the case of the boundary condition (1.4) we have 2 X i¼1
3 Aa cðx; tÞ Dg X dx P 2 cðx; tÞ a i¼1
@cðx; tÞ cosðn; xi Þ @xi
Cð0;T
@cðt; xÞ ¼ ¼0 @n Cð0;T
and in the case of the boundary condition (1.5) for a > 12 we have
c2a1 ðx; tÞjCð0;T ¼ 0:
Z
1 c1 a2
i¼1
X
@ @cðx; tÞ dx ¼ 0: DðxÞc2a1 ðx; tÞ @xi @xi
X
2 @cðx; tÞ dx @xi X a 2 3 Z 1X @c ðx; tÞ DðxÞ dx ¼ 2 @xi a i¼1 X
3 X Aa cðx; tÞ dx ¼ cðx; tÞ i¼1
Z
DðxÞc2a2 ðx; tÞ
and for a ¼ 0:
Z X
X
Aa cðx; tÞ dx ¼ cðx; tÞ
i¼1
Aa cðx; tÞ dx P 0: cðx; tÞ (1) Using GðhÞðtÞ 6 kS for all t 2 ½0; T we obtain the following inequality
dln c P Sðq kÞ qbln c dt which proves the estimate 1. (2) By virtue of the Friedrichs inequality [28]
Z X
u2 ðxÞdx 6 c1
1
Z Z
Sa2 c1
2 ca ðx; tÞdx X
2 ln cðx; tÞdx
X
where c1 depends only on X. Thus the estimate 2 is proven. h Theorem 3. For the boundary condition (1.4) with a P 0 and the condition (1.5) with a > 12 and q; S; k; clim as above, the following estimate holds
Uðu; TÞ P
Sðq kÞ ð1 eqbT Þ þ ln c0 eqbT : qb
Proof. For the boundary conditions (1.4) and (1.5) and q; S; k; clim as above the following estimate holds due to Theorem 2:
dln c P Sðq kÞ qln c: dt Let us consider the following initial condition problem
dz ¼ Sðq kÞ qbz; dt
zð0Þ ¼ ln c0 :
z¼
Sðq kÞ ð1 eqbt Þ þ ln c0 eqbt : qb
Using the comparison Lemma [19, p. 11] and Theorem 2 we get the result. h Note 2. In the case of the boundary condition (1.5) for a > 12 we can get a better lower estimate for Uðu; TÞ. In this case we deal with a Riccati equation
dz Dg 2 ¼ Sðq kÞ qbz þ z :¼ PðzÞ: dt c1 Sa2 D
2 @ ln cðx; tÞ DðxÞ dx: @xi X
3 Z X
In both cases we obtain
Z
X
1 Sc1 a2
The solution of this problem is
Thus we have for a – 0:
Z
c2a ðx; tÞdx P P
In both cases we get 3 Z X
2 Z a Z @c ðx; tÞ Dg dx P 2 c2a ðx; tÞdx: @xi a c1 X X
Moreover, in the view of the Cauchy–Bunyakovsky–Schwarz inequality we get
Aa cðx; tÞ dx: cðx; tÞ
X
Z X
Now let us estimate
Z
for an arbitrary function u 2 W 21 ðXÞ and some constants c1 ; c2 P 0 depending only on X, we obtain the following estimate in the case of the boundary condition (1.5) for a > 12:
2 Z 3 Z X @uðx; tÞ dx þ c2 u2 ðsÞds @xi C i¼1 X
Let us denote a0 :¼ Sðq kÞ; a1 :¼ qb; a2 :¼ c Sga2 . We obtain the 1 differential equation
dz ¼ a0 þ a1 z þ a2 z2 dt
with zð0Þ ¼ ln c0 :
Let us note first that only the situation
4a0 a2 þ a21 ¼ 4
Dg ðq kÞ þ q2 b 2 > 0 a2 c1
can occur in our case as we will see below in Section 4 due to the condition q < k. The roots of PðzÞ
z1;2 ¼
a1
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4a0 a2 þ a21 2a2
;
z1 < 0 < z2 ;
are particular solutions of the Riccati equation. As we see below, z1 is an attractor and z2 is a repeller. We set
y¼
1 : z z1
92
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
We get a linear differential equation
dy ¼ ða1 þ 2a2 z1 Þy a2 ; dt
yð0Þ ¼
1 ln c0 z1
:
The solution of this equation in the case a21 > 4a0 a2 (in particular, it means that a1 þ 2a2 z1 – 0) is given by
yðtÞ ¼
a2 a2 1 eða1 þ2a2 z1 Þt : þ þ a1 þ 2a2 z1 a1 þ 2a2 z1 ln c0 z1
Thus the solution of the equation is given by
1
þ 1 eða1 þ2a2 z1 Þt ln c0 z1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a21 4a0 a2 ¼ z1 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a2 4a a a2 þ a2 þ 1 0 2 e 4a0 a2 þa1 t
zðtÞ ¼ z1 þ
a2 þ a1 þ2a 2 z1
a2 a1 þ2a2 z1
consider the cases a ¼ 0; a ¼ 0:05 (for the Neumann boundary condition), and ah ¼ 0 and respectively ah ¼ 0:1. We need an estimate for the dissipation rate of the chemotherapeutic agent. Studying the pharmakinetics description of a drug with the most simple application schedule we note that the halfdecay time for men is 20 days. Thus we get ech 20 ¼ 12 or, respec1 tively, ch ¼ 0:0347 day . From the pharmakinetics description of the same drug we can also learn that a single dose administered in glioma case is 10 mg/kg. So we can estimate the parameter q in the case uðx; tÞ ¼ uðtÞ 2 K 2 as q2 ¼ 0:8. We assume for the control functions of all types that the product of the application area of the medicine and q is the same in all three cases. Thus we get
q2 S ¼ q3 SC ¼ q1 102 and finally q1 ¼ 2880; q3 ¼ 2:ð2Þ.
ln c0 z1
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 4a a þa2 with z1 ¼ 1 2a2 0 2 1 < 0. It follows for the boundary condition (1.5) and a > 12 in view of the comparison Lemma [19, p.11] that
Uðu; TÞ P zðTÞ: Evidently,
lim Uðu; TÞ P lim zðTÞ ¼ z1 :
T!1
T!1
However, the numerical calculations show that this more precise lower estimate is only marginally better than the estimate from Theorem 3.
Numerical Results A tumor can be usually diagnosed if its radius is equal to 15 mm. The maximum tumor extent radius is 37.5 mm ([20,6]). Thus we take as a simulation space the square with the side length a ¼ 6cm thinking as a 2-dimensional projection of the domain X to visualize the results in the 3-dimensional space. Further we represent the square as a 60 60-mesh consisting of equal squares having the area 1mm2 . As a time unit we take 1 day. We consider the following initial spatial density of the tumor cells ([26, p.547]):
( c0 ðx; yÞ :¼
ðxx0 Þ2 ðyy0 Þ2 B
Ae
;
1;
ðx x0 Þ2 þ ðy y0 Þ2 6 R20 else
3
4. Alternative control Formulation of the alternative control task
with A ¼ 2 10 ; R0 ¼ 1:8; B ¼ 0:6. Below we use for the numerical calculations ðx0 ; y0 Þ ¼ ð3; 3Þ: We consider for our numerical simulations the therapy function (see Fig. 1)
We will search alternative control functions in each functional space K i ; i ¼ 1; 2; 3 on the base of the following consideration: ~ ðx; tÞ 2 K i ; i ¼ 1; 2; 3, satisfying we look for such a function u (1.6), where the inequality (1.2) holds and for the functional
GðhÞ ¼
Fðu; TÞ ¼
Z
c0 ðxÞGðhðx; TÞÞdx
ð4:1Þ
X
kh ; k 2 Rþ 1þh
and calculate k in each of the three cases depending on a treatment intensity k is considered from
kS ¼ kmax
the equality
t2½0;T
Z X
h dx: 1þh
holds.
~ ðtÞ (see (4.1)) in the class We will search the proper functions u of the piecewise constant functions on the interval ½0; T dividing the whole interval into equal parts:
Parameter estimates
~ ðtÞ ¼ ui ; u
It will be useful in what follows to have some parameter values. We adopt the parameter estimates from [27,26,32]:
As a result of numerical calculations we obtain the following alternative control function in all three cases considered
~ ; TÞ ¼ sup Fðu; TÞ Fðu u
– Diffusion coefficient (Gray matter-White matter) 1
Dg ¼ 0:0013 cm2 day ;
~ ðx; tÞ ¼ u 1
Dw ¼ 0:005 cm2 day ;
– Maximal density of the tumor cells
– Growth rate of the tumor cells 1
q ¼ 0:012 day ; 1
In the case ah ¼ 0 we set dh ¼ 0:0864 cm2 day equal to the apparent diffusion coefficient (ADC) for the water. We will
ði 1Þ T i T ; ; n n
q;
0 6 t < t0 ;
0;
t0 6 t 6 T
i 2 f1; . . . ; ng;
ui 6 q:
where the parameter t 0 can be found from the equation
Z
clim ¼ 105 cells cm2 ) ln clim ¼ 11:5; b ¼ 0:0869
t2
t0
hðtÞdt ¼ Q 0 :
ð4:2Þ
0
Note that the moment of time t 0 and, as a consequence, the ~ ðx; tÞ depend on the parameter q alternative control function u and Q 0 as well as on the chosen time of the end of the therapy T. Let us also add that for the ODE systems in a similar case it is proven that the control function given above is the optimal control ([10,5]).
93
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
2000
c
1500
1000
500
0 6 5
6
4
5 3
4 3
2
2
1 0
y
1 0
x
Fig. 1. Initial conditions (distribution of cancer cells density) for the control task.
160 140
k=0 u=u(t)
∫ Ω ln c(x,T) dx
120 100 80 u=δ(x*)u(t)
60
u=χR’(x*)u(t) u=χ (x*)u(t) R’
40 20
0
10
20
30
40
50
60
70
80
90
T Fig. 2. The numerical upper and analytical lower (the lower black line) estimates for three control functions types in the case of the Neumann boundary conditions: for uðx; tÞ ¼ dðx x ÞuðtÞ (blue line), for uðx; tÞ ¼ vC ðxÞuðtÞ (green and yellow lines, respectively for different compact sets) and for uðx; tÞ ¼ uðtÞ (red line), a ¼ 0; ah ¼ 0; k ¼ 0:0196. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Let us compare the results shown in the pictures 2–5. The pictures 2–4 illustrate the case of the usual Laplace operator (a ¼ ah ¼ 0). The difference between them is in the increasing middle intensity of the therapy (k) from picture 2 to picture 4. The upper black line shows the growth of the cost function without any therapy administered for comparison aims in every picture. Surprisingly, we can see that while the increasing therapy intensity has little effect on the results in the cases uðx; tÞ 2 K 2 or uðx; tÞ 2 K 3 , in the case uðx; tÞ 2 K 1 it looks significantly better. Fig. 4 shows, in particular, a relative deceleration of the tumor reduction in the case of uðx; tÞ 2 K 1 . In this case the time moment t0 (compare 4.2) of the end of the active therapy for uðx; tÞ 2 K 1 occurs much earlier than the intended end of the therapy T due to the higher value than for uðx; tÞ 2 K 2 ; K 3 of the parameter q for this kind of control. The effect is more visible than in Fig. 2,3 due to the higher intensity of the application of the chemotherapeutic agent. The results obtained show, in particular, that the restriction on the cumulative amount as well as the intensity of the application of
the chemotherapeutic agent have a big influence on the therapy process which emphasises the need for further medical investigations in this area (see Fig. 5). The picture 5 shows the case of a small non-linearity in the diffusion (a ¼ 0:05; ah ¼ 0:1). In this case the results of the therapy in all cases are significantly worse than in the case of the linear diffusion because of a very extensive growth of the tumor. Let us also note that the lower estimate, calculated analytically because of its generality is relatively rough. The Fig. 6(a)–(c) shows numerically simulated change of the tumor after 90 days with the application of different therapy protocols. It is important to compare the scales on the axis c in Fig. (a)–(c). While the maximum value is still at about 1000 units after an application of the protocol 1, in the last picture we have the maximum value of only about 3 units. If we compare the initial volume of the tumor (with its maximum value at approximately 1500 units) and its volume after 90 days of therapy with three protocols considered, we can see that the
94
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
200 k=0
150
∫ ln c(x,T) dx
u=u(t) 100
50
Ω
u=δ(x*)u(t)
u=χR’(x*)u(t) u=χ (x*)u(t) R’
0
−50
0
10
20
30
40
50
60
70
80
90
T Fig. 3. The numerical upper and analytical lower estimates for three control functions types in the case of the Neumann boundary conditions: for uðx; tÞ ¼ dðx x ÞuðtÞ (blue line), for uðx; tÞ ¼ vC ðxÞuðtÞ (green and yellow lines, respectively for different compact sets) and for uðx; tÞ ¼ uðtÞ (red line), a ¼ 0; ah ¼ 0; k ¼ 2 0:0196. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
200 150
k=0
u=u(t)
∫ Ω ln c(x,T) dx
100 50 0
u=δ(x*)u(t)
u=χR’(x*)u(t) u=χ (x*)u(t) R’
−50 −100 −150
0
10
20
30
40
50
60
70
80
90
T Fig. 4. The numerical upper and analytical lower (the lower black line) estimates for three control functions types in the case of the Neumann boundary conditions: for uðx; tÞ ¼ dðx x ÞuðtÞ ((blue line)), for uðx; tÞ ¼ vC ðxÞuðtÞ (green and yellow lines, respectively for different compact sets) and for uðx; tÞ ¼ uðtÞ (red line), a ¼ 0; ah ¼ 0; k ¼ 3 0:0196. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
180 α=0.05; α =0.1
160
k=0
u(x,t)=u(t) u(x,t)=χR’(x*)u(t)
h
∫ Ω ln c(x,T) dx
140 120 100 u(x,t)=χΩ’(x*)u(t)
80 60
u(x,t)=δ(x*)u(t)
40 20
0
10
20
30
40
50
60
70
80
90
T Fig. 5. The numerical upper and analytical lower (the lower black line) estimates for three control functions types in the case of the Neumann boundary conditions: for uðx; tÞ ¼ dðx x ÞuðtÞ (continuous line), for uðx; tÞ ¼ vC ðxÞuðtÞ (dotted line) and for uðx; tÞ ¼ uðtÞ (dashed line), a ¼ 0:05; ah ¼ 0:1; k ¼ 0:0196.
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
(a)
1000
c
T
800 600 400 200 0 4 4
3 3
2
2
1
(b)
1
0 0
y
x
250 200 150 100 50 0 0 1 2 3 4 5
x
1
6 0
4
3
2
5
6
y
95
moment the amount of the chemotherapeutic agent applied is regulated by a control function with a bounded maximum. Additionally, the total quantity of chemotherapeutic agent which can be used during the treatment process is also bound. The main goal of the work is to compare the quality of the optimal strategy of treatment with the quality of some other treatment control strategies, proposed by the authors, called the alternative. As the criterion of the quality of the treatment (the cost-function), the amount of cancer cells at the end of the therapy is chosen.The alternative control strategy is presented for three types of control function. In the first case the dose of the chemotherapeutic agent varies depending only on time. In the second case the application of the chemotherapeutic agent is targeted on the tumor and the third one models the case of very precise and intensive administered application of the chemotherapeutic agent comparable with an injection. The results obtained show the preciser the delivery of the chemotherapeutical agent is targeted on the tumor the better the result of the therapy is. In the work the lower boundary of the cost-function is analytically estimated. Thus it becomes possible to compare the quality of the optimal treatment strategy with the quality of the alternative treatment strategy without solving the optimal control problem. Thus, in the work a corridor for the cost-function has been constructed which enables observation of the behavior of the tumor during the therapy. The lower boundary of the corridor is an analytical estimate. The upper boundary of this corridor is determined by the proposed suboptimal control functions. Since the proposed alternative control strategy is essentially easier to calculate as the optimal one, it enables qualitative assessments and forecasts for the results of the treatment to be made. Acknowledgments The authors thank Dr.Yu.S. Semenov (MIIT, Moscow, Russia), whose expertise provided valuable quality control. The authors are also grateful to Dr. G. Schmalz (UNE, Armidale, Australia) for his valuable and constructive comments.
(c) c(x,y,T=90)
3 2.5
References
2 1.5 6
1 0
4 2
2
4
x
6 0
y
Fig. 6. The graphs for cðx; TÞ with T ¼ 90 days and ah ¼ 0; a ¼ 0. (a) u ¼ uðtÞ ðq ¼ 0:8Þ (b) u ¼ vC ðx; x ÞuðtÞ ðq ¼ 2:8294Þ (c) u ¼ dðx ÞuðtÞ ðq ¼ 2880Þ.
biggest positive change by far is obtained after an application of the protocol 3 (Fig. 6(c)).
5. Conclusion Presented in the work is a mathematical spatial model for glioma therapy consisting of two non-linear partial differential equations of parabolic type with non-linear spatially varying diffusion depending on a parameter. The case of the linear spatially-varying diffusion arose as a special case for a particular value of the parameter. The first equation describes the dynamics of cancer cell growth under Gompertzian law and the interaction between a controlled quantity of the chemotherapeutic agent and cells given in terms of a therapy function. The second differential equation describes the dynamics of the chemotherapeutic agent. At any given
[1] E.K. Afenya, C.P. Calderón, Growth kinetic of cancer cells prior to detection and treatment an alternative view, Discrete Contin. Dyn. Syst. – Ser. B 4 (1) (2004). [2] E.K. Afenya, C.P. Calderón, Diverse ideas on the growth kinetics of disseminated cancer cells, Bull. Math. Biol. (2000). [3] M. Assefa, R. Rockne, M. Szeto, K.R. Swanson, Mathematical modeling of Glioma proliferation and diffusion, Ethnicity and Disease 19 (2, Suppl. 3) (2009) 60–61. Pubmed ID: 19554787. [4] G.I. Barenblatt, On some unsteady fluid and gas motions in a porous medium, Prikladnaya Matematika i Mekhanika (Applied Mathematics and Mechanics (PMM)) 16 (1) (1952) 67 (in Russian). [5] A.S. Bratus, E.S. Chumerina, Optimal control in therapy of solid tumor growth, Comput. Math. Math. Phys. 48 (6) (2008) 892–911. [6] P.K. Burgess, P.M. Kulesa, J.D. Murray, E.C. Alvord Jr., The interaction of growth rates and diffusion coefficients in a three-dimensional mathematical model of Gliomas, J. Neuropathol. Exp. Neurol. 56 (6) (1997) 704. [7] A.C. Burton, Rate of growth of solid tumor as a problem of diffusion, Growth 30 (1966) 157. [8] S.P. Chakrabarty, F.B. Hanson, Optimal control of drug delivery to brain tumors for a distributed parameters model, in: Proceedings of American Control Conference, 2005, p. 1. [9] M. El-Borai, Application of stochastic control theory in modeling of brain cancer, in: Casablanca International Workshop on Mathematical Biology: Analysis & Control, June 20–24, 2011. [10] E. Fimmel, Y.S. Semenov, A.S. Bratus, On optimal and suboptimal treatment strategies for a mathematical model of Leukemia, Math. Biosci. Eng. 10 (1) (2013) 151. [11] K.R. Fister, J.C. Panetta, Optimal control applied to competing chemotherapeutic cell-kill strategies, SIAM J. Appl. Math. 63 (6) (2003) 1954. [12] A. Friedman, Therapeutic approach to brain cancer, in: Abstracts of 8th European Conference on Mathematical and Theoretical Biology, Krakov, June 28–July 2, 2011, p. 318. [13] R.A. Gatenby, E.T. Gawliski, The glycolytic phenotype in carcinogenesis and tumor invasion: insight thought mathematical models, Cancer Res. 63 (15) (2003) 3847.
96
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
[14] V. Greenspan, Models for the growth of solid tumor as a problem by diffusion, Stud. Appl. Math. 52 (1972) 317. [15] G.H. Hardy, J.E. Littlewood, G. Polya, Inequalities, Cambridge University Press, 1934. [16] H.L.P. Harpold, E.C. Alvord Jr., K.R. Swanson, The evolution of mathematical modeling of glioma growth and invasion, J. Neuropathol. Exp. Neurol. 66 (1) (2007) 1. [17] T. Hines, Mathematically modeling the mass-effect of invasive brain tumors, SIURO 3 (2010). Available from: . [18] M. Itika, M.U. Salamcib, S.P. Banksa, Optimal control of drug therapy in cancer treatment, Nonlinear Anal. Theory Methods Appl. 71 (12) (2009) e1473. [19] E. Kamke, Differentialgleichungen Lösungsmethoden und Lösungen, Band 1, Gewöhnliche Differentialgleichungen, Chelsea Publishing Company, New York, 1971. [20] A.R. Kansal, S. Torquato, G.R. Harsh IV, E.A. Chiocca, T.S. Deisboeck, Simulated brain tumor growth dynamics using a three-dimensional cellular automaton, J. Theor. Biol. 203 (2000) 367. [21] E. Khain, L.M. Sander, A.M. Stein, A model for glioma growth, Complexity 11 (2005) 53. [22] U. Ledzewicz, H. Schaettler, Optimal and suboptimal protocols for a class of mathematical models for tumor anti-angiogenesis, J. Theor. Biol. 252 (2008) 295. [23] X. Li, Y. Yong, Optimal control Theory for Infinite Dimensional Systems, Birkhauser, Boston, 1995. [24] D. Mackenzie, Mathematical modeling and cancer, SIAM News 37 (1) (2004). [25] A. Martínez-González, G.F. Calvo, L.A. Pérez Romasanta, V.M. Pérez-García, Hypoxic cell waves around Necrotic cores in Glioblastoma: a biomathematical model and its therapeutic implications, Bull. Math. Biol. 62 (3) (2012) 527.
[26] Murray, Mathematical biology, vol. 2, Spatial Models and Biomedical Applications, 3rd ed., Springer, 2003. [27] G. Powathil, M. Kohandel, S. Sivaloganathan, A. Oza, M. Milosevic, Mathematical modeling of brain tumors: effects of radiotherapy and chemotherapy, Phys. Med. Biol. 52 (2007) 3291. [28] K. Rektoris, Variational Methods in Mathematics: Science and Engineering, Mir, Moscow, 1985 (in Russian). [29] A.A. Samarski, V.A. Galaktionov, S.P. Kurdyumov, A.P. Mikhailov, Blow-up in quasilinear parabolic equations, de Gruyter Expositions in Mathematics, vol. 19, de Gruyter, Berlin and Hawthorne, NY, 1995. [30] D.L. Silbergeld, M.R. Chicoine, Isolation and characterization of human malignant glioma cells from histologically normal brain, J. Neuroseurg. 86 (1997) 525. [31] 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. [32] 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. [33] K.R. Swanson, E.C. Alvord Jr., J.D. Murray, Quantifying efficacy of chemotherapy of brain tumors (Gliomas) with homogeneous and heterogeneous drug delivery, Acta Biotheor. 50 (4) (2002) 223. [34] K.R. Swanson, H.L.P. Harpold, D.L. Peacock, R. Rockne, C. Pennington, L. Kilbride, R. Grant, J. Wardlaw, E.C. Alvord Jr., Velocity of radial expansion of contrast-enhancing Gliomas and effectiveness of radiotherapy in individual patients: a proof of principle, Clin. Oncol. 20 (2008) 301. [35] 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.
Contents lists available at ScienceDirect
Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs
On assessing quality of therapy in non-linear distributed mathematical models for brain tumor growth dynamics A.S. Bratus b, E. Fimmel a,⇑, S.Yu. Kovalenko c a
Mannheim University of Applied Sciences, Paul-Wittsack-Str. 10, 68163 Mannheim, Germany Lomonosov Moscow State University, Faculty of Computational Mathematics and Cybernetics, GSP-1, 1/52, Leninskie Gory, 119991 Moscow, Russia c Federal Science and Clinical Center of the Federal Medical and Biological Agency, 28 Orehovuy boulevard, 115682 Moscow, Russia b
a r t i c l e
i n f o
Article history: Received 15 June 2012 Received in revised form 12 December 2013 Accepted 16 December 2013 Available online 31 December 2013 Keywords: Optimal therapy control Chemotherapy Cancer model
a b s t r a c t In this paper a mathematical model for glioma therapy based on the Gompertzian law of cell growth is presented. In the common case the model is considered with non-linear spatially varying diffusion depending on a parameter. The case of the linear spatially-varying diffusion arose as a special case for a particular value of the parameter. Effectiveness of the medicine is described in terms of a therapy function. At any given moment the amount of the applied chemotherapeutic agent is regulated by a control function with a bounded maximum. Additionally, the total quantity of chemotherapeutic agent which can be used during the treatment process is bounded. The main goal of the work is to compare the quality of the optimal strategy of treatment with the quality of another one, proposed by the authors and called the alternative strategy. As the criterion of the quality of the treatment, the amount of the cancer cells at the end of the therapy is chosen. The authors concentrate their efforts on finding a good estimate for the lower bound of the cost-function. Thus it becomes possible to compare the quality of the optimal treatment strategy with the quality of the alternative treatment strategy without explicitly finding the optimal control function. Ó 2013 Elsevier Inc. All rights reserved.
Introduction Mathematical models of brain tumor growth under consideration of diffusion have been developed and studied in many papers (see for example [31–33,16,34,3,25,27,35,17,20,21,24]). An optimal-control approach to cancer treatment, especially to brain cancer treatment, is well-established (see [8,9,11,18]). In particular, a comparison of optimal and suboptimal treatment protocols was already considered for some models of tumor anti-angiogenesis in [22]. The basic model we are developing further was suggested in [31,33] and is based on the model of Murray [26]. In these models the authors assume a linear spatially-varying diffusion and the exponential growth law for the tumor cells:
@cðx; tÞ ¼ rðDðxÞrcðx; tÞÞ þ qcðx; tÞ; @t
ðsee ½29; from p:536Þ
@cðx; tÞ ¼ rðDðxÞrcðx; tÞÞ þ qcðx; tÞ GðtÞcðx; tÞ; @t
ðsee ½36; 34Þ
⇑ Corresponding author. Tel.: +49 6212926243; fax: +49 6212926237. E-mail address: e.fi[email protected] (E. Fimmel). 0025-5564/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mbs.2013.12.007
where the function c denotes the density of tumor cells, DðxÞ has only two values DðxÞ ¼ DG , for x in gray matter and DðxÞ ¼ DW for x in white matter with DW > DG ; q 2 Rþ is a constant and GðtÞ ¼ k 2 R when the chemotherapy is being administered, and GðtÞ ¼ 0 otherwise. The model suggested in Section 1 is distinguished from the Murray et al. model mainly by the assumption of the Gompertzian law of cell growth and the non-linear spatially-varying diffusion depending on a parameter. The dynamics of the amount of the chemotherapeutic agent is in contrast to the model of Murray et al. a part of the model too. The class of the therapy functions considered is wider than in [31]. In the model suggested in this article we will adopt the Gompertzian growth law for the tumor cells since the Gompertzian models ensure that the growth rate of the tumor cells reduces as the maximal density of the tumor cells is approached while the exponential growth assumes unlimited growth of the density of tumor cells. Another significant difference between the model suggested and the model from [31,33] is a non-linear spatially-varying diffusion as in [1,29]. The case of the linear spatially-varying diffusion arose as a special case for a particular parameter value. The description of the diffusion of the malignant glioma cells is still a subject of discussion (see [1,2,13,7,14,30]). In [30] it is argued that
89
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
there is no significant difference in motility between cultures of a gross tumor and a histologically normal brain. On the other hand many authors stress an extensive non-linear character of the tumor cells distribution [12,1,13]. In [1] a non-linear diffusion in the form rðC k rCÞ depending on a parameter k is considered. For k > 0 the tissue is modeled rather as a porous medium and for k ¼ 0 as a special case with a linear diffusion studied in detail in [2]. Analogous law is usually used for the description of the diffusion in porous medium (see, for instance, [4,29]). We assume this description of the diffusion for the model suggested too. The nature of the interaction between the tumor size and the prescribed treatment is still up for debate. In the present work a wide class of increasing concave bounded therapy functions is examined. The main goal of the work is to compare the quality of the optimal strategy of treatment with the quality of another one, proposed by the authors in Section 4 and called ‘the alternative’. As the criterion of the quality of the treatment, the amount of the cancer cells at the end of the therapy is chosen. At any given moment the amount of the applied chemotherapeutic agent is regulated by a control function with a bounded maximum. Additionally, the total quantity of chemotherapeutic agent which can be used during the treatment process is bounded. Despite the essential progress achieved in the solving of optimization problems for distributed systems (see [23]) the search for the optimal control for problems described in terms of non-linear equations of parabolic type is still mathematically very difficult. Therefore the authors concentrate their efforts in Section 3 on finding a good estimate for the lower boundary of the cost-function. Thus it becomes possible to compare the quality of the optimal treatment strategy with the quality of the alternative treatment strategy without explicitly finding the optimal control function. The corridor between the analytical lower estimate and the value of the cost-function on the alternative control enables qualitative assessments and forecasts for the results of the treatment to be made.
notes the growth rate of tumor cells, clim the maximum density of the tumor cells, dh is the diffusion coefficient of the medicine, ch the dissipation rate of the medicine, GðhÞ describes the influence of the chemotherapeutic (called therapy function). We will consider the class of increasing concave bounded functions with Gð0Þ ¼ 0 e.g.,
1. Model and optimization problem formulation
(see Section 3 of the present article).
Let X R3 be a bounded domain of volume S with a smooth boundary C; n be the outer normal unit vector to C. For f 2 L1 ðXÞ we denote onward
The actual optimal control problem
f :¼
Z
GðhÞ ¼
kh ; 1þh
k 2 Rþ :
The control function uðx; tÞ 2 L1 ðX ½0; TÞ denotes the quantity of the chemotherapeutic agent applied to a patient at the moment t in x 2 X. For the cumulative amount of the chemotherapeutic agent the following inequality
Z
T
0
Z
hðx; tÞdxdt 6 Q 0
ð1:2Þ
X
with Q 0 > 0 holds. Note 1. For a ¼ 0 (respectively ah ¼ 0) the diffusions become linear
A0 cðx; tÞ ¼ rðDðxÞrcðx; tÞÞ;
A0 hðx; tÞ ¼ dh Dhðx; tÞ:
The model description is completed with the following initial conditions
cðx; 0Þ ¼ c0 ðxÞ > 0;
hðx; 0Þ ¼ 0
ð1:3Þ
and some boundary conditions which prohibit the migration of the tumor cells outside the brain boundary. It makes sense to consider two kinds of boundary conditions: either if a P 0
@cðx; tÞ j ¼ 0; @n Cð0;T
@hðx; tÞ j ¼ 0; @n Cð0;T
ð1:4Þ
or if a > 12
cðx; tÞjCð0;T ¼ 0;
@hðx; tÞ j ¼0 @n Cð0;T
ð1:5Þ
Let K # L1 ðX ½0; TÞ. In what follows we will consider three possibilities to define K:
f ðxÞdx:
X 1
Let cðx; tÞ 2 C ðX ½0; TÞ be the density of tumor cells and hðx; tÞ the amount of the chemotherapeutic in x 2 X at the moment t. We will consider the following model describing the growth of the tumor cells under the influence of chemotherapy:
– K 1 : fuðx; tÞ 2 L1 ðX ½0; TÞ j uðx; tÞ ¼ dðx x ÞuðtÞ; x 2 Xg where d be the Dirac function, – K 2 : fuðx; tÞ 2 L1 ðX ½0; TÞ j uðx; tÞ ¼ uðtÞg, – K 3 : fuðx; tÞ 2 L1 ðX ½0; TÞ j uðx; tÞ ¼ vC ðxÞuðtÞg, where
vC ðxÞ :¼
@cðx; tÞ ¼ qcðx; tÞð1 b ln cðx; tÞÞ þ Aa cðx; tÞ cðx; tÞGðhÞ @t @hðx; tÞ ¼ ch hðx; tÞ þ Aah hðx; tÞ þ uðx; tÞ; 0 6 t 6 T; @t
ð1:1Þ
where q; b; dh ; ch ; clim ; a; ah 2 Rþ are constants and b ¼ ln 1c ,
1;
x 2 C;
0;
x R C
is the characteristic function of a compact set C X of area SC . In each of the three cases the inequality
0 6 uðtÞ 6 q
ð1:6Þ
lim
Aa cðx; tÞ :¼
3 X @ @cðx; tÞ ; DðxÞ cðx; tÞ2a @xi @xi i¼1
Aa h hðx; tÞ :¼ dh
3 X @ @hðx; tÞ hðx; tÞ2ah @xi @xi i¼1
with DðxÞ ¼ Dw for white matter and DðxÞ ¼ Dg (see [26]) for gray matter. The parameters a; ah characterize the degree of the non-linearity of the respectively diffusions. The parameter q (see [26]) de-
must hold with q > 0. Consider the following optimal control problem. We need to find the time T 2 ð0; 1Þ and the control function uðx; tÞ 2 K such that uðx; tÞ; T provide the infimum to the cost-function:
Uðu; TÞ ¼
Z
ln cðx; TÞdx ¼: ln cðTÞ
X
Let
inf Uðu; TÞ ¼ Uðu ; TÞ: u2K
ð1:7Þ
90
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
In what follows (see Section 3) we will find a lower estimate U for Uðu ; TÞ. It is evident that we get an upper estimate U for Uðu ; TÞ if we substitute an arbitrary non-optimal control function u in Uðu; TÞ. Thus the inequality holds.
U 6 Uðu ; TÞ 6 U: In what follows (see Sections 4 and 5) we will suggest a suboptimal control function for each functional space K i , i ¼ 1; 2; 3, called the alternative control function so that the corresponding value of the cost-function is in some sense sufficiently close to the obtained lower estimate of the cost-function. If this holds, the results achieved can be useful in the medical prognosis.
The following Theorem helps us to understand the relationship between the important constant k (see Theorem 2) and the restrictive constant q (see (1.6)) for the dose of the drug for each one of three functional spaces considered K i ; i ¼ 1; 2; 3 : Theorem 1. Let X; h; C X be as above, G : Rþ ! R be an increasing concave bounded function. Then
GðhÞ 6 S k and (1) k :¼ GScq for uðx; tÞ 2 K 1 , h (2) k :¼ Gcq for uðx; tÞ 2 K 2 , h SC q (3) k :¼ G Sc for uðx; tÞ 2 K 3 .
2. Relations between parameters
h
Lemma 1. The following inequalities take place: hðtÞ 6 cq
(1) For uðx; tÞ 2 K 1 we have
Z
T
Z
0
hðx; tÞdxdt ¼
0
X
h
T Z
Z
Proof. For a concave integrable function f ðxÞ the Jensen inequality [15] has the form
and
t ech ðtsÞ uðsÞds dt 6 Q 0
f
0
Z
Z
T
0
Z
hðx; tÞdxdt ¼ S
hðtÞ 6 cSq
0
X
h
t
and
Z
T
0
Z
h G S
ech ðtsÞ uðsÞds dt 6 Q 0
0
hðx; tÞdxdt ¼ SC
Z
h
t
ch ðtsÞ
e
uðsÞds dt 6 Q 0
! ¼G
Z Z 1 1 1 hðxÞdx P GðhðxÞÞdx ¼ GðhÞ; S X S S X
3. A lower estimate for the cost-function
0
SC qT 2 6 : 2
Our aim in this section is to obtain a lower estimate for the objective function considered which does not depend on the choice of a control function.
Proof. Integrating the second equation from (1.1) over X and using Stokes’s formula we get
Theorem 2. Let Gðhðx; tÞÞ; c1 ðx; tÞ 2 L1 ðXÞ,
dhðtÞ ¼ ch hðtÞ þ uðtÞ dt
ln c 2 L2 ðXÞ for a ¼ 0: Also we assume for a > 0 that
with hð0Þ ¼ 0. Hence
hðtÞ ¼
Z
t
kðxÞdx X
S is the volume of the domain X. Since G is increasing and the estimate of h from Lemma 1 is valid, we obtain the upper estimates for GðhÞ. h
hðtÞ 6 ScC q and
T Z
0
X
X
Z
In our case
SqT 2 6 ; 2 (3) For uðx; tÞ 2 K 3 we have
kðxÞf ðhðxÞÞdx if kðxÞ P 0 and
¼ 1:
T Z
Z
Z P
X
qT 2 6 : 2 (2) For uðx; tÞ 2 K 2 we have
kðxÞhðxÞdx
ca ðx; tÞ 2 W 21 ðXÞ and
ech ðtsÞ uðsÞds 6
0
M
ch
ð1 ech t Þ 6
M
ch
if
uðtÞ 6 M ðÞ
2 Rþ :
Z
ln cðx; tÞdx < 1: X
Let k 2 Rþ 1 with GðhÞðtÞ 6 kS for all t 2 ½0; T and S being the volume of the domain X. Then the following estimates are valid:
Consequently we obtain: for uðx; tÞ 2 K 1 for uðx; tÞ 2 K 2 for uðx; tÞ 2 K 3
that that that
(1) for the boundary condition (1.4) with arbitrary a P 0
uðtÞ ¼ uðtÞ 6 q, uðtÞ ¼ S uðtÞ 6 Sq, uðtÞ ¼ SC uðtÞ 6 SC q.
dln c P Sðq kÞ qbln c; dt (2) for the boundary condition (1.5) and a > 12
This proves the first inequality in each case. Substituting in each RT R case uðtÞ in () we get the formulas to calculate 0 X hðx; tÞdxdt.
2 dln c Dg P Sðq kÞ qbln c þ ðln cÞ dt c1 Sa2
We get upper estimates for Q 0 in each case using the inequality
Z 0
T
hðtÞdt 6
Z 0
T
M
ch
ð1 ech t Þdt ¼
M
c
2 h
ð1 þ ech T þ T ch Þ 6
MT 2 2
if ch is sufficiently small and, in addition, Q 0 should be chosen to be as small as possible. h
where c1 is the constant from Friedrichs inequality which depends only on X.2 1 2
The number k can be viewed as a characteristic of treatment intensity. 2 2 For a rectangle X with sides a and b, for instance, c1 ¼ 24a2 b 2 [28, p. 193]. p ða þb Þ
91
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
Proof. Dividing the first equation of (1.1) by cðx; tÞ and integrating over X we get
Z
dln c ¼ qðS bln cÞ þ dt
X
Aa cðx; tÞ dx GðhÞ: cðx; tÞ
Aa cðx; tÞ 1 @ @cðx; tÞ dx dx ¼ DðxÞc2a ðx; tÞ cðx; tÞ cðx; tÞ @xi @xi i¼1 X 3 Z X @ @cðx; tÞ dx ¼ DðxÞc2a1 ðx; tÞ @xi @xi i¼1 X 3 Z X @ 1 @cðx; tÞ DðxÞc2a ðx; tÞ dx: @xi cðx; tÞ @xi i¼1 X
Z
3 Z X
X
In view of the Stokes’s theorem
@ @cðx; tÞ dx DðxÞc2a1 ðx; tÞ @xi @xi i¼1 X 3 Z X @cðx; tÞ DðxÞc2a1 ðx; tÞ cosðn; xi Þds: ¼ @xi i¼1 C
3 Z X
In the case of the boundary condition (1.4) we have 2 X i¼1
3 Aa cðx; tÞ Dg X dx P 2 cðx; tÞ a i¼1
@cðx; tÞ cosðn; xi Þ @xi
Cð0;T
@cðt; xÞ ¼ ¼0 @n Cð0;T
and in the case of the boundary condition (1.5) for a > 12 we have
c2a1 ðx; tÞjCð0;T ¼ 0:
Z
1 c1 a2
i¼1
X
@ @cðx; tÞ dx ¼ 0: DðxÞc2a1 ðx; tÞ @xi @xi
X
2 @cðx; tÞ dx @xi X a 2 3 Z 1X @c ðx; tÞ DðxÞ dx ¼ 2 @xi a i¼1 X
3 X Aa cðx; tÞ dx ¼ cðx; tÞ i¼1
Z
DðxÞc2a2 ðx; tÞ
and for a ¼ 0:
Z X
X
Aa cðx; tÞ dx ¼ cðx; tÞ
i¼1
Aa cðx; tÞ dx P 0: cðx; tÞ (1) Using GðhÞðtÞ 6 kS for all t 2 ½0; T we obtain the following inequality
dln c P Sðq kÞ qbln c dt which proves the estimate 1. (2) By virtue of the Friedrichs inequality [28]
Z X
u2 ðxÞdx 6 c1
1
Z Z
Sa2 c1
2 ca ðx; tÞdx X
2 ln cðx; tÞdx
X
where c1 depends only on X. Thus the estimate 2 is proven. h Theorem 3. For the boundary condition (1.4) with a P 0 and the condition (1.5) with a > 12 and q; S; k; clim as above, the following estimate holds
Uðu; TÞ P
Sðq kÞ ð1 eqbT Þ þ ln c0 eqbT : qb
Proof. For the boundary conditions (1.4) and (1.5) and q; S; k; clim as above the following estimate holds due to Theorem 2:
dln c P Sðq kÞ qln c: dt Let us consider the following initial condition problem
dz ¼ Sðq kÞ qbz; dt
zð0Þ ¼ ln c0 :
z¼
Sðq kÞ ð1 eqbt Þ þ ln c0 eqbt : qb
Using the comparison Lemma [19, p. 11] and Theorem 2 we get the result. h Note 2. In the case of the boundary condition (1.5) for a > 12 we can get a better lower estimate for Uðu; TÞ. In this case we deal with a Riccati equation
dz Dg 2 ¼ Sðq kÞ qbz þ z :¼ PðzÞ: dt c1 Sa2 D
2 @ ln cðx; tÞ DðxÞ dx: @xi X
3 Z X
In both cases we obtain
Z
X
1 Sc1 a2
The solution of this problem is
Thus we have for a – 0:
Z
c2a ðx; tÞdx P P
In both cases we get 3 Z X
2 Z a Z @c ðx; tÞ Dg dx P 2 c2a ðx; tÞdx: @xi a c1 X X
Moreover, in the view of the Cauchy–Bunyakovsky–Schwarz inequality we get
Aa cðx; tÞ dx: cðx; tÞ
X
Z X
Now let us estimate
Z
for an arbitrary function u 2 W 21 ðXÞ and some constants c1 ; c2 P 0 depending only on X, we obtain the following estimate in the case of the boundary condition (1.5) for a > 12:
2 Z 3 Z X @uðx; tÞ dx þ c2 u2 ðsÞds @xi C i¼1 X
Let us denote a0 :¼ Sðq kÞ; a1 :¼ qb; a2 :¼ c Sga2 . We obtain the 1 differential equation
dz ¼ a0 þ a1 z þ a2 z2 dt
with zð0Þ ¼ ln c0 :
Let us note first that only the situation
4a0 a2 þ a21 ¼ 4
Dg ðq kÞ þ q2 b 2 > 0 a2 c1
can occur in our case as we will see below in Section 4 due to the condition q < k. The roots of PðzÞ
z1;2 ¼
a1
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4a0 a2 þ a21 2a2
;
z1 < 0 < z2 ;
are particular solutions of the Riccati equation. As we see below, z1 is an attractor and z2 is a repeller. We set
y¼
1 : z z1
92
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
We get a linear differential equation
dy ¼ ða1 þ 2a2 z1 Þy a2 ; dt
yð0Þ ¼
1 ln c0 z1
:
The solution of this equation in the case a21 > 4a0 a2 (in particular, it means that a1 þ 2a2 z1 – 0) is given by
yðtÞ ¼
a2 a2 1 eða1 þ2a2 z1 Þt : þ þ a1 þ 2a2 z1 a1 þ 2a2 z1 ln c0 z1
Thus the solution of the equation is given by
1
þ 1 eða1 þ2a2 z1 Þt ln c0 z1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a21 4a0 a2 ¼ z1 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a2 4a a a2 þ a2 þ 1 0 2 e 4a0 a2 þa1 t
zðtÞ ¼ z1 þ
a2 þ a1 þ2a 2 z1
a2 a1 þ2a2 z1
consider the cases a ¼ 0; a ¼ 0:05 (for the Neumann boundary condition), and ah ¼ 0 and respectively ah ¼ 0:1. We need an estimate for the dissipation rate of the chemotherapeutic agent. Studying the pharmakinetics description of a drug with the most simple application schedule we note that the halfdecay time for men is 20 days. Thus we get ech 20 ¼ 12 or, respec1 tively, ch ¼ 0:0347 day . From the pharmakinetics description of the same drug we can also learn that a single dose administered in glioma case is 10 mg/kg. So we can estimate the parameter q in the case uðx; tÞ ¼ uðtÞ 2 K 2 as q2 ¼ 0:8. We assume for the control functions of all types that the product of the application area of the medicine and q is the same in all three cases. Thus we get
q2 S ¼ q3 SC ¼ q1 102 and finally q1 ¼ 2880; q3 ¼ 2:ð2Þ.
ln c0 z1
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 4a a þa2 with z1 ¼ 1 2a2 0 2 1 < 0. It follows for the boundary condition (1.5) and a > 12 in view of the comparison Lemma [19, p.11] that
Uðu; TÞ P zðTÞ: Evidently,
lim Uðu; TÞ P lim zðTÞ ¼ z1 :
T!1
T!1
However, the numerical calculations show that this more precise lower estimate is only marginally better than the estimate from Theorem 3.
Numerical Results A tumor can be usually diagnosed if its radius is equal to 15 mm. The maximum tumor extent radius is 37.5 mm ([20,6]). Thus we take as a simulation space the square with the side length a ¼ 6cm thinking as a 2-dimensional projection of the domain X to visualize the results in the 3-dimensional space. Further we represent the square as a 60 60-mesh consisting of equal squares having the area 1mm2 . As a time unit we take 1 day. We consider the following initial spatial density of the tumor cells ([26, p.547]):
( c0 ðx; yÞ :¼
ðxx0 Þ2 ðyy0 Þ2 B
Ae
;
1;
ðx x0 Þ2 þ ðy y0 Þ2 6 R20 else
3
4. Alternative control Formulation of the alternative control task
with A ¼ 2 10 ; R0 ¼ 1:8; B ¼ 0:6. Below we use for the numerical calculations ðx0 ; y0 Þ ¼ ð3; 3Þ: We consider for our numerical simulations the therapy function (see Fig. 1)
We will search alternative control functions in each functional space K i ; i ¼ 1; 2; 3 on the base of the following consideration: ~ ðx; tÞ 2 K i ; i ¼ 1; 2; 3, satisfying we look for such a function u (1.6), where the inequality (1.2) holds and for the functional
GðhÞ ¼
Fðu; TÞ ¼
Z
c0 ðxÞGðhðx; TÞÞdx
ð4:1Þ
X
kh ; k 2 Rþ 1þh
and calculate k in each of the three cases depending on a treatment intensity k is considered from
kS ¼ kmax
the equality
t2½0;T
Z X
h dx: 1þh
holds.
~ ðtÞ (see (4.1)) in the class We will search the proper functions u of the piecewise constant functions on the interval ½0; T dividing the whole interval into equal parts:
Parameter estimates
~ ðtÞ ¼ ui ; u
It will be useful in what follows to have some parameter values. We adopt the parameter estimates from [27,26,32]:
As a result of numerical calculations we obtain the following alternative control function in all three cases considered
~ ; TÞ ¼ sup Fðu; TÞ Fðu u
– Diffusion coefficient (Gray matter-White matter) 1
Dg ¼ 0:0013 cm2 day ;
~ ðx; tÞ ¼ u 1
Dw ¼ 0:005 cm2 day ;
– Maximal density of the tumor cells
– Growth rate of the tumor cells 1
q ¼ 0:012 day ; 1
In the case ah ¼ 0 we set dh ¼ 0:0864 cm2 day equal to the apparent diffusion coefficient (ADC) for the water. We will
ði 1Þ T i T ; ; n n
q;
0 6 t < t0 ;
0;
t0 6 t 6 T
i 2 f1; . . . ; ng;
ui 6 q:
where the parameter t 0 can be found from the equation
Z
clim ¼ 105 cells cm2 ) ln clim ¼ 11:5; b ¼ 0:0869
t2
t0
hðtÞdt ¼ Q 0 :
ð4:2Þ
0
Note that the moment of time t 0 and, as a consequence, the ~ ðx; tÞ depend on the parameter q alternative control function u and Q 0 as well as on the chosen time of the end of the therapy T. Let us also add that for the ODE systems in a similar case it is proven that the control function given above is the optimal control ([10,5]).
93
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
2000
c
1500
1000
500
0 6 5
6
4
5 3
4 3
2
2
1 0
y
1 0
x
Fig. 1. Initial conditions (distribution of cancer cells density) for the control task.
160 140
k=0 u=u(t)
∫ Ω ln c(x,T) dx
120 100 80 u=δ(x*)u(t)
60
u=χR’(x*)u(t) u=χ (x*)u(t) R’
40 20
0
10
20
30
40
50
60
70
80
90
T Fig. 2. The numerical upper and analytical lower (the lower black line) estimates for three control functions types in the case of the Neumann boundary conditions: for uðx; tÞ ¼ dðx x ÞuðtÞ (blue line), for uðx; tÞ ¼ vC ðxÞuðtÞ (green and yellow lines, respectively for different compact sets) and for uðx; tÞ ¼ uðtÞ (red line), a ¼ 0; ah ¼ 0; k ¼ 0:0196. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Let us compare the results shown in the pictures 2–5. The pictures 2–4 illustrate the case of the usual Laplace operator (a ¼ ah ¼ 0). The difference between them is in the increasing middle intensity of the therapy (k) from picture 2 to picture 4. The upper black line shows the growth of the cost function without any therapy administered for comparison aims in every picture. Surprisingly, we can see that while the increasing therapy intensity has little effect on the results in the cases uðx; tÞ 2 K 2 or uðx; tÞ 2 K 3 , in the case uðx; tÞ 2 K 1 it looks significantly better. Fig. 4 shows, in particular, a relative deceleration of the tumor reduction in the case of uðx; tÞ 2 K 1 . In this case the time moment t0 (compare 4.2) of the end of the active therapy for uðx; tÞ 2 K 1 occurs much earlier than the intended end of the therapy T due to the higher value than for uðx; tÞ 2 K 2 ; K 3 of the parameter q for this kind of control. The effect is more visible than in Fig. 2,3 due to the higher intensity of the application of the chemotherapeutic agent. The results obtained show, in particular, that the restriction on the cumulative amount as well as the intensity of the application of
the chemotherapeutic agent have a big influence on the therapy process which emphasises the need for further medical investigations in this area (see Fig. 5). The picture 5 shows the case of a small non-linearity in the diffusion (a ¼ 0:05; ah ¼ 0:1). In this case the results of the therapy in all cases are significantly worse than in the case of the linear diffusion because of a very extensive growth of the tumor. Let us also note that the lower estimate, calculated analytically because of its generality is relatively rough. The Fig. 6(a)–(c) shows numerically simulated change of the tumor after 90 days with the application of different therapy protocols. It is important to compare the scales on the axis c in Fig. (a)–(c). While the maximum value is still at about 1000 units after an application of the protocol 1, in the last picture we have the maximum value of only about 3 units. If we compare the initial volume of the tumor (with its maximum value at approximately 1500 units) and its volume after 90 days of therapy with three protocols considered, we can see that the
94
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
200 k=0
150
∫ ln c(x,T) dx
u=u(t) 100
50
Ω
u=δ(x*)u(t)
u=χR’(x*)u(t) u=χ (x*)u(t) R’
0
−50
0
10
20
30
40
50
60
70
80
90
T Fig. 3. The numerical upper and analytical lower estimates for three control functions types in the case of the Neumann boundary conditions: for uðx; tÞ ¼ dðx x ÞuðtÞ (blue line), for uðx; tÞ ¼ vC ðxÞuðtÞ (green and yellow lines, respectively for different compact sets) and for uðx; tÞ ¼ uðtÞ (red line), a ¼ 0; ah ¼ 0; k ¼ 2 0:0196. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
200 150
k=0
u=u(t)
∫ Ω ln c(x,T) dx
100 50 0
u=δ(x*)u(t)
u=χR’(x*)u(t) u=χ (x*)u(t) R’
−50 −100 −150
0
10
20
30
40
50
60
70
80
90
T Fig. 4. The numerical upper and analytical lower (the lower black line) estimates for three control functions types in the case of the Neumann boundary conditions: for uðx; tÞ ¼ dðx x ÞuðtÞ ((blue line)), for uðx; tÞ ¼ vC ðxÞuðtÞ (green and yellow lines, respectively for different compact sets) and for uðx; tÞ ¼ uðtÞ (red line), a ¼ 0; ah ¼ 0; k ¼ 3 0:0196. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
180 α=0.05; α =0.1
160
k=0
u(x,t)=u(t) u(x,t)=χR’(x*)u(t)
h
∫ Ω ln c(x,T) dx
140 120 100 u(x,t)=χΩ’(x*)u(t)
80 60
u(x,t)=δ(x*)u(t)
40 20
0
10
20
30
40
50
60
70
80
90
T Fig. 5. The numerical upper and analytical lower (the lower black line) estimates for three control functions types in the case of the Neumann boundary conditions: for uðx; tÞ ¼ dðx x ÞuðtÞ (continuous line), for uðx; tÞ ¼ vC ðxÞuðtÞ (dotted line) and for uðx; tÞ ¼ uðtÞ (dashed line), a ¼ 0:05; ah ¼ 0:1; k ¼ 0:0196.
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
(a)
1000
c
T
800 600 400 200 0 4 4
3 3
2
2
1
(b)
1
0 0
y
x
250 200 150 100 50 0 0 1 2 3 4 5
x
1
6 0
4
3
2
5
6
y
95
moment the amount of the chemotherapeutic agent applied is regulated by a control function with a bounded maximum. Additionally, the total quantity of chemotherapeutic agent which can be used during the treatment process is also bound. The main goal of the work is to compare the quality of the optimal strategy of treatment with the quality of some other treatment control strategies, proposed by the authors, called the alternative. As the criterion of the quality of the treatment (the cost-function), the amount of cancer cells at the end of the therapy is chosen.The alternative control strategy is presented for three types of control function. In the first case the dose of the chemotherapeutic agent varies depending only on time. In the second case the application of the chemotherapeutic agent is targeted on the tumor and the third one models the case of very precise and intensive administered application of the chemotherapeutic agent comparable with an injection. The results obtained show the preciser the delivery of the chemotherapeutical agent is targeted on the tumor the better the result of the therapy is. In the work the lower boundary of the cost-function is analytically estimated. Thus it becomes possible to compare the quality of the optimal treatment strategy with the quality of the alternative treatment strategy without solving the optimal control problem. Thus, in the work a corridor for the cost-function has been constructed which enables observation of the behavior of the tumor during the therapy. The lower boundary of the corridor is an analytical estimate. The upper boundary of this corridor is determined by the proposed suboptimal control functions. Since the proposed alternative control strategy is essentially easier to calculate as the optimal one, it enables qualitative assessments and forecasts for the results of the treatment to be made. Acknowledgments The authors thank Dr.Yu.S. Semenov (MIIT, Moscow, Russia), whose expertise provided valuable quality control. The authors are also grateful to Dr. G. Schmalz (UNE, Armidale, Australia) for his valuable and constructive comments.
(c) c(x,y,T=90)
3 2.5
References
2 1.5 6
1 0
4 2
2
4
x
6 0
y
Fig. 6. The graphs for cðx; TÞ with T ¼ 90 days and ah ¼ 0; a ¼ 0. (a) u ¼ uðtÞ ðq ¼ 0:8Þ (b) u ¼ vC ðx; x ÞuðtÞ ðq ¼ 2:8294Þ (c) u ¼ dðx ÞuðtÞ ðq ¼ 2880Þ.
biggest positive change by far is obtained after an application of the protocol 3 (Fig. 6(c)).
5. Conclusion Presented in the work is a mathematical spatial model for glioma therapy consisting of two non-linear partial differential equations of parabolic type with non-linear spatially varying diffusion depending on a parameter. The case of the linear spatially-varying diffusion arose as a special case for a particular value of the parameter. The first equation describes the dynamics of cancer cell growth under Gompertzian law and the interaction between a controlled quantity of the chemotherapeutic agent and cells given in terms of a therapy function. The second differential equation describes the dynamics of the chemotherapeutic agent. At any given
[1] E.K. Afenya, C.P. Calderón, Growth kinetic of cancer cells prior to detection and treatment an alternative view, Discrete Contin. Dyn. Syst. – Ser. B 4 (1) (2004). [2] E.K. Afenya, C.P. Calderón, Diverse ideas on the growth kinetics of disseminated cancer cells, Bull. Math. Biol. (2000). [3] M. Assefa, R. Rockne, M. Szeto, K.R. Swanson, Mathematical modeling of Glioma proliferation and diffusion, Ethnicity and Disease 19 (2, Suppl. 3) (2009) 60–61. Pubmed ID: 19554787. [4] G.I. Barenblatt, On some unsteady fluid and gas motions in a porous medium, Prikladnaya Matematika i Mekhanika (Applied Mathematics and Mechanics (PMM)) 16 (1) (1952) 67 (in Russian). [5] A.S. Bratus, E.S. Chumerina, Optimal control in therapy of solid tumor growth, Comput. Math. Math. Phys. 48 (6) (2008) 892–911. [6] P.K. Burgess, P.M. Kulesa, J.D. Murray, E.C. Alvord Jr., The interaction of growth rates and diffusion coefficients in a three-dimensional mathematical model of Gliomas, J. Neuropathol. Exp. Neurol. 56 (6) (1997) 704. [7] A.C. Burton, Rate of growth of solid tumor as a problem of diffusion, Growth 30 (1966) 157. [8] S.P. Chakrabarty, F.B. Hanson, Optimal control of drug delivery to brain tumors for a distributed parameters model, in: Proceedings of American Control Conference, 2005, p. 1. [9] M. El-Borai, Application of stochastic control theory in modeling of brain cancer, in: Casablanca International Workshop on Mathematical Biology: Analysis & Control, June 20–24, 2011. [10] E. Fimmel, Y.S. Semenov, A.S. Bratus, On optimal and suboptimal treatment strategies for a mathematical model of Leukemia, Math. Biosci. Eng. 10 (1) (2013) 151. [11] K.R. Fister, J.C. Panetta, Optimal control applied to competing chemotherapeutic cell-kill strategies, SIAM J. Appl. Math. 63 (6) (2003) 1954. [12] A. Friedman, Therapeutic approach to brain cancer, in: Abstracts of 8th European Conference on Mathematical and Theoretical Biology, Krakov, June 28–July 2, 2011, p. 318. [13] R.A. Gatenby, E.T. Gawliski, The glycolytic phenotype in carcinogenesis and tumor invasion: insight thought mathematical models, Cancer Res. 63 (15) (2003) 3847.
96
A.S. Bratus et al. / Mathematical Biosciences 248 (2014) 88–96
[14] V. Greenspan, Models for the growth of solid tumor as a problem by diffusion, Stud. Appl. Math. 52 (1972) 317. [15] G.H. Hardy, J.E. Littlewood, G. Polya, Inequalities, Cambridge University Press, 1934. [16] H.L.P. Harpold, E.C. Alvord Jr., K.R. Swanson, The evolution of mathematical modeling of glioma growth and invasion, J. Neuropathol. Exp. Neurol. 66 (1) (2007) 1. [17] T. Hines, Mathematically modeling the mass-effect of invasive brain tumors, SIURO 3 (2010). Available from: . [18] M. Itika, M.U. Salamcib, S.P. Banksa, Optimal control of drug therapy in cancer treatment, Nonlinear Anal. Theory Methods Appl. 71 (12) (2009) e1473. [19] E. Kamke, Differentialgleichungen Lösungsmethoden und Lösungen, Band 1, Gewöhnliche Differentialgleichungen, Chelsea Publishing Company, New York, 1971. [20] A.R. Kansal, S. Torquato, G.R. Harsh IV, E.A. Chiocca, T.S. Deisboeck, Simulated brain tumor growth dynamics using a three-dimensional cellular automaton, J. Theor. Biol. 203 (2000) 367. [21] E. Khain, L.M. Sander, A.M. Stein, A model for glioma growth, Complexity 11 (2005) 53. [22] U. Ledzewicz, H. Schaettler, Optimal and suboptimal protocols for a class of mathematical models for tumor anti-angiogenesis, J. Theor. Biol. 252 (2008) 295. [23] X. Li, Y. Yong, Optimal control Theory for Infinite Dimensional Systems, Birkhauser, Boston, 1995. [24] D. Mackenzie, Mathematical modeling and cancer, SIAM News 37 (1) (2004). [25] A. Martínez-González, G.F. Calvo, L.A. Pérez Romasanta, V.M. Pérez-García, Hypoxic cell waves around Necrotic cores in Glioblastoma: a biomathematical model and its therapeutic implications, Bull. Math. Biol. 62 (3) (2012) 527.
[26] Murray, Mathematical biology, vol. 2, Spatial Models and Biomedical Applications, 3rd ed., Springer, 2003. [27] G. Powathil, M. Kohandel, S. Sivaloganathan, A. Oza, M. Milosevic, Mathematical modeling of brain tumors: effects of radiotherapy and chemotherapy, Phys. Med. Biol. 52 (2007) 3291. [28] K. Rektoris, Variational Methods in Mathematics: Science and Engineering, Mir, Moscow, 1985 (in Russian). [29] A.A. Samarski, V.A. Galaktionov, S.P. Kurdyumov, A.P. Mikhailov, Blow-up in quasilinear parabolic equations, de Gruyter Expositions in Mathematics, vol. 19, de Gruyter, Berlin and Hawthorne, NY, 1995. [30] D.L. Silbergeld, M.R. Chicoine, Isolation and characterization of human malignant glioma cells from histologically normal brain, J. Neuroseurg. 86 (1997) 525. [31] 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. [32] 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. [33] K.R. Swanson, E.C. Alvord Jr., J.D. Murray, Quantifying efficacy of chemotherapy of brain tumors (Gliomas) with homogeneous and heterogeneous drug delivery, Acta Biotheor. 50 (4) (2002) 223. [34] K.R. Swanson, H.L.P. Harpold, D.L. Peacock, R. Rockne, C. Pennington, L. Kilbride, R. Grant, J. Wardlaw, E.C. Alvord Jr., Velocity of radial expansion of contrast-enhancing Gliomas and effectiveness of radiotherapy in individual patients: a proof of principle, Clin. Oncol. 20 (2008) 301. [35] 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.
