- Email: [email protected]
Determination of an optimal control strategy for drug administration in tumor treatment using multi-objective optimization differential evolution
Accepted Manuscript Title: Determination of an optimal control strategy for drug administration in tumor treatment using multi-objective optimization differential evolution Author: Fran Sérgio Lobato, Vinicius Silvério Machado, Valder Steffen Jr PII: DOI: Reference:
S0169-2607(15)30460-0 http://dx.doi.org/doi: 10.1016/j.cmpb.2016.04.004 COMM 4130
To appear in:
Computer Methods and Programs in Biomedicine
Received date: Revised date: Accepted date:
25-12-2015 17-3-2016 6-4-2016
Please cite this article as: Fran Sérgio Lobato, Vinicius Silvério Machado, Valder Steffen Jr, Determination of an optimal control strategy for drug administration in tumor treatment using multi-objective optimization differential evolution, Computer Methods and Programs in Biomedicine (2016), http://dx.doi.org/doi: 10.1016/j.cmpb.2016.04.004. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Determination of an Optimal Control Strategy for Drug Administration in Tumor Treatment using Multi-objective Optimization Differential Evolution Fran Sérgio Lobato1, Vinicius Silvério Machado2, Valder Steffen Jr3 1
School of Chemical Engineering, FEQUI,
2
School of Medicine, FAMED,
3
School of Mechanical Engineering, FEMEC
Universidade Federal de Uberlândia, UFU, P.O. Box 593, 38400-902, Uberlândia, Brazil. URL: http://www.feq.ufu.br/ (Fran Sérgio Lobato) Highlights •
To formulate and to solve a multiobjective optimal control problem applied to treatment
of tumors. •
Objectives: to minimize the cancerous cells concentration and the drug concentration.
•
The results demonstrated that an optimal protocol can be obtained considering both the
objectives. Abstract The mathematical modeling of physical and biologic systems represents an interesting alternative to study the behavior of these phenomena. In this context, the development of mathematical models to simulate the dynamic behavior of tumors is configured as an important theme in the current days. Among the advantages resulting from using these models is their application to optimization and inverse problem approaches. Traditionally, the formulated Optimal Control Problem (OCP) has the objective of minimizing the size of tumor cells by the end of the treatment. In this case an important aspect is not considered, namely, the optimal concentrations of drugs may affect the patients’ health significantly. In this sense, the present work has the objective of obtaining an optimal protocol for drug administration to patients with cancer, through the minimization of both the cancerous cells concentration and the prescribed drug
1
[email protected] (Corresponding author) [email protected] 3 [email protected] 2
Page 1 of 22
concentration. The resolution of this multi-objective problem is obtained through the Multiobjective Optimization Differential Evolution (MODE) algorithm. The Pareto’s Curve obtained supplies a set of optimal protocols from which an optimal strategy for drug administration can be chosen, according to a given criterion. Keywords: Mathematical Modeling of Tumors, Differential Evolution Algorithm, Optimal Control Problem, Multi-objective Optimization. 1.
Introduction In the last decades, despite of numerous efforts dedicated to the development of new
technologies for tumors diagnosis and treatment, the incidence of cancer has increased significantly. The ideal treatment for this disease has not yet been found: to destroy tumor cells without affecting normal cells, and inhibit its aggressive manifestation [1]. As mentioned by Pillis and Radunskaya [2], a tumor’s response depends on many factors, including the severity of the disease, the efficiency of the treatment, and the strength of patient’s own immune response, to cite only a few. In this context, the development of mathematical models to represent the tumor growth is configured as an interesting research area. The objective is to find the best protocol for drug administration, although the growth of a cancerous tumor is a very complicated process involving complex interactions [2–4]. In the literature, various mathematical models for cancer treatment have been proposed due to better modeling capabilities and due to the refinement of the techniques that can be used to estimate the necessary control parameters. Besides, various contributions are focused on analyzing the problems including chemotherapy treatment as a mean to deplete cancer cells, immunotherapy as a way to boost the immune system, as well as a combination of the previous mentioned treatments [2–16]. The determination of optimal protocol for drug administration characterizes an Optimal Control Problem (OCP). In the literature, the OCP is solved by using the Optimal Control Theory, with applications in different fields of science and engineering, such as aerospace, process control, robotics, bioengineering, economics, finance, and management science [17–22]. The OCP consists in the determination of the control variable profiles that minimize (maximize) a given performance index as characterized by a differential index, i.e., the number of times required to differentiate the Differential Algebraic-Equation system (DAE) to obtain a system of Ordinary Differential Equations (ODE) [17; 22]. For this problem, when the differential index fluctuation occurs as due to the activation and deactivation of inequality constraints or when the
Page 2 of 22
control variable behaves linearly in the Hamiltonian function, this is called Singular Optimal Control Problem (SOCP) [17–22]. Several numerical methods have been proposed to solve SOCP [17–19], which are classified according to three broad categories: direct optimization methods, Pontryagin’s Maximum Principle (PMP) based methods and HJB based (Hamilton-Jacob-Bellman) methods. Among these methods, the direct approach has been preferentially used in the last several years. This approach uses control parameterization or state and control parameterizations, transforming the original problem into a finite dimensional optimization problem. Recently, algorithms based on heuristic approaches have been used to solve the SOCP due to their success observed in solving general optimization problems. In this context, the Differential Evolution Algorithm (DE), proposed by Storn and Price [23; 24], has been applied successfully. Kapadi and Gudi [25] determined the substrate concentration profiles in a feed-batch reactor with singular arc. Lobato et al. [26] presented a new algorithm for dealing with control optimization problems. The proposed methodology consists in the extension of the DE to problems with multiple objectives, through the incorporation of special mechanisms such as rank ordering and neighborhood potential solution-candidates exploration. This algorithm is applied to determine the switching times (events) and operation time of the lysine fermentation process. Wang and Chiou [27] proposed an algorithm based on DE to determine the optimal control and optimal time location problems of differential-algebraic systems. Chowdhury et al. [28] presented a hybrid evolutionary direct search technique based on DE to solve optimal control problems. Lobato et al. [29] proposed a strategy for the dynamic updating of the population size to reduce the number of objective function evaluations. This strategy is based on the definition of convergence rate to evaluate the homogeneity of the population in the evolutionary process. It has been applied to the solution of singular optimal control problems in chemical and mechanical engineering. Traditionally, the main objective considered in the determination of optimal protocol for drug administration is the minimization of the size of tumor cells at the end of the treatment [2– 4]. In this case, the high drug concentrations in the organism, which affects significantly the patient’s health, is not taken into account. In this contribution, the main goal is to introduce a systematic methodology to minimize simultaneously the cancerous cells concentration and the drug concentration during the treatment by using the Multi-objective Optimization Differential Evolution (MODE) algorithm [30]. The work is organized as follows. Section 2 presents the
Page 3 of 22
mathematical model of tumor growth. In Section 3 the general aspects regarding the formulation and solution of SOCP is presented. A definition of the Multi-objective Optimization Problem (MOP) is presented in Section 4. A review about DE and its extension to deal with multi-criteria optimization are presented in Sections 5 and 6, respectively. In Section 7, the proposed methodology is discussed. The results obtained from the application of the methodology conveyed as based on two test-cases are presented in Section 9. Finally, the conclusions are outlined in Section 9. 2. Mathematical Modeling of Tumor Growth In this work, the mathematical model proposed by Pillis and Radunskaya [2] is adopted. This model has the following characteristics: Immune Response: includes immune cells whose growth may be stimulated by the presence of the tumor and that can destroy tumor cells through a kinetic process; Competition Terms: normal cells and tumor cells compete for available resources, while immune cells and tumor cells compete in a predator-prey fashion; Optimal Control Theory for Chemotherapy: a set of optimal drug therapies is calculated through the minimization of cancerous cells concentration during the treatment, while keeping the normal cells in a condition above a required level; Logistic Growth Law. In the model proposed by Pillis and Radunskaya [2], I denotes the number of immune cells at time t, T represents the number of tumor cells at time t, N stands for the number of normal (host) cells at time t, and u is the control strategy. The source of the immune cells is considered to be outside the system, so that its reasonable to assume a constant influx rate s (steady source rate for immune cells in the absence of a tumor, 0 ≤ s ≤ 0.5). In the absence of any tumor, the cells will die off at a per capita rate dl, resulting in a long-term population size of s/d1 cells. The presence of tumor cells stimulates the immune response as represented by the positive nonlinear growth term for the immune cells ρ IT α T
where ρ and α are positive constants. The tumor cells as well as the normal cells are modeled by a logistic growth law, with parameters ri (r1 may be bigger or smaller than r2, depending on the type of cancer or stage of growth) and bi ( b1 1
1
b2
1)
representing the per capita growth rates
Page 4 of 22
and reciprocal carrying capacities of the two types of cells (i = 1, 2 and 3 identifies those associated with the cells tumor, normal tissue and immune, respectively). Mathematically, this model is described by the following system of differential equation: N r2 N (1 b 2 N ) c 4 T N a 3 u
N (0 ) N
T r1T (1 b 1T ) c 2 I T c 3 T N a 2 u I s
ρ IT α T
c1 I T d 1 I a 1u
o
(1)
T (0 ) To
I (0 ) I o
(2)
(3)
where the parameters ai are the cell death rates (0 ≤ ai ≤ 0.5, with a3 ≤ a1 ≤ a2), ci are the competition rates (ci are considered positive), d1 is the mortality rate, α is related with inverse declivity of immune response (0 ≤ α ≤ 0.5) and ρ is the immune response rate (0 ≤ ρ ≤ 1). In this model, the term T N represents the probability of finding tumor and normal cells, and the term IT describes the probability of finding T and I [2]. It is important to emphasize that, since the model used is qualitative and does not focus on a particular tumor type, its not immediately apparent how to measure which growth law is preferable in this context. However, as mentioned by the authors, the choice of growth law does not significantly affect the qualitative behavior of the model. To better understand the dynamics of the phenomenon, Pillis and Radunskaya [2] analyzed the system without and with drug input. For the first case, the authors demonstrated that depending on the values of the parameters, there could be zero, one, two or three equilibria states. In this analysis tumor-free equilibrium, dead equilibria and coexisting equilibria were considered. More details about the corresponding stability analysis can be found in [2]. 3.
Formulation of the Singular Optimal Control Problem Mathematically, the SOCP can be formulated as follows [17–19]: tf
m in J ( z ( t t ) , t f
u ( t ),t f
L ( z, u , t)dt
(4)
t0
where z is the vector of state variables and u is the vector of control variables. Ψ and L are the first and second terms of the performance index, respectively. The objective is subject to the implicit Differential-Algebraic Equations (DAE) system as given by: f (z, z, u, t) 0 g (z, u, t) 0
(5)
(6)
Page 5 of 22
p (u , t ) 0
q(z, u, t)
(7) 0
ttf
(8)
with consistent initial conditions given by: φ ( z (t0 ), z (t0 ), u (t0 ), t0 ) 0
where J(.), L(.),
.
; f (.),
(9) φ ( .)
mz
;
z
mz
;
u
mu
;
g
mg
;
p
m
p
and
q
mq
.
According to the Optimal Control Theory [17; 18], the solution of the SOCP, defined by Eq.(4) to Eq.(9), is satisfied by the co-state equations and the stationary condition given, respectively, by: λ
T
H u
H
λ (t f )=
z
0
z
(10) ttf
(11)
where H is the Hamiltonian function defined by: H L
T
f
(12)
The system formed by Eq.(10) to Eq.(12) is known as the Euler-Lagrange equations, which are characterized as Boundary Value Problems (BVPs). According to Bryson and Ho [17] and Feehery [18], the main difficulties associated with the SOCP solution are the following: the existence of end-point conditions (or region constraints) implies multipliers and associated complementary conditions that significantly increase the difficulty of solving the BVP by the indirect method; the existence of constraints involving the state variables and the application of slack variables method may originate DAE of higher index, regardless the constraint activation status, even in problems where the number of inequality constraints is equal to the number of control variables; the Lagrange multipliers can be very sensitive to the initial conditions. 3.1.
Singular Arcs The solution of SOCP presents special challenges since it demands the knowledge of the
sequence and the number of constraint activations and deactivations (events) along the trajectory. When the amount of constraints is reduced, its usually possible to determine this sequence by examining the solution of the problem without constraints. However, the presence of a large number of restrictions brings a problem of combinatorial nature [18].
Page 6 of 22
In the events, due to discontinuities in the state and/or the co-state variables, changes in the functional form of the DAE and/or in the trajectories of the control variable in each phase may occur. As a consequence, the differential index of the system can change throughout the solution trajectory, increasing when the inequality becomes active. The existence of sections of fluctuating index leads to different errors in these sections, demanding the application of adjusted numerical strategies for each section. Therefore, its necessary to know previously the moments of activation and deactivation of the restrictions in order to solve the problem adequately. Another difficulty is the presence of singular arcs, where the second derivative matrix of the Hamiltonian with respect to the control is only positive semi-definite [17; 18]. A particular case of great interest is the one that appears when the control variable behaves linearly in the Hamiltonian function. In general, no minimum optimal solution will exist for such problems unless inequality constraints in the state and/or control are specified. If the inequality constraints are linear with respect to the control variable, its reasonable to expect that the minimum solution, if it exists, will always impose that the control variables are located at a point belonging to the border of the viable region of control [17; 19; 22]. Consider the following system of equations: z F1 ( z ) F 2 ( z ) u
z (to ) z o
(13)
with the control variable given by: u m in u u m a x
(14)
The Hamiltonian function is defined as: T
H λ ( F1 ( z ) F 2 ( z ) u )
(15)
For this class of control we have: u m ax u u m in
T
λ F2 0 T
λ F2 0
(16)
T
λ F2 0
where ℑ is the Switching Function [18; 19]. Next section presents basic concepts about multi-objective optimization, as well general aspects about DE and its extension to the multi-objective context. 4.
Multi-objective Optimization Problem
Page 7 of 22
The Multi-objective Optimization Problem (MOP) is very common in different areas such as mathematics, engineering and sciences in general. This problem is different from that of a single-objective optimization since MOP usually has not only one but a set of non-inferior optimal solutions, known as Pareto’s solutions (that form the so-called Pareto’s curve). Besides, the optimality concept for this kind of problem is different from the one used in single optimization problems since in multi-objective optimization its not possible to find a single optimal solution that satisfies all the goals, simultaneously [31]. The notion of optimality for MOP is different from the one used for single optimization problems. The most common idea about multi-objective optimization found in the literature was originally proposed by Francis Ysidro Edgeworth [32] and later generalized by Vilfredo Pareto [33]. This idea can be described as follows [33]: “a solution is optimal if its dominated by no other feasible solution, which means that there exists no other solution that is superior at least in case of one objective function value, and equal or superior with respect to the other values of the objective functions”. This definition leads to finding a set of solutions that is called the Pareto optimal set, whose corresponding elements are called non-dominated or non-inferior solutions. In the literature, there are several methods available in the literature for solving MOPs [31; 34]. These methods follow a preference-based approach, in which a relative preference vector is used to scalarize multiple objectives. As classical optimization methods use a point to point approach, in which a solution found for a given iteration is modified to obtain a new solution, the outcome of using a classical optimization method is a single optimized solution. On the other hand, the evolutionary algorithms can find multiple optimal solutions in one single simulation run due to their population-based search approach, being this way well adapted for multi-objective optimization problems. More details regarding the corresponding methodologies can be found in the literature [31; 34]. Mathematically, the MOP is defined as [31]: m in [ f 1 ( x ) f 2 ( x )
f k ( x )]
(17)
subject to m inequality and l equality constraints g i ( x ) 0 , i 1, .., m
(18)
h j ( x ) 0 , j 1, .., l
(19)
and design space defined according to the decision (or design) variables
Page 8 of 22
U
L
x n x n x n , n 1,
, N.
(20)
where k is the number of objective functions fi:
N
.
In the last decades, several algorithms based on evolutionary mechanisms have been proposed to find approximations to the Pareto optimal solutions. As reported in the literature, the first Multi-objective Evolutionary Algorithm (MEA) was the Vector Evaluated Genetic Algorithm [35]. On one hand there is a first group known as the first-generation, which includes all the early MEAs: Vector Evaluated Genetic Algorithm [35], Non-dominated Pareto Genetic Algorithm [37], Non-dominated Sorting Genetic Algorithm [38]. On the other hand, there is a second group named the second-generation MEAs, which comprises very efficient optimizers like Strength Pareto Evolutionary Algorithm [36], Non-dominated Sorting Genetic Algorithm II [39], and Multi-objective Optimization Differential Evolution [30], among others. Basically, the main features that distinguish the second-generation MEAs from the first-generation group is the mechanism of adaptation assignment in terms of dominance and the incorporation of elitism [31; 39]. 5.
Differential Evolution - A Brief Description Differential Evolution (DE) is an optimization technique that belongs to the family of
evolutionary computation, which differs from other evolutionary algorithms in the mutation and recombination schemes [23; 24; 30]. DE executes its mutation operation by adding a weighted difference vector between two individuals to a third individual. Then, the mutated individuals will perform discrete crossover and greedy selection with the corresponding individuals from the last generation to produce the offspring. A classical DE algorithm is presented next [24].
Page 9 of 22
NP is the population size, P is the population of the current generation, and P’ is the population to be formed for the next generation, C[i] is the candidate solution with population index i, C[i][j] is the j-th entry in the solution vector of C[i], r is a random number between 0 and 1, CR is the crossover constant, and F is the weight applied to the random differential (scaling factor). Storn and co-workers [24] have given some simple rules for choosing the key parameters of DE for general applications. Normally, NP should be about 5 to 10 times the dimension (number of parameters in a vector) of the problem. As for F, it lies in the range 0.4 to 1.0. Initially F = 0.5 can be tried, then F and/or NP is increased if the population converges prematurely. These authors proposed various mutation schemes for the generation of new vectors
Page 10 of 22
(candidate solutions) by combining the vectors that are randomly chosen from the current population as shown:
6.
rand/1:
x xk F ( xk xk )
rand/2:
x xk F ( xk xk xk xk )
best/1:
x x best F ( x k x k )
best/2:
x x best F ( x k x k x k x k )
1
2
1
3
2
3
2
2
4
5
3
3
4
5
rand/best/1:
x x k F ( x b est x k x k x k )
rand/best/2:
x x k F ( x b est x k ) F ( x k x k x k x k )
1
1
1
1
1
2
1
2
3
4
Multi-objective Optimization Differential Evolution In this section the Multi-objective Optimization Differential Evolution (MODE)
algorithm proposed by Lobato and Steffen Jr [30] to solve multi-objective optimization problems is presented. This evolutionary strategy differs from other algorithms by the incorporation of two operators to the original DE algorithm, namely the mechanisms of rank ordering [31; 36] and exploration of the neighborhood potential solution candidates [40]. The general structure of the proposed algorithm for MOOP using DE is briefly described in the following [30]. An initial population of size NP is randomly generated. All dominated solutions are removed from the population through the operator Fast Non-Dominated Sorting [31]. This procedure is repeated until each vector is a member of a front. Three parents are selected at random in the population. A child is generated from the three parents (this process continues until NP children are generated). Starting from population P1 of size 2NP, neighbors are generated to each one of the individuals of the population in the following way [40]: χ ( x) [ x D k ( g ) / 2, x D k ( g ) / 2]
(21)
where Dk (g )
k
[U L ]
(22)
R
Dk(g) is a vector in
n
and a function of the generation counter g. R is the number of pseudo
fronts defined by the user and the initial maximum neighborhood size in a population is Dk(0) = [U L], where L and U represent the lower and upper bounds of the variables. The pre-defined number of individuals in each pseudo front is given by [40]:
Page 11 of 22
n k r n k 1
k 2,
, R
(23)
where nk is the number of individuals in the k-th front and r (< 1) is the reduction rate. For a given population with N individuals, nk can be calculated as: nk N
1 r 1 r
r
R
k 1
(24)
According to Hu and co-workers [40], if r < 1, the number of individuals in the first pseudo front is the highest and each pseudo front has an exponentially reducing number of solutions, this emphasizing a local search. On the contrary, a greater value for r results in more solutions in the last pseudo front and hence emphasizes the global search. In this way, the neighbors generated are classified according to the dominance criterion and only the non-dominated neighbors (P2) will be put together with P1 to form P3. The population P3 is then classified according to the dominance criterion. If the number of individuals of the population P3 is larger than a number defined by the user, its truncated according to the criterion of the Crowding Distance [31]. The crowding distance describes the density of solutions surrounding a vector. To compute the crowding distance for a set of population members the vectors are sorted according to their objective function values for each objective function. To the vectors with the smallest or largest values an infinite crowding distance (or an arbitrarily large number for practical purposes) is assigned. For all other vectors the crowding distance is calculated according to: m 1
d is t x i
j0
j ,i 1
f
j ,i 1
j ,m ax
f
j , m in
f f
where fj corresponds to the j-th objective function and m equals the number of objective functions. 7.
Methodology As mentioned by Ledzewicz and Schattler [6], Krabs and von Wolfersdorf [9] and by
Kim and Maurer [41], one difficulty in applying the considered optimal control problems in cancer chemotherapy is the occurance of optimal solutions which are seldom or never used in medical practice. In this case, a desired optimal solution by the physician is the bang-bang control consisting of a starting interval with maximal dose of drug (u = umax) followed by an interval zero-therapy (u = umin = 0) until the end of treatment.
Page 12 of 22
In this context, the methodology proposed consists in transforming the original SOCP into a nonlinear optimization problem with constant differential index (equal to one). Let the time interval t [to, … , tf] be discretized using Nt time nodes, ti, such that t0 = to < t1 < … < tNt = tf. In each subinterval t [ti ti+1], i = 1, 2, … , Nt, let the control input be approximated by (bang-bang control) u u i ( u m in o r u m a x )
fo r
ti t ti1
(26)
With the piecewise linear approximation the unknown control input u(t) is replaced by Nt unknown parameters u1, u2, … , uNt. Besides, in this formulation, the localization of each event ti also is unknown and should be calculated, resulting in 2Nt-1 design variables. The resulting nonlinear optimization problem (with differential index equal to one, i.e., bang-bang control) is solved by using the MODE algorithm proposed. 8.
Results and Discussion In order to evaluate the performance of the MODE algorithm, some practical points
regarding the application of the procedure should be emphasized: As mentioned earlier by Pillis and Radunskaya [2], the mathematical model proposed is qualitative and does not focus on a particular tumor type. Objectives: to minimize both the cancer cells concentration and the drug concentration administered to the patient, defined respectively as m in
T dt
m in
udt
To compare the results obtained by using the methodology proposed (MODE), the Non-dominated Sorting Genetic Algorithm (NSGA II) is also considered. The parameters considered by the two evolutionary approaches are presented in the Table 1. It should be emphasized that the preference for the use of the presented parameter setting is due to the variety of complex test-cases that were successfully solved by the MODE and NSGA II algorithms [30; 31]. To discretize the control variable, 15 elements (Nt) were considered. The stopping criterion used: maximum number of objective function evaluations for each evolutionary algorithm (50 + 400 × 50).
Page 13 of 22
It is worth mentioning that the system of differential equations is solved by using the well known Runge-Kutta Method 4-5th order. The parameters used in the test-cases were chosen by Pillis and Radunskaya [2] to guarantee the system stability. A detailed description of the parameters used in the mathematical model and their meaning can be found in [2]. As mentioned by Pillis and Radunskaya [2], the initial immune level is quite small. Besides, the healthy steady-state immune level is given in this model by I = s/d1. We will start at two different immune levels: N(0) = 0.9, T (0) = 0.25 and I(0) = 0.25. This corresponds to a tumor with approximately 0.25 × 1011 cells, or a sphere of radius between 1.8 and 3.9 centimeters. 8.1.
Test-Case 1 This first test-case considers the following parameters [2]: a1 = 0.2, a2 = 0.3, a3 = 0.1, b1
= b2 = 1, c1 = c3 = c4 = 1, c2 = 0.5, d1 = 0.2, r1 = 1.5, r2 = 1, s = 0.33, α = 0.3 and ρ = 0.01. For this set of parameters, Pillis and Radunskaya [2] demonstrated two unstable dead equilibria, one stable coexisting equilibrium (I = 0.44, T = 0.56 and N = 0.44), and a stable tumor free equilibrium (N =1.0, T = 0.0 and I = 1.65) exist. Figure 1 present the Pareto’s Curve (a), the control variable profile (b) and the cells concentration with (c) and without treatment (d). In terms of the Pareto’s Curve, the results obtained by using MODE are equivalent to those obtained by using NSGA II, i.e., the shape of the Pareto’s Curve found is the same for both approaches. It should be observed that the Pareto’s Curve presents the non-dominated solutions, i.e., there exists no other solution that is superior (at least) in the case of cells concentration with tumor, and equal or superior with respect to drug concentration. The profiles presented in Figs. 1(b)-1(d) are defined by reference points C (MODE) and F (NSGA II) (see Table 2). These points are calculated by using the point (3.5,10.5) as a reference of the coordinated axes used to obtain the distance between this point and each one (C or F) of the points (solutions) along the Pareto’s Curve. Thus, the smallest distance obtained was defined as being the choice criterion. In Figure 1(b) its possible to observe the activation (or deactivation) of the control variable, i.e., if the drug is administered or not for each evolutionary strategy. Besides, in both results obtained, the action of the treatment is initially verified in the organism during a larger interval of time in the beginning of the drug adminstration, then the administration is realized in
Page 14 of 22
shorter intervals of time, which preserves the patient’s health from the clinical point of view. The difference observed in these two protocols is due to the difference between the points C and F, i.e., to points calculated by using the MODE and NSGA II, respectively. In Figures 1(c) and 1(d) the cells concentration profiles of normal, immune and tumor are presented with and without the action of the control variable, respectively. In these figures we can visualize the importance of the control strategy used. In this case, if no drug is administered to patient, the cells concentration with tumor increases, and the normal and immune cells concentrations decrease with time. On the other hand, if drug is administered to patient, the cells concentration with tumor decrease (close to zero), and the normal and immune cells concentrations increase with time. In terms of the results obtained by using MODE and NSGA II, the cells concentration profiles are similar, although there is a difference between the control profiles represented by points C and F, as observed in Fig. 1(b). Table 2 presents the best point (A - MODE and D - NSGA II), in terms of the minimization of cells concentration with tumor and the best point (B - MODE and E - NSGA II) in terms of the minimization of drug concentration. The points C (MODE) and F (NSGA II) represents the smallest distance between the reference point (3.5,10.5) and any other point belonging to the Pareto’s Curve, as defined earlier. 8.2.
Test-Case 2 In the second test-case, the same parameters used in the previous test-case are considered,
except the value of s (0.30), which characterizes the steady source rate for immune cells in the absence of a tumor. Figure 2 presents the Pareto’s Curve (a), the control variable profile (b) and the cells concentration with (c) and without treatment (d). As observed in the previous test-case, the shape of the Pareto’s Curve obtained by using MODE and NSGA II are similar. The profiles presented in Figs. 2(b)-2(d) are defined by the reference points C (MODE) and F (NSGA II) (see Table 3), calculated as described above. Figure 2(b) presents the control variable profile considering these points (C and F). The protocols, suggested by the optimal control algorithms, claim that the drug should be administered initially by a long time period, then following basically without drug administration. In Figures 2(c) and 2(d) the cells concentration profiles of normal, immune with tumor are presented with and without action of the control variable, respectively. In these figures we can visualize the effect of the control strategy used. Regarding the cells concentration profiles obtained by using MODE and NSGA II, they are considered as
Page 15 of 22
being similar, although there are differences between the control profiles represented by points C and F (see Fig. 2(b)). Table 3 presents the best point (A - MODE and D - NSGA II) in terms of the minimization of the cells concentration with tumor and the best point (B - MODE and E - NSGA II) regarding the minimization of drug concentration. The points C and F represents the smallest distance between the reference point (0,0) and any other point belonging to the Pareto’s Curve, as defined earlier. In this table we can observe that the difference between each set of points is more significant as compared to the first test-case. 8.3.
Test-Case 3 In the third test-case, the same parameters used in the first test-case are considered,
except the value of ρ (0.02), which characterizes the immune response rate. Figure 3 presents the Pareto’s Curve (a), the control variable profile (b) and the cells concentration with (c) and without treatment (d). As observed in the previous test-cases, the shape of the Pareto’s Curve obtained by using MODE and NSGA II are considered similar, although there are differences between the points (extreme) of the curve (see Fig. 3(a)). The profiles presented in Figs. 3(b)-3(d) are defined by reference points C (MODE) and F (NSGA II) (see Table 4), calculated as described above. Figure 3(b) presents the control variable profile obtained by using MODE and NSGA II considering the points C and F, respectively. As observed in the previous test-cases, each optimal protocol claims that the drug should be administered initially by a long time period, then following basically without drug administration. In practical terms, this result is very interesting, since the tumor cells concentration is minimized while a smaller amount of drug is administered to the patient. In Figures 3(c) and 3(d) the cells concentration profiles of normal, immune with tumor are presented with and without action of the control variable, respectively. As mentioned earlier and in these figures we can visualize the effect of the control strategy used with respect to cells concentration profiles of normal, and immune with tumor during the treatment. Besides, the profiles obtained by using the MODE and NSGA II are considered as being similar. Table 4 presents the best point (A - MODE and D - NSGA II) in terms of the minimization of the cells concentration with tumor and the best point (B - MODE and E - NSGA II) regarding the minimization of drug concentration. The points C (MODE) and F (NSGA II)
Page 16 of 22
represents the smallest distance between the reference point (2.5,0) and any other point belonging to the Pareto’s Curve, as defined earlier. 9.
Conclusions In this contribution, a new multi-objective optimization problem considering both the
minimization of the cancerous cells concentration and the prescribed drug concentration to determine an optimal protocol for drug administration in the treatment of tumors was proposed. To solve this problem, the Multi-objective Optimization Differential Evolution (MODE) algorithm [30] was used. In addition, the results obtained from the application of MODE for three test-cases are compared with those obtained by the Non-dominated Sorting Genetic Algorithm (NSGA II) demonstrating that an optimal protocol can be obtained for the drug administration that satisfactorily fulfills both the objectives. It is important to mention that the results obtained are valid for the parameters used in the differential model, i.e., this work is purely theoretical, however a proof of the concept was demonstrated. Consequently, the use of a mathematical model associated with optimization tools may contribute in the future to the development of optimal protocols to be used in real patients. As mentioned by Pillis and Radunskaya [2], these solutions are not guaranteed to be globally optimal protocols and it may be possible to find other protocols that are just as effective but use less total drug administration. The authors are aware of the difficulties of transferring the numerical results to real world conditions. However, its important to recognize the potential of the present study for the development of new protocols for cancer treatment. It is worth mentioning that the problem formulated in this work is not normally considered in the specialized literature (only the minimization of the cancerous cells concentration is normally proposed). In this context, the formulation of the multi-objective optimization problem and its solution by using the MODE algorithm represents the main contribution of this research effort. As a sequence of this work, we intend to evaluate other models of carcinomas, and compare the obtained optimal protocols with medical protocols used in traditional treatments. Acknowledgments The authors acknowledge the financial support provided by FAPEMIG (Fundação de Amparo a Pesquisa do Estado de Minas Gerais) and CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico), research agencies in Brazil, through the INCT-EIE. References
Page 17 of 22
[1] C. C. Parker, D. P. Dearnaley, Radical radiotherapy for prostate cancer, Cancer Treatment Reviews, 29 (2003), 161. [2] L. G. Pillis, A. Radunskaya A., The dynamics of an optimally controlled tumor model: a case study, Mathematical and Computer Modelling, 37 (2003), 1221. [3] L. G. Pillis, A. E. Radunskaya, A model of tumor growth with optimal control, Argonne National Labs, Technical Report, (2000). [4] L. G. Pillis, A. E. Radunskaya, A mathematical tumor model with immune resistance and drug therapy: an optimal control approach, Jornal of Theoretical Medicine, 3 (2001), 79. [5] K. R. Fister, J. C. Panetta, Optimal control applied to cell-cycle-specific cancer chemotherapy, SIAMJ. Appl. Mathematic, 60 (2000), 1059. [6] U. Ledzewicz, H. Schattler, Optimal bang-bang controls for a 2-compartment model of cancer chemotherapy, Journal of Optimization Theory Applied, 114 (2002a), 609. [7] U. Ledzewicz, H. Schattler, Analysis of a cell-cycle specific model for cancer chemotherapy, Journal of Biological Systems, 10 (2002b), 183. [8] A. Swierniak, U. Ledzewicz, H. Schattler, Optimal control for a class of compartimental models in cancer chemotherapy, International Journal Applied Mathematical Computing Science, 13 (2003), 357. [9] W. Krabs, L. von Wolfersdorf, Two optimal control problems in cancer chemotherapy with drug resistance, Annals of the Academy of Romanian Scientists Series on Mathematics and its Applications, ISSN 2066-6594, 3 (2011). [10] P. Dua, V. Dua, E. N. Pistikopoulos, Optimal delivery of chemotherapeutic agents in cancer, Computers and Chemical Engeneering, 32 (2009), 99. [11] F. Castiglione, B. Piccoli, Cancer immunotherapy, mathematical modelling and optimal control, Journal of Theoretical Biology, 247 (2009), 723. [12] L. G. de Pillis, W. Gu, A. E. Radunskaya, Mixed immunotherapy and chemotherapy of tumours: modelling, applications and biological interpretations, Journal of Theoretical Biology, 238 (2006), 841. [13] A. d’ Onofrio, U. Ledzewicz and H. Maurer, On optimal delivery of combination therapy for tumours, Mathematical Biosciences, 222 (2009), 13. [14] D. Dingli, M. D. Cascino, K. Josic, S. J. Russell, Z. Bajzer, Mathematical modeling of cancer radiovirotherapy, Mathematical Biosciences, 199 (2006), 55.
Page 18 of 22
[15] M. Engelhart, D. Lebiedz, S. Sager, Optimal control for selected cancer chemotherapy ODE models: A view on the potential of optimal schedules and choice of objective function, Mathematical Biosciences, 229 (2011), 123. [16] A. Cappuccio, F. Castiglione, B. Piccoli, Determination of the optimal therapeutic protocols in cancer immunotherapy, Mathematical Biosciences, 209 (2007), 1. [17] A. E. Bryson, Y. C. Ho, Applied Optimal Control, Hemisphere Publishing, Washington (1975). [18] W. F. Feehery, Dynamic Optimization with Path Constraints, PhD thesis, MIT (1998). [19] F. S. Lobato, Hybrid Approach for Dynamic Optimization Problems, M. Sc. Thesis, FEQUI/UFU, Uberlândia, Brazil, (2004) (in portuguese). [20] O. von Stryk, R. Bulirsch, Direct and indirect methods for trajectory optimization, Annals of Operations Research, 37 (1992), 357. [21] O. von Stryk, User’s guide for DIRCOL - a direct collocation method for the numerical solution of optimal control problems, Technische Universitat Darmastadt, Fachgebiet Simulation and Systemoptimierung (SIM), (1999). [22] L. T. Biegler, A. M. Cervantes, A. Wachter, Advances in simultaneous strategies for dynamic process optimization, Chemical Engineering Science, 57 (2002), 575. [23] R. Storn, K. V. Price, Differential evolution: a simple and efficient adaptive scheme for global optimization over continuous spaces, International Computer Science Institute, 12 (1995), 1. [24] R. Storn, K. Price, J. A. Lampinen, Differential Evolution - A Practical Approach to Global Optimization, Springer-Natural Computing Series, (2005). [25] M. D. Kapadi, R. D. Gudi, Optimal control of fed-batch fermentation involving multiple feeds using differential evolution, Process Biochemistry, 39 (2001), 1709. [26] F. S. Lobato, L. C. Oliveira-Lopes, V. V. Murata, V. Steffen Jr., Solution of multi-objective optimal control problems with index fluctuation using differential evolution, 6th Brazilian Conference on Dynamics, Control and Applications-DINCON, (2007). [27] F. Wang, J. Chiou, Optimal control and optimal time location problems of differentialalgebraic systems by differential evolution, Industrial Engineering Chemical Research, 36 (1997), 5348.
Page 19 of 22
[28] A. Chowdhury, A. Ghosh, S. Das, A. Abraham, A hybrid evolutionary direct search technique for solving optimal control problems, 10th International Conference on Hybrid Intelligent Systems (HIS 2010), Atlanta, USA, (2010). [29] F. S. Lobato, A. J. Silva-Neto, V. Steffen Jr, Solution of singular optimal control problems using the improved differential evolution algorithm, Journal of Artificial Intelligence and Soft Computing Research, 1 (2011), 195. [30] F. S. Lobato, V. Steffen Jr, A new multi-objective optimization algorithm based on differential evolution and neighborhood exploring evolution strategy, Journal of Artificial Intelligence and Soft Computing Research, 1 (2011), 259. [31] K. Deb, Multi-objective Optimization using Evolutionary Algorithms, John Wiley & Sons, Chichester-UK, ISBN 0-471-87339-X, (2001). [32] F. Y. Edgeworth, Mathematical Physics, P. Keagan, London, England, (1881). [33] V. Pareto, Cours D’Economie Politique, Vol. I and II, F. Rouge, Lausanne, (1896). [34] F. S. Lobato, Multi-objective Optimization for Engineering System Design, Thesis, Federal University of Uberlândia, Brazil, (2008) (in portuguese). [35] J. D. Schaffer, Multiple Objective Optimization with Vector Evaluated Genetic Algorithms, Ph.D. thesis, Vanderbilt University, (1984). [36] E. Zitzler, L. Thiele, Multi-objective evolutionary algorithms: a comparative case study and the strength pareto approach, IEEE Trans. on Evolutionary Computation, 3 (1999), 257. [37] J. Horn, N. Nafpliotis, Multiobjective optimization using the niched pareto genetic algorithm, IlliGAL Report 93005, Illinois Genetic Algorithms Laboratory, University of Illinois, Champaign-IL, (1993). [38] N. Srinivas, K. Deb, Multiobjective optimization using nondominated sorting in genetic algorithms, Evolutionary Computational, 2 (1994), 221. [39] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA II, IEEE Transactions on Evolutionary Computation, 6 (2002), 182. [40] X. Hu, C. A. C. Coello, Z. Huang, A new multi-objective evolutionary algorithm: neighborhood exploring evolution strategy, Available from: http://www.lania.mx/ ccoello/EMOO, (2006). [41] J. H. R. Kim, H. Maurer, Sensitivity analysis of optimal control problems with bang-bang controls, Proc. of the 42nd IEEE Conf. on Decision and Control, Maui, USA, (2003), 3281.
Page 20 of 22
Figure 1. Pareto’s Curve, Control Variable and Cells Concentration for test-case 1. Figure 2. Pareto’s Curve, Control Variable and Cells Concentration for test-case 2. Figure 3. Pareto’s Curve, Control Variable and Cells Concentration for test-case 3.
Table 1. Parameters considered by the NSGA II and the MODE algorithms. NSGA II [31]
MODE [30]
Population Size
50
50
Selection Strategy
Binary Tournament
DE/rand/1/bin
Crossover Probability
0.8
0.8
Mutation Probability
0.02
-
Number of Generations
400
200
Reduction Rate
-
0.9
Number of Pseudo-Curves
-
10
Table 2. Points of Pareto’s Curve for test-case 1. Point
Cells Concentration
Drug Concentration
with Tumor MODE
NSGA II
A
3.51
49.89
B
12.21
9.92
C
4.67
13.66
D
3.52
47.50
E
11.67
9.67
F
4.60
14.61
Page 21 of 22
Table 3. Points of Pareto’s Curve for test-case 2. Point
Cells Concentration
Drug Concentration
with Tumor MODE
NSGA II
A
5.21
67.81
B
85.04
4.21
C
12.79
21.87
D
5.22
59.43
E
86.87
2.42
F
12.94
22.08
Cells Concentration
Drug Concentration
Table 4. Points of Pareto’s Curve for test-case 3. Point
with Tumor MODE
NSGA II
A
3.46
66.55
B
19.89
4.49
C
6.21
10.23
D
3.46
63.29
E
16.38
5.07
F
8.58
6.73
Page 22 of 22
S0169-2607(15)30460-0 http://dx.doi.org/doi: 10.1016/j.cmpb.2016.04.004 COMM 4130
To appear in:
Computer Methods and Programs in Biomedicine
Received date: Revised date: Accepted date:
25-12-2015 17-3-2016 6-4-2016
Please cite this article as: Fran Sérgio Lobato, Vinicius Silvério Machado, Valder Steffen Jr, Determination of an optimal control strategy for drug administration in tumor treatment using multi-objective optimization differential evolution, Computer Methods and Programs in Biomedicine (2016), http://dx.doi.org/doi: 10.1016/j.cmpb.2016.04.004. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Determination of an Optimal Control Strategy for Drug Administration in Tumor Treatment using Multi-objective Optimization Differential Evolution Fran Sérgio Lobato1, Vinicius Silvério Machado2, Valder Steffen Jr3 1
School of Chemical Engineering, FEQUI,
2
School of Medicine, FAMED,
3
School of Mechanical Engineering, FEMEC
Universidade Federal de Uberlândia, UFU, P.O. Box 593, 38400-902, Uberlândia, Brazil. URL: http://www.feq.ufu.br/ (Fran Sérgio Lobato) Highlights •
To formulate and to solve a multiobjective optimal control problem applied to treatment
of tumors. •
Objectives: to minimize the cancerous cells concentration and the drug concentration.
•
The results demonstrated that an optimal protocol can be obtained considering both the
objectives. Abstract The mathematical modeling of physical and biologic systems represents an interesting alternative to study the behavior of these phenomena. In this context, the development of mathematical models to simulate the dynamic behavior of tumors is configured as an important theme in the current days. Among the advantages resulting from using these models is their application to optimization and inverse problem approaches. Traditionally, the formulated Optimal Control Problem (OCP) has the objective of minimizing the size of tumor cells by the end of the treatment. In this case an important aspect is not considered, namely, the optimal concentrations of drugs may affect the patients’ health significantly. In this sense, the present work has the objective of obtaining an optimal protocol for drug administration to patients with cancer, through the minimization of both the cancerous cells concentration and the prescribed drug
1
[email protected] (Corresponding author) [email protected] 3 [email protected] 2
Page 1 of 22
concentration. The resolution of this multi-objective problem is obtained through the Multiobjective Optimization Differential Evolution (MODE) algorithm. The Pareto’s Curve obtained supplies a set of optimal protocols from which an optimal strategy for drug administration can be chosen, according to a given criterion. Keywords: Mathematical Modeling of Tumors, Differential Evolution Algorithm, Optimal Control Problem, Multi-objective Optimization. 1.
Introduction In the last decades, despite of numerous efforts dedicated to the development of new
technologies for tumors diagnosis and treatment, the incidence of cancer has increased significantly. The ideal treatment for this disease has not yet been found: to destroy tumor cells without affecting normal cells, and inhibit its aggressive manifestation [1]. As mentioned by Pillis and Radunskaya [2], a tumor’s response depends on many factors, including the severity of the disease, the efficiency of the treatment, and the strength of patient’s own immune response, to cite only a few. In this context, the development of mathematical models to represent the tumor growth is configured as an interesting research area. The objective is to find the best protocol for drug administration, although the growth of a cancerous tumor is a very complicated process involving complex interactions [2–4]. In the literature, various mathematical models for cancer treatment have been proposed due to better modeling capabilities and due to the refinement of the techniques that can be used to estimate the necessary control parameters. Besides, various contributions are focused on analyzing the problems including chemotherapy treatment as a mean to deplete cancer cells, immunotherapy as a way to boost the immune system, as well as a combination of the previous mentioned treatments [2–16]. The determination of optimal protocol for drug administration characterizes an Optimal Control Problem (OCP). In the literature, the OCP is solved by using the Optimal Control Theory, with applications in different fields of science and engineering, such as aerospace, process control, robotics, bioengineering, economics, finance, and management science [17–22]. The OCP consists in the determination of the control variable profiles that minimize (maximize) a given performance index as characterized by a differential index, i.e., the number of times required to differentiate the Differential Algebraic-Equation system (DAE) to obtain a system of Ordinary Differential Equations (ODE) [17; 22]. For this problem, when the differential index fluctuation occurs as due to the activation and deactivation of inequality constraints or when the
Page 2 of 22
control variable behaves linearly in the Hamiltonian function, this is called Singular Optimal Control Problem (SOCP) [17–22]. Several numerical methods have been proposed to solve SOCP [17–19], which are classified according to three broad categories: direct optimization methods, Pontryagin’s Maximum Principle (PMP) based methods and HJB based (Hamilton-Jacob-Bellman) methods. Among these methods, the direct approach has been preferentially used in the last several years. This approach uses control parameterization or state and control parameterizations, transforming the original problem into a finite dimensional optimization problem. Recently, algorithms based on heuristic approaches have been used to solve the SOCP due to their success observed in solving general optimization problems. In this context, the Differential Evolution Algorithm (DE), proposed by Storn and Price [23; 24], has been applied successfully. Kapadi and Gudi [25] determined the substrate concentration profiles in a feed-batch reactor with singular arc. Lobato et al. [26] presented a new algorithm for dealing with control optimization problems. The proposed methodology consists in the extension of the DE to problems with multiple objectives, through the incorporation of special mechanisms such as rank ordering and neighborhood potential solution-candidates exploration. This algorithm is applied to determine the switching times (events) and operation time of the lysine fermentation process. Wang and Chiou [27] proposed an algorithm based on DE to determine the optimal control and optimal time location problems of differential-algebraic systems. Chowdhury et al. [28] presented a hybrid evolutionary direct search technique based on DE to solve optimal control problems. Lobato et al. [29] proposed a strategy for the dynamic updating of the population size to reduce the number of objective function evaluations. This strategy is based on the definition of convergence rate to evaluate the homogeneity of the population in the evolutionary process. It has been applied to the solution of singular optimal control problems in chemical and mechanical engineering. Traditionally, the main objective considered in the determination of optimal protocol for drug administration is the minimization of the size of tumor cells at the end of the treatment [2– 4]. In this case, the high drug concentrations in the organism, which affects significantly the patient’s health, is not taken into account. In this contribution, the main goal is to introduce a systematic methodology to minimize simultaneously the cancerous cells concentration and the drug concentration during the treatment by using the Multi-objective Optimization Differential Evolution (MODE) algorithm [30]. The work is organized as follows. Section 2 presents the
Page 3 of 22
mathematical model of tumor growth. In Section 3 the general aspects regarding the formulation and solution of SOCP is presented. A definition of the Multi-objective Optimization Problem (MOP) is presented in Section 4. A review about DE and its extension to deal with multi-criteria optimization are presented in Sections 5 and 6, respectively. In Section 7, the proposed methodology is discussed. The results obtained from the application of the methodology conveyed as based on two test-cases are presented in Section 9. Finally, the conclusions are outlined in Section 9. 2. Mathematical Modeling of Tumor Growth In this work, the mathematical model proposed by Pillis and Radunskaya [2] is adopted. This model has the following characteristics: Immune Response: includes immune cells whose growth may be stimulated by the presence of the tumor and that can destroy tumor cells through a kinetic process; Competition Terms: normal cells and tumor cells compete for available resources, while immune cells and tumor cells compete in a predator-prey fashion; Optimal Control Theory for Chemotherapy: a set of optimal drug therapies is calculated through the minimization of cancerous cells concentration during the treatment, while keeping the normal cells in a condition above a required level; Logistic Growth Law. In the model proposed by Pillis and Radunskaya [2], I denotes the number of immune cells at time t, T represents the number of tumor cells at time t, N stands for the number of normal (host) cells at time t, and u is the control strategy. The source of the immune cells is considered to be outside the system, so that its reasonable to assume a constant influx rate s (steady source rate for immune cells in the absence of a tumor, 0 ≤ s ≤ 0.5). In the absence of any tumor, the cells will die off at a per capita rate dl, resulting in a long-term population size of s/d1 cells. The presence of tumor cells stimulates the immune response as represented by the positive nonlinear growth term for the immune cells ρ IT α T
where ρ and α are positive constants. The tumor cells as well as the normal cells are modeled by a logistic growth law, with parameters ri (r1 may be bigger or smaller than r2, depending on the type of cancer or stage of growth) and bi ( b1 1
1
b2
1)
representing the per capita growth rates
Page 4 of 22
and reciprocal carrying capacities of the two types of cells (i = 1, 2 and 3 identifies those associated with the cells tumor, normal tissue and immune, respectively). Mathematically, this model is described by the following system of differential equation: N r2 N (1 b 2 N ) c 4 T N a 3 u
N (0 ) N
T r1T (1 b 1T ) c 2 I T c 3 T N a 2 u I s
ρ IT α T
c1 I T d 1 I a 1u
o
(1)
T (0 ) To
I (0 ) I o
(2)
(3)
where the parameters ai are the cell death rates (0 ≤ ai ≤ 0.5, with a3 ≤ a1 ≤ a2), ci are the competition rates (ci are considered positive), d1 is the mortality rate, α is related with inverse declivity of immune response (0 ≤ α ≤ 0.5) and ρ is the immune response rate (0 ≤ ρ ≤ 1). In this model, the term T N represents the probability of finding tumor and normal cells, and the term IT describes the probability of finding T and I [2]. It is important to emphasize that, since the model used is qualitative and does not focus on a particular tumor type, its not immediately apparent how to measure which growth law is preferable in this context. However, as mentioned by the authors, the choice of growth law does not significantly affect the qualitative behavior of the model. To better understand the dynamics of the phenomenon, Pillis and Radunskaya [2] analyzed the system without and with drug input. For the first case, the authors demonstrated that depending on the values of the parameters, there could be zero, one, two or three equilibria states. In this analysis tumor-free equilibrium, dead equilibria and coexisting equilibria were considered. More details about the corresponding stability analysis can be found in [2]. 3.
Formulation of the Singular Optimal Control Problem Mathematically, the SOCP can be formulated as follows [17–19]: tf
m in J ( z ( t t ) , t f
u ( t ),t f
L ( z, u , t)dt
(4)
t0
where z is the vector of state variables and u is the vector of control variables. Ψ and L are the first and second terms of the performance index, respectively. The objective is subject to the implicit Differential-Algebraic Equations (DAE) system as given by: f (z, z, u, t) 0 g (z, u, t) 0
(5)
(6)
Page 5 of 22
p (u , t ) 0
q(z, u, t)
(7) 0
ttf
(8)
with consistent initial conditions given by: φ ( z (t0 ), z (t0 ), u (t0 ), t0 ) 0
where J(.), L(.),
.
; f (.),
(9) φ ( .)
mz
;
z
mz
;
u
mu
;
g
mg
;
p
m
p
and
q
mq
.
According to the Optimal Control Theory [17; 18], the solution of the SOCP, defined by Eq.(4) to Eq.(9), is satisfied by the co-state equations and the stationary condition given, respectively, by: λ
T
H u
H
λ (t f )=
z
0
z
(10) ttf
(11)
where H is the Hamiltonian function defined by: H L
T
f
(12)
The system formed by Eq.(10) to Eq.(12) is known as the Euler-Lagrange equations, which are characterized as Boundary Value Problems (BVPs). According to Bryson and Ho [17] and Feehery [18], the main difficulties associated with the SOCP solution are the following: the existence of end-point conditions (or region constraints) implies multipliers and associated complementary conditions that significantly increase the difficulty of solving the BVP by the indirect method; the existence of constraints involving the state variables and the application of slack variables method may originate DAE of higher index, regardless the constraint activation status, even in problems where the number of inequality constraints is equal to the number of control variables; the Lagrange multipliers can be very sensitive to the initial conditions. 3.1.
Singular Arcs The solution of SOCP presents special challenges since it demands the knowledge of the
sequence and the number of constraint activations and deactivations (events) along the trajectory. When the amount of constraints is reduced, its usually possible to determine this sequence by examining the solution of the problem without constraints. However, the presence of a large number of restrictions brings a problem of combinatorial nature [18].
Page 6 of 22
In the events, due to discontinuities in the state and/or the co-state variables, changes in the functional form of the DAE and/or in the trajectories of the control variable in each phase may occur. As a consequence, the differential index of the system can change throughout the solution trajectory, increasing when the inequality becomes active. The existence of sections of fluctuating index leads to different errors in these sections, demanding the application of adjusted numerical strategies for each section. Therefore, its necessary to know previously the moments of activation and deactivation of the restrictions in order to solve the problem adequately. Another difficulty is the presence of singular arcs, where the second derivative matrix of the Hamiltonian with respect to the control is only positive semi-definite [17; 18]. A particular case of great interest is the one that appears when the control variable behaves linearly in the Hamiltonian function. In general, no minimum optimal solution will exist for such problems unless inequality constraints in the state and/or control are specified. If the inequality constraints are linear with respect to the control variable, its reasonable to expect that the minimum solution, if it exists, will always impose that the control variables are located at a point belonging to the border of the viable region of control [17; 19; 22]. Consider the following system of equations: z F1 ( z ) F 2 ( z ) u
z (to ) z o
(13)
with the control variable given by: u m in u u m a x
(14)
The Hamiltonian function is defined as: T
H λ ( F1 ( z ) F 2 ( z ) u )
(15)
For this class of control we have: u m ax u u m in
T
λ F2 0 T
λ F2 0
(16)
T
λ F2 0
where ℑ is the Switching Function [18; 19]. Next section presents basic concepts about multi-objective optimization, as well general aspects about DE and its extension to the multi-objective context. 4.
Multi-objective Optimization Problem
Page 7 of 22
The Multi-objective Optimization Problem (MOP) is very common in different areas such as mathematics, engineering and sciences in general. This problem is different from that of a single-objective optimization since MOP usually has not only one but a set of non-inferior optimal solutions, known as Pareto’s solutions (that form the so-called Pareto’s curve). Besides, the optimality concept for this kind of problem is different from the one used in single optimization problems since in multi-objective optimization its not possible to find a single optimal solution that satisfies all the goals, simultaneously [31]. The notion of optimality for MOP is different from the one used for single optimization problems. The most common idea about multi-objective optimization found in the literature was originally proposed by Francis Ysidro Edgeworth [32] and later generalized by Vilfredo Pareto [33]. This idea can be described as follows [33]: “a solution is optimal if its dominated by no other feasible solution, which means that there exists no other solution that is superior at least in case of one objective function value, and equal or superior with respect to the other values of the objective functions”. This definition leads to finding a set of solutions that is called the Pareto optimal set, whose corresponding elements are called non-dominated or non-inferior solutions. In the literature, there are several methods available in the literature for solving MOPs [31; 34]. These methods follow a preference-based approach, in which a relative preference vector is used to scalarize multiple objectives. As classical optimization methods use a point to point approach, in which a solution found for a given iteration is modified to obtain a new solution, the outcome of using a classical optimization method is a single optimized solution. On the other hand, the evolutionary algorithms can find multiple optimal solutions in one single simulation run due to their population-based search approach, being this way well adapted for multi-objective optimization problems. More details regarding the corresponding methodologies can be found in the literature [31; 34]. Mathematically, the MOP is defined as [31]: m in [ f 1 ( x ) f 2 ( x )
f k ( x )]
(17)
subject to m inequality and l equality constraints g i ( x ) 0 , i 1, .., m
(18)
h j ( x ) 0 , j 1, .., l
(19)
and design space defined according to the decision (or design) variables
Page 8 of 22
U
L
x n x n x n , n 1,
, N.
(20)
where k is the number of objective functions fi:
N
.
In the last decades, several algorithms based on evolutionary mechanisms have been proposed to find approximations to the Pareto optimal solutions. As reported in the literature, the first Multi-objective Evolutionary Algorithm (MEA) was the Vector Evaluated Genetic Algorithm [35]. On one hand there is a first group known as the first-generation, which includes all the early MEAs: Vector Evaluated Genetic Algorithm [35], Non-dominated Pareto Genetic Algorithm [37], Non-dominated Sorting Genetic Algorithm [38]. On the other hand, there is a second group named the second-generation MEAs, which comprises very efficient optimizers like Strength Pareto Evolutionary Algorithm [36], Non-dominated Sorting Genetic Algorithm II [39], and Multi-objective Optimization Differential Evolution [30], among others. Basically, the main features that distinguish the second-generation MEAs from the first-generation group is the mechanism of adaptation assignment in terms of dominance and the incorporation of elitism [31; 39]. 5.
Differential Evolution - A Brief Description Differential Evolution (DE) is an optimization technique that belongs to the family of
evolutionary computation, which differs from other evolutionary algorithms in the mutation and recombination schemes [23; 24; 30]. DE executes its mutation operation by adding a weighted difference vector between two individuals to a third individual. Then, the mutated individuals will perform discrete crossover and greedy selection with the corresponding individuals from the last generation to produce the offspring. A classical DE algorithm is presented next [24].
Page 9 of 22
NP is the population size, P is the population of the current generation, and P’ is the population to be formed for the next generation, C[i] is the candidate solution with population index i, C[i][j] is the j-th entry in the solution vector of C[i], r is a random number between 0 and 1, CR is the crossover constant, and F is the weight applied to the random differential (scaling factor). Storn and co-workers [24] have given some simple rules for choosing the key parameters of DE for general applications. Normally, NP should be about 5 to 10 times the dimension (number of parameters in a vector) of the problem. As for F, it lies in the range 0.4 to 1.0. Initially F = 0.5 can be tried, then F and/or NP is increased if the population converges prematurely. These authors proposed various mutation schemes for the generation of new vectors
Page 10 of 22
(candidate solutions) by combining the vectors that are randomly chosen from the current population as shown:
6.
rand/1:
x xk F ( xk xk )
rand/2:
x xk F ( xk xk xk xk )
best/1:
x x best F ( x k x k )
best/2:
x x best F ( x k x k x k x k )
1
2
1
3
2
3
2
2
4
5
3
3
4
5
rand/best/1:
x x k F ( x b est x k x k x k )
rand/best/2:
x x k F ( x b est x k ) F ( x k x k x k x k )
1
1
1
1
1
2
1
2
3
4
Multi-objective Optimization Differential Evolution In this section the Multi-objective Optimization Differential Evolution (MODE)
algorithm proposed by Lobato and Steffen Jr [30] to solve multi-objective optimization problems is presented. This evolutionary strategy differs from other algorithms by the incorporation of two operators to the original DE algorithm, namely the mechanisms of rank ordering [31; 36] and exploration of the neighborhood potential solution candidates [40]. The general structure of the proposed algorithm for MOOP using DE is briefly described in the following [30]. An initial population of size NP is randomly generated. All dominated solutions are removed from the population through the operator Fast Non-Dominated Sorting [31]. This procedure is repeated until each vector is a member of a front. Three parents are selected at random in the population. A child is generated from the three parents (this process continues until NP children are generated). Starting from population P1 of size 2NP, neighbors are generated to each one of the individuals of the population in the following way [40]: χ ( x) [ x D k ( g ) / 2, x D k ( g ) / 2]
(21)
where Dk (g )
k
[U L ]
(22)
R
Dk(g) is a vector in
n
and a function of the generation counter g. R is the number of pseudo
fronts defined by the user and the initial maximum neighborhood size in a population is Dk(0) = [U L], where L and U represent the lower and upper bounds of the variables. The pre-defined number of individuals in each pseudo front is given by [40]:
Page 11 of 22
n k r n k 1
k 2,
, R
(23)
where nk is the number of individuals in the k-th front and r (< 1) is the reduction rate. For a given population with N individuals, nk can be calculated as: nk N
1 r 1 r
r
R
k 1
(24)
According to Hu and co-workers [40], if r < 1, the number of individuals in the first pseudo front is the highest and each pseudo front has an exponentially reducing number of solutions, this emphasizing a local search. On the contrary, a greater value for r results in more solutions in the last pseudo front and hence emphasizes the global search. In this way, the neighbors generated are classified according to the dominance criterion and only the non-dominated neighbors (P2) will be put together with P1 to form P3. The population P3 is then classified according to the dominance criterion. If the number of individuals of the population P3 is larger than a number defined by the user, its truncated according to the criterion of the Crowding Distance [31]. The crowding distance describes the density of solutions surrounding a vector. To compute the crowding distance for a set of population members the vectors are sorted according to their objective function values for each objective function. To the vectors with the smallest or largest values an infinite crowding distance (or an arbitrarily large number for practical purposes) is assigned. For all other vectors the crowding distance is calculated according to: m 1
d is t x i
j0
j ,i 1
f
j ,i 1
j ,m ax
f
j , m in
f f
where fj corresponds to the j-th objective function and m equals the number of objective functions. 7.
Methodology As mentioned by Ledzewicz and Schattler [6], Krabs and von Wolfersdorf [9] and by
Kim and Maurer [41], one difficulty in applying the considered optimal control problems in cancer chemotherapy is the occurance of optimal solutions which are seldom or never used in medical practice. In this case, a desired optimal solution by the physician is the bang-bang control consisting of a starting interval with maximal dose of drug (u = umax) followed by an interval zero-therapy (u = umin = 0) until the end of treatment.
Page 12 of 22
In this context, the methodology proposed consists in transforming the original SOCP into a nonlinear optimization problem with constant differential index (equal to one). Let the time interval t [to, … , tf] be discretized using Nt time nodes, ti, such that t0 = to < t1 < … < tNt = tf. In each subinterval t [ti ti+1], i = 1, 2, … , Nt, let the control input be approximated by (bang-bang control) u u i ( u m in o r u m a x )
fo r
ti t ti1
(26)
With the piecewise linear approximation the unknown control input u(t) is replaced by Nt unknown parameters u1, u2, … , uNt. Besides, in this formulation, the localization of each event ti also is unknown and should be calculated, resulting in 2Nt-1 design variables. The resulting nonlinear optimization problem (with differential index equal to one, i.e., bang-bang control) is solved by using the MODE algorithm proposed. 8.
Results and Discussion In order to evaluate the performance of the MODE algorithm, some practical points
regarding the application of the procedure should be emphasized: As mentioned earlier by Pillis and Radunskaya [2], the mathematical model proposed is qualitative and does not focus on a particular tumor type. Objectives: to minimize both the cancer cells concentration and the drug concentration administered to the patient, defined respectively as m in
T dt
m in
udt
To compare the results obtained by using the methodology proposed (MODE), the Non-dominated Sorting Genetic Algorithm (NSGA II) is also considered. The parameters considered by the two evolutionary approaches are presented in the Table 1. It should be emphasized that the preference for the use of the presented parameter setting is due to the variety of complex test-cases that were successfully solved by the MODE and NSGA II algorithms [30; 31]. To discretize the control variable, 15 elements (Nt) were considered. The stopping criterion used: maximum number of objective function evaluations for each evolutionary algorithm (50 + 400 × 50).
Page 13 of 22
It is worth mentioning that the system of differential equations is solved by using the well known Runge-Kutta Method 4-5th order. The parameters used in the test-cases were chosen by Pillis and Radunskaya [2] to guarantee the system stability. A detailed description of the parameters used in the mathematical model and their meaning can be found in [2]. As mentioned by Pillis and Radunskaya [2], the initial immune level is quite small. Besides, the healthy steady-state immune level is given in this model by I = s/d1. We will start at two different immune levels: N(0) = 0.9, T (0) = 0.25 and I(0) = 0.25. This corresponds to a tumor with approximately 0.25 × 1011 cells, or a sphere of radius between 1.8 and 3.9 centimeters. 8.1.
Test-Case 1 This first test-case considers the following parameters [2]: a1 = 0.2, a2 = 0.3, a3 = 0.1, b1
= b2 = 1, c1 = c3 = c4 = 1, c2 = 0.5, d1 = 0.2, r1 = 1.5, r2 = 1, s = 0.33, α = 0.3 and ρ = 0.01. For this set of parameters, Pillis and Radunskaya [2] demonstrated two unstable dead equilibria, one stable coexisting equilibrium (I = 0.44, T = 0.56 and N = 0.44), and a stable tumor free equilibrium (N =1.0, T = 0.0 and I = 1.65) exist. Figure 1 present the Pareto’s Curve (a), the control variable profile (b) and the cells concentration with (c) and without treatment (d). In terms of the Pareto’s Curve, the results obtained by using MODE are equivalent to those obtained by using NSGA II, i.e., the shape of the Pareto’s Curve found is the same for both approaches. It should be observed that the Pareto’s Curve presents the non-dominated solutions, i.e., there exists no other solution that is superior (at least) in the case of cells concentration with tumor, and equal or superior with respect to drug concentration. The profiles presented in Figs. 1(b)-1(d) are defined by reference points C (MODE) and F (NSGA II) (see Table 2). These points are calculated by using the point (3.5,10.5) as a reference of the coordinated axes used to obtain the distance between this point and each one (C or F) of the points (solutions) along the Pareto’s Curve. Thus, the smallest distance obtained was defined as being the choice criterion. In Figure 1(b) its possible to observe the activation (or deactivation) of the control variable, i.e., if the drug is administered or not for each evolutionary strategy. Besides, in both results obtained, the action of the treatment is initially verified in the organism during a larger interval of time in the beginning of the drug adminstration, then the administration is realized in
Page 14 of 22
shorter intervals of time, which preserves the patient’s health from the clinical point of view. The difference observed in these two protocols is due to the difference between the points C and F, i.e., to points calculated by using the MODE and NSGA II, respectively. In Figures 1(c) and 1(d) the cells concentration profiles of normal, immune and tumor are presented with and without the action of the control variable, respectively. In these figures we can visualize the importance of the control strategy used. In this case, if no drug is administered to patient, the cells concentration with tumor increases, and the normal and immune cells concentrations decrease with time. On the other hand, if drug is administered to patient, the cells concentration with tumor decrease (close to zero), and the normal and immune cells concentrations increase with time. In terms of the results obtained by using MODE and NSGA II, the cells concentration profiles are similar, although there is a difference between the control profiles represented by points C and F, as observed in Fig. 1(b). Table 2 presents the best point (A - MODE and D - NSGA II), in terms of the minimization of cells concentration with tumor and the best point (B - MODE and E - NSGA II) in terms of the minimization of drug concentration. The points C (MODE) and F (NSGA II) represents the smallest distance between the reference point (3.5,10.5) and any other point belonging to the Pareto’s Curve, as defined earlier. 8.2.
Test-Case 2 In the second test-case, the same parameters used in the previous test-case are considered,
except the value of s (0.30), which characterizes the steady source rate for immune cells in the absence of a tumor. Figure 2 presents the Pareto’s Curve (a), the control variable profile (b) and the cells concentration with (c) and without treatment (d). As observed in the previous test-case, the shape of the Pareto’s Curve obtained by using MODE and NSGA II are similar. The profiles presented in Figs. 2(b)-2(d) are defined by the reference points C (MODE) and F (NSGA II) (see Table 3), calculated as described above. Figure 2(b) presents the control variable profile considering these points (C and F). The protocols, suggested by the optimal control algorithms, claim that the drug should be administered initially by a long time period, then following basically without drug administration. In Figures 2(c) and 2(d) the cells concentration profiles of normal, immune with tumor are presented with and without action of the control variable, respectively. In these figures we can visualize the effect of the control strategy used. Regarding the cells concentration profiles obtained by using MODE and NSGA II, they are considered as
Page 15 of 22
being similar, although there are differences between the control profiles represented by points C and F (see Fig. 2(b)). Table 3 presents the best point (A - MODE and D - NSGA II) in terms of the minimization of the cells concentration with tumor and the best point (B - MODE and E - NSGA II) regarding the minimization of drug concentration. The points C and F represents the smallest distance between the reference point (0,0) and any other point belonging to the Pareto’s Curve, as defined earlier. In this table we can observe that the difference between each set of points is more significant as compared to the first test-case. 8.3.
Test-Case 3 In the third test-case, the same parameters used in the first test-case are considered,
except the value of ρ (0.02), which characterizes the immune response rate. Figure 3 presents the Pareto’s Curve (a), the control variable profile (b) and the cells concentration with (c) and without treatment (d). As observed in the previous test-cases, the shape of the Pareto’s Curve obtained by using MODE and NSGA II are considered similar, although there are differences between the points (extreme) of the curve (see Fig. 3(a)). The profiles presented in Figs. 3(b)-3(d) are defined by reference points C (MODE) and F (NSGA II) (see Table 4), calculated as described above. Figure 3(b) presents the control variable profile obtained by using MODE and NSGA II considering the points C and F, respectively. As observed in the previous test-cases, each optimal protocol claims that the drug should be administered initially by a long time period, then following basically without drug administration. In practical terms, this result is very interesting, since the tumor cells concentration is minimized while a smaller amount of drug is administered to the patient. In Figures 3(c) and 3(d) the cells concentration profiles of normal, immune with tumor are presented with and without action of the control variable, respectively. As mentioned earlier and in these figures we can visualize the effect of the control strategy used with respect to cells concentration profiles of normal, and immune with tumor during the treatment. Besides, the profiles obtained by using the MODE and NSGA II are considered as being similar. Table 4 presents the best point (A - MODE and D - NSGA II) in terms of the minimization of the cells concentration with tumor and the best point (B - MODE and E - NSGA II) regarding the minimization of drug concentration. The points C (MODE) and F (NSGA II)
Page 16 of 22
represents the smallest distance between the reference point (2.5,0) and any other point belonging to the Pareto’s Curve, as defined earlier. 9.
Conclusions In this contribution, a new multi-objective optimization problem considering both the
minimization of the cancerous cells concentration and the prescribed drug concentration to determine an optimal protocol for drug administration in the treatment of tumors was proposed. To solve this problem, the Multi-objective Optimization Differential Evolution (MODE) algorithm [30] was used. In addition, the results obtained from the application of MODE for three test-cases are compared with those obtained by the Non-dominated Sorting Genetic Algorithm (NSGA II) demonstrating that an optimal protocol can be obtained for the drug administration that satisfactorily fulfills both the objectives. It is important to mention that the results obtained are valid for the parameters used in the differential model, i.e., this work is purely theoretical, however a proof of the concept was demonstrated. Consequently, the use of a mathematical model associated with optimization tools may contribute in the future to the development of optimal protocols to be used in real patients. As mentioned by Pillis and Radunskaya [2], these solutions are not guaranteed to be globally optimal protocols and it may be possible to find other protocols that are just as effective but use less total drug administration. The authors are aware of the difficulties of transferring the numerical results to real world conditions. However, its important to recognize the potential of the present study for the development of new protocols for cancer treatment. It is worth mentioning that the problem formulated in this work is not normally considered in the specialized literature (only the minimization of the cancerous cells concentration is normally proposed). In this context, the formulation of the multi-objective optimization problem and its solution by using the MODE algorithm represents the main contribution of this research effort. As a sequence of this work, we intend to evaluate other models of carcinomas, and compare the obtained optimal protocols with medical protocols used in traditional treatments. Acknowledgments The authors acknowledge the financial support provided by FAPEMIG (Fundação de Amparo a Pesquisa do Estado de Minas Gerais) and CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico), research agencies in Brazil, through the INCT-EIE. References
Page 17 of 22
[1] C. C. Parker, D. P. Dearnaley, Radical radiotherapy for prostate cancer, Cancer Treatment Reviews, 29 (2003), 161. [2] L. G. Pillis, A. Radunskaya A., The dynamics of an optimally controlled tumor model: a case study, Mathematical and Computer Modelling, 37 (2003), 1221. [3] L. G. Pillis, A. E. Radunskaya, A model of tumor growth with optimal control, Argonne National Labs, Technical Report, (2000). [4] L. G. Pillis, A. E. Radunskaya, A mathematical tumor model with immune resistance and drug therapy: an optimal control approach, Jornal of Theoretical Medicine, 3 (2001), 79. [5] K. R. Fister, J. C. Panetta, Optimal control applied to cell-cycle-specific cancer chemotherapy, SIAMJ. Appl. Mathematic, 60 (2000), 1059. [6] U. Ledzewicz, H. Schattler, Optimal bang-bang controls for a 2-compartment model of cancer chemotherapy, Journal of Optimization Theory Applied, 114 (2002a), 609. [7] U. Ledzewicz, H. Schattler, Analysis of a cell-cycle specific model for cancer chemotherapy, Journal of Biological Systems, 10 (2002b), 183. [8] A. Swierniak, U. Ledzewicz, H. Schattler, Optimal control for a class of compartimental models in cancer chemotherapy, International Journal Applied Mathematical Computing Science, 13 (2003), 357. [9] W. Krabs, L. von Wolfersdorf, Two optimal control problems in cancer chemotherapy with drug resistance, Annals of the Academy of Romanian Scientists Series on Mathematics and its Applications, ISSN 2066-6594, 3 (2011). [10] P. Dua, V. Dua, E. N. Pistikopoulos, Optimal delivery of chemotherapeutic agents in cancer, Computers and Chemical Engeneering, 32 (2009), 99. [11] F. Castiglione, B. Piccoli, Cancer immunotherapy, mathematical modelling and optimal control, Journal of Theoretical Biology, 247 (2009), 723. [12] L. G. de Pillis, W. Gu, A. E. Radunskaya, Mixed immunotherapy and chemotherapy of tumours: modelling, applications and biological interpretations, Journal of Theoretical Biology, 238 (2006), 841. [13] A. d’ Onofrio, U. Ledzewicz and H. Maurer, On optimal delivery of combination therapy for tumours, Mathematical Biosciences, 222 (2009), 13. [14] D. Dingli, M. D. Cascino, K. Josic, S. J. Russell, Z. Bajzer, Mathematical modeling of cancer radiovirotherapy, Mathematical Biosciences, 199 (2006), 55.
Page 18 of 22
[15] M. Engelhart, D. Lebiedz, S. Sager, Optimal control for selected cancer chemotherapy ODE models: A view on the potential of optimal schedules and choice of objective function, Mathematical Biosciences, 229 (2011), 123. [16] A. Cappuccio, F. Castiglione, B. Piccoli, Determination of the optimal therapeutic protocols in cancer immunotherapy, Mathematical Biosciences, 209 (2007), 1. [17] A. E. Bryson, Y. C. Ho, Applied Optimal Control, Hemisphere Publishing, Washington (1975). [18] W. F. Feehery, Dynamic Optimization with Path Constraints, PhD thesis, MIT (1998). [19] F. S. Lobato, Hybrid Approach for Dynamic Optimization Problems, M. Sc. Thesis, FEQUI/UFU, Uberlândia, Brazil, (2004) (in portuguese). [20] O. von Stryk, R. Bulirsch, Direct and indirect methods for trajectory optimization, Annals of Operations Research, 37 (1992), 357. [21] O. von Stryk, User’s guide for DIRCOL - a direct collocation method for the numerical solution of optimal control problems, Technische Universitat Darmastadt, Fachgebiet Simulation and Systemoptimierung (SIM), (1999). [22] L. T. Biegler, A. M. Cervantes, A. Wachter, Advances in simultaneous strategies for dynamic process optimization, Chemical Engineering Science, 57 (2002), 575. [23] R. Storn, K. V. Price, Differential evolution: a simple and efficient adaptive scheme for global optimization over continuous spaces, International Computer Science Institute, 12 (1995), 1. [24] R. Storn, K. Price, J. A. Lampinen, Differential Evolution - A Practical Approach to Global Optimization, Springer-Natural Computing Series, (2005). [25] M. D. Kapadi, R. D. Gudi, Optimal control of fed-batch fermentation involving multiple feeds using differential evolution, Process Biochemistry, 39 (2001), 1709. [26] F. S. Lobato, L. C. Oliveira-Lopes, V. V. Murata, V. Steffen Jr., Solution of multi-objective optimal control problems with index fluctuation using differential evolution, 6th Brazilian Conference on Dynamics, Control and Applications-DINCON, (2007). [27] F. Wang, J. Chiou, Optimal control and optimal time location problems of differentialalgebraic systems by differential evolution, Industrial Engineering Chemical Research, 36 (1997), 5348.
Page 19 of 22
[28] A. Chowdhury, A. Ghosh, S. Das, A. Abraham, A hybrid evolutionary direct search technique for solving optimal control problems, 10th International Conference on Hybrid Intelligent Systems (HIS 2010), Atlanta, USA, (2010). [29] F. S. Lobato, A. J. Silva-Neto, V. Steffen Jr, Solution of singular optimal control problems using the improved differential evolution algorithm, Journal of Artificial Intelligence and Soft Computing Research, 1 (2011), 195. [30] F. S. Lobato, V. Steffen Jr, A new multi-objective optimization algorithm based on differential evolution and neighborhood exploring evolution strategy, Journal of Artificial Intelligence and Soft Computing Research, 1 (2011), 259. [31] K. Deb, Multi-objective Optimization using Evolutionary Algorithms, John Wiley & Sons, Chichester-UK, ISBN 0-471-87339-X, (2001). [32] F. Y. Edgeworth, Mathematical Physics, P. Keagan, London, England, (1881). [33] V. Pareto, Cours D’Economie Politique, Vol. I and II, F. Rouge, Lausanne, (1896). [34] F. S. Lobato, Multi-objective Optimization for Engineering System Design, Thesis, Federal University of Uberlândia, Brazil, (2008) (in portuguese). [35] J. D. Schaffer, Multiple Objective Optimization with Vector Evaluated Genetic Algorithms, Ph.D. thesis, Vanderbilt University, (1984). [36] E. Zitzler, L. Thiele, Multi-objective evolutionary algorithms: a comparative case study and the strength pareto approach, IEEE Trans. on Evolutionary Computation, 3 (1999), 257. [37] J. Horn, N. Nafpliotis, Multiobjective optimization using the niched pareto genetic algorithm, IlliGAL Report 93005, Illinois Genetic Algorithms Laboratory, University of Illinois, Champaign-IL, (1993). [38] N. Srinivas, K. Deb, Multiobjective optimization using nondominated sorting in genetic algorithms, Evolutionary Computational, 2 (1994), 221. [39] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA II, IEEE Transactions on Evolutionary Computation, 6 (2002), 182. [40] X. Hu, C. A. C. Coello, Z. Huang, A new multi-objective evolutionary algorithm: neighborhood exploring evolution strategy, Available from: http://www.lania.mx/ ccoello/EMOO, (2006). [41] J. H. R. Kim, H. Maurer, Sensitivity analysis of optimal control problems with bang-bang controls, Proc. of the 42nd IEEE Conf. on Decision and Control, Maui, USA, (2003), 3281.
Page 20 of 22
Figure 1. Pareto’s Curve, Control Variable and Cells Concentration for test-case 1. Figure 2. Pareto’s Curve, Control Variable and Cells Concentration for test-case 2. Figure 3. Pareto’s Curve, Control Variable and Cells Concentration for test-case 3.
Table 1. Parameters considered by the NSGA II and the MODE algorithms. NSGA II [31]
MODE [30]
Population Size
50
50
Selection Strategy
Binary Tournament
DE/rand/1/bin
Crossover Probability
0.8
0.8
Mutation Probability
0.02
-
Number of Generations
400
200
Reduction Rate
-
0.9
Number of Pseudo-Curves
-
10
Table 2. Points of Pareto’s Curve for test-case 1. Point
Cells Concentration
Drug Concentration
with Tumor MODE
NSGA II
A
3.51
49.89
B
12.21
9.92
C
4.67
13.66
D
3.52
47.50
E
11.67
9.67
F
4.60
14.61
Page 21 of 22
Table 3. Points of Pareto’s Curve for test-case 2. Point
Cells Concentration
Drug Concentration
with Tumor MODE
NSGA II
A
5.21
67.81
B
85.04
4.21
C
12.79
21.87
D
5.22
59.43
E
86.87
2.42
F
12.94
22.08
Cells Concentration
Drug Concentration
Table 4. Points of Pareto’s Curve for test-case 3. Point
with Tumor MODE
NSGA II
A
3.46
66.55
B
19.89
4.49
C
6.21
10.23
D
3.46
63.29
E
16.38
5.07
F
8.58
6.73
Page 22 of 22