- Email: [email protected]
Determination of the optimal therapeutic protocols in cancer immunotherapy
Mathematical Biosciences 209 (2007) 1–13 www.elsevier.com/locate/mbs
Determination of the optimal therapeutic protocols in cancer immunotherapy Antonio Cappuccio *, Filippo Castiglione, Benedetto Piccoli Istituto Applicazioni Calcolo – CNR, V.le del Policlinico 137, 00161 Roma, Italy Received 27 April 2006; received in revised form 12 December 2006; accepted 23 February 2007 Available online 6 March 2007
Abstract Cancer immunotherapy aims at eliciting an immune system response against the tumor. However, it is often characterized by toxic side-effects. Limiting the tumor growth and, concurrently, avoiding the toxicity of a drug, is the problem of protocol design. We formulate this question as an optimization problem and derive an algorithm for its solution. Unlike the standard optimal control approach, the algorithm simulates impulse-like drug administrations. It relies on an exact computation of the gradient of the cost function with respect to any protocol by means of the variational equations, that can be solved in parallel with the system. In comparison with previous versions of this method [F. Castiglione, B. Piccoli, Optimal control in a model of dendritic cell transfection cancer immunotherapy, Bull. Math. Biol. 68 (2006) 255; B. Piccoli, F. Castiglione, Optimal vaccine scheduling in cancer immunotherapy, Physica A. 370 (2) (2007) 672], we optimize both the timing and the dosage of each administration and introduce a penalty term to avoid clustering of subsequent injections, a requirement consistent with the clinical practice. In addition, we implement the optimization scheme to simulate the case of multitherapies. The procedure works for any ODE system describing the pharmacokinetics and pharmacodynamics of an arbitrary number of therapeutic agents. In this work, it was tested for a well known model of the tumor– immune system interaction [D. Kirschner, J.C. Panetta, Modeling immunotherapy of tumor–immune interaction, J. Math. Biol. 37 (1998) 235]. Exploring three immunotherapeutic scenarios (CTL therapy, IL-2
*
Corresponding author. Tel.: +39 0688470241; fax: +39 064404306. E-mail addresses: [email protected] (A. Cappuccio), [email protected] (F. Castiglione), [email protected]. it (B. Piccoli). 0025-5564/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2007.02.009
2
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
therapy and combined therapy), we display the stability and efficacy of the optimization method, obtaining protocols that are successful compromises between various clinical requirements. 2007 Elsevier Inc. All rights reserved. Keywords: Optimization; Protocol design; Immunotherapy; Cancer; Drug holidays
1. Introduction Most cancer patients are treated with a combination of surgery, radiotherapy and chemotherapy. Surgery alone hardly suffices for a complete and long term eradication of cancerous cells. On the other hand, radiation and chemotherapy affect both malignant and healthy cells, thus causing severe side-effects. As a rule, the immune system of the patient scarcely supplements such antitumor treatments. Indeed, the host immune system eliminates antigens that recognizes as non-self, while generally ignores what it considers self. The tumor is composed of malignant self-cells and therefore it can elude the immune surveillance. Over the past years, cancer immunotherapy has gained interest as a way to enhance the anticancer immune system response. In clinical practice, this is achieved by using immunostimulants (e.g., interleukins), cell growth factors, vaccines, or transfecting monoclonal antibodies and T cells [4]. Planning an immunotherapeutic treatment is a complex matter. It depends on various aspects like the specific tumor line, its stage of development and the health condition of the patient. The choice of the therapeutic agents is the first step of the problem. In different contexts, a cure may involve the application of one drug only (i.e., monotherapy), or of more drugs (i.e., multitherapy). Once the proper factors are selected, the issue is designing the protocols of administration. The tumor usually follows a rapid growth trajectory in the absence of immune response, and, if inadequate quantities of immunostimulants are applied, the immune system fails to control its growth in the long run; to the contrary, intensive treatments are able in theory to challenge the tumor proliferation, but also induce toxic side-effects. Efficient cures require a balance between the contradictory targets of controlling the tumor growth and of limiting toxicities harmful to the patient. The realization of such a balance is of crucial importance to fully exploit the potential of a therapy. In clinical practice, side-effects are lowered by a proper calibration of the dosage and of the timing for each administration and, on top of that, by inserting temporary interruptions of the treatment (i.e., drug holidays). Given the complexity of the problem, criteria suggesting the optimal dosage and the convenient allocation of drug holidays are mainly driven by empirical methods. Thus, computational methods may provide a useful support to address the problem of selecting the best protocols. A convenient framework to formulate this class of problems is provided by the theory of optimal control. This, in few words, is a tool to sort out the optimal strategies to influence the dynamics of the considered system, i.e., minimizing or maximizing certain performance indexes usually called ‘cost functions’. Originally developed to handle problems of aerospace and cybernetics, optimal control is helping in designing novel applications also in the biomedical research [5–7].
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
3
Recently, an algorithm was introduced to find the optimal schedule consisting of a sequence of inoculations of an immunotherapeutic agent [1]. In the present approach, both the time and the dose of each ’drug’ injection (disregard the nature of the treatment) are viewed as variables to optimize. The cost function is chosen as a combination of four components to match a number of clinical requirements: (i) a final and a (ii) running cost on the tumor dynamics, (iii) a cost including the total amount of drug applied during the treatment and (iv) a cost related to the time-separation between two successive inoculations. The last term corresponds to the medical constraint of inserting drug holidays during the cure to lessen the toxicity due to the accumulation. Furthermore, the method was extended to perform the optimization in the case of multitherapies. The procedure works for any dynamical system describing the pharmacokinetics and pharmacodynamics of an arbitrary number of therapeutic agents. In this study, we applied the method to the Panetta–Kirschner (PK) model [3]. Such model captures some relevant features of the tumor– immune system interaction and, technically relevant, its oscillatory behavior allowed to test the numerical stability and robustness of the algorithm. The results of the simulations indicate that the optimal drug administration satisfies all the required clinical constraints. Moreover, the numerical stability of the optimization method proved to be satisfactory even for oscillatory dynamics exhibited in the PK-model. The paper is organized as follows: in Section 2, we present a brief description of the Panetta– Kirshner model; in Section 3, we introduce the optimization problem and its immunotherapeutic application in Section 4; in Section 5, the results are shown and their discussion is reported in Section 6.
2. The mathematical model and the problem of protocol design The model of Panetta and Kirschner [3] is a well known mathematical description of the tumor– immune system interaction [6,8]. Despite its simplicity, it exhibits rich dynamics that are in qualitative agreement with experimental findings, such as tumor recurrence [9,10] as well as short and long term tumor oscillations, [11,12]. The model represents the reciprocal interactions between the effector cells (E), the tumor cells (T) and the concentration of Interleukin-2 (IL). It consists of the following differential equations: p EI L E_ ¼ cT l2 E þ 1 þ s1 g1 þ I L aET T_ ¼ r2 T ð1 bT Þ ; g2 þ T p ET I_ L ¼ 2 l3 I L þ s2 : g3 þ T
ð1Þ ð2Þ ð3Þ
In brief, the effectors decay at rate l2 and are stimulated by the interaction with the tumor as well as by the presence of Interleukin-2; the tumor follows a logistic growth and is reduced by the effectors; Interleukin-2 is produced when the effectors interact with the tumor and decays at rate l3. The parameters s1 and s2 appearing in Eqs. (1) and (3) incorporate the therapeutic factors. The
4
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
model accounts for three treatment strategies: (I) immunotherapy by CTLs (s150, s2 = 0); (II) immunotherapy by IL-2 (s1 = 0, s250); (III) combined immunotherapy (s150, s250). A detailed stability analysis of the system (1)–(3) can be found in [3]. For completeness, we summarize the main results of that study. In the untreated case (s1 = s2 = 0), tumor eradication is impossible and the presence of oscillations is a typical evolution of the system (Fig. 1). The antigenicity of the tumor (parameter c), and the strength of the immune response (parameter a), determine the period of the oscillations. In case of treatment, significant correlations appear between the antigenicity of the tumor and the type and the intensity of treatment. For CTL therapy, if the intensity of the treatment and the tumor antigenicity are below a certain threshold dependent on the parameters, the tumor size sets to a large value. For larger values of c, the system shows bistability, that is, the system converges to either the large-tumor state or the tumor-free state. For IL-2 immunotherapy, setting s2 below a threshold is substantially equivalent to having no therapy. The tumor size reaches a large steady state, nearly its carrying capacity. In contrast, higher values of s2 are able to completely eradicate the tumor. For IL-2 therapy, the system admits a stable tumor-free steady state: if the treatment is strong enough it can clear the tumor regardless of the values of c, a. However, high amounts of drug are in most clinical cases accompanied by an increasing risk of side-effects that are very dangerous for the patient. Common side-effects of immunotherapeutic treatments include, for instance, nausea, vomit, decreased blood cell counts and fluid retention. These physiological complications pose severe constraints when determining realistic protocols. Reductions of the doses and interruptions of the cure are necessary remedies to avoid toxicities. Therefore, the question of protocol design is to incorporate into a therapy those constraints and concurrently control the tumor growth. 0.004
Tumor 0.00012 0.0001 8e-05 6e-05 4e-05 2e-05 0
tumor cells (arbitrary scale)
0.0035 0.003 0.0025
0
0.002
0.0009 0.0008 0.0007 0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0
CTL
250
500
IL-2
0
time (arbitrary)
250
500
time (arbitrary)
0.0015 0.001 0.0005 0 0
100
200
300
400
500
time (arbitrary)
Fig. 1. Typical tumor oscillations of the PK-model. Parameter values: c = 1.009, l2 = 0.0378, p1 = 0.044, g1 = 0.02, r2 = 0.123, b = 1, a = 0.018, g2 = 104, p2 = 0.9, g3 = 105, l2 = 1.8.
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
5
Techniques of continuous optimal control serves the purpose of finding the best therapeutic strategies and have been already applied to the Panetta–Kirshner model [6]. In clinical applications, the administrative procedures may not take place through continuous infusions (as in [6]), but rather by applying a finite number of relatively short inoculations. In the following, we describe an algorithm to optimize therapies making use of impulse-like drug administrations that is sensible to the above mentioned clinical requirements.
3. The optimization problem Before describing the method applied in this work, we recall a few concepts from the theory. We consider control systems of the following type x_ ¼ F ðx; uÞ ¼ f ðxÞ þ u;
ð4Þ
xð0Þ ¼ x0 :
ð5Þ
where x is the n-state, f(x) a field and u the control function. To fix the ideas, let us assume that the control acts on the first h components of the system, each representing the pharmacokinetics of a specific drug type. A multitherapeutic protocol consists of N administrations of h drug types at some times ti, i = 0, . . . , N 1. This can be represented by a control function uðtÞ ¼
N1 X h X i¼0
ui;k ðt ti Þv½ti ;ti þg ek :
ð6Þ
k¼1
where g is the duration of the administration and ek is the kth vector of the standard Rn basis. In previous works [1], we showed that if g is very short (in comparison with a characteristic time of the model), the control function u in (6) can be approximated by considering the limit g ! 0. In this limit, u can be formally written as uðtÞ ¼
N1 X h X i¼0
~ ui;k dti ðtÞ:
k¼1
where d(Æ) is the Dirac delta function and Z ti þg ~ ui;k ¼ ui;k ðt ti Þdt ti
Fixing a final time s, the space of schedules S is defined as the set S :¼ fti : t0 < t1 < < tN1 < s;
ui;k P 0; i ¼ 0; . . . ; N 1; k ¼ 1; . . . ; hg:
For every element s 2 S, we denote by xs(t) the corresponding trajectory with initial conditions as in (5). Furthermore, we introduce a final cost function U ¼ Uðxs ðsÞ; sÞ: We focus on the following optimization problem (P): minðUðxs ðsÞ; sÞÞ: s2S
6
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
Note that, in general, an analytic dependence of the cost function upon its variables is not available. Our goal is computing the derivatives of U w.r.t. its variables (ti, ui,k), i = 1, . . . , n; k = 1, . . . , h. Then, well known optimization schemes can be implemented to search for the minimum of U. Bearing in mind that the cost function depends on the schedule both explicitly and implicitly through the final point of the trajectory, we can write DUðxs ðsÞ; sÞ oxs ðsÞ oUðxs ðsÞ; sÞ ¼ rx Uðxs ðsÞ; sÞ þ ; Dti oti oti DUxs ðsÞ; sÞ oxs ðsÞ oUðxs ðsÞ; sÞ ¼ rx Uðxs ðsÞ; sÞ þ : Dui;k oui;k oui;k
ð7Þ ð8Þ
and oxous ðsÞ can be given [1,2,13]. Recently, it was proved that an exact computation of the terms oxots ðsÞ i i;k More precisely, the following propositions hold: Proposition 1. Consider an injection at time ti of h drug types with intensity ui,k, k = 1, . . . , h. If, for t > ti, the n-vector vi(t) solves the equation, v_ i ðtÞ ¼ rx f ðxs ðtÞÞ vi ðtÞ; vi ðti Þ ¼ f ðxs ðti ÞÞ f xs ðti Þ þ
h X
! ui;k ek
ð9Þ ð10Þ
k¼1
then vi ðsÞ ¼
oxs ðsÞ : oti
ð11Þ
Proposition 2. Consider an injection at time ti of drug type k with intensity ui,k. If, for t > ti and k = 1, . . . , h, the n-vectors zi,k(t) solve the equations ð12Þ z_ i;k ðtÞ ¼ rx f ðxs ðtÞÞ zi;k ðtÞ; zi;k ðti Þ ¼ ei;k ;
ð13Þ
then zi;k ðsÞ ¼
oxs ðsÞ : oui;k
ð14Þ
Remark 1. Propositions 1 and 2 hold even if the cost function includes an integral part. Indeed, if Z s U ¼ U1 ðxs ðsÞ; sÞ þ U2 ðxs ðtÞÞdt 0
it is sufficient to add a new variable xn+1 satisfying x_ nþ1 ¼ U2 ðxÞ with initial condition xn+1(0) = 0 and then consider the minimization of the extended final cost [2] min ðU1 ðxs ðsÞ; sÞ þ xnþ1 ðsÞÞ: s2S
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
7
4. Applications to cancer immunotherapy In our specific case, the variable x represents the vector (E, T, IL), the field f is given by the equations (1)–(3), while the control function u is given by the vector (s1, 0, 0) for CTL therapy, (0, 0, s2) for IL-2 therapy and (s1, 0, s2) for the combined therapy. The function U to optimize is chosen as a weighted sum of four terms: Uðxs ðsÞ; sÞ ¼
4 X
wh uh :
h¼1
with the following meaning (1) wh, h = 1,. . .,4 are constant weights of the combination; (2) u1 = T(x R s s(s)) is the final tumor mass; (3) u2 ¼ 0 ð½T ðxs ðtÞÞ T max þ Þ2 dt is an integral cost of the tumor mass exceeding a certain max threshold PN1TPh (see [2]), the symbol [Æ]+ indicating the positive part; (4) u3 ¼ i¼0 k¼1 ui;k is a cost associated to the total amount of drug delivered; (5) u4 = u4(t1 t0, . . . , tN1 tN2) is a term penalizing the interval between two following inoculations and is discussed below. The problem (P) turns then into ! 4 X min wh uk : u2S
h¼1
In order to apply the optimization algorithm, we need to compute the derivatives as in (7,8). Propositions 1 and 2 provide the derivatives of u1 and u2 with respect to every element s. P The function u3 ¼ N1 i¼0 ui has derivatives ou3 ¼ 1: oui;k
ð15Þ
for i = 0, . . . , N 1 and k = 0, . . . , h. The term u4 is written as a sum of functions of the time intervals nk = tk+1 tk between two successive administrations u4 ðt1 t0 ; . . . ; tN1 tN2 Þ ¼
N2 X
gðnk Þ:
k¼0
with g(Æ) a decreasing function. The partial derivatives of u4 w.r.t. variables ti are easily calculated. For i = 1, . . . , N 2, we get ou4 og og ¼ ð16Þ oni1 ni1 ¼ti ti1 oni ni ¼tiþ1 ti oti Moreover, we have
8
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
ou4 og ¼ ; on0 n0 ¼t1 t0 ot0 ou4 og ¼ on ot N1
ð17Þ ð18Þ
N2 nN2 ¼tN1 tN2
In the following, we choose gðnk Þ ¼ ednk : where d is a constant whose reciprocal value gives a measure of the drug holiday. In conclusion, we can write an analytic expression of (7,8) DU oxs ðT Þ ou ¼ ðw1 ru1 ðxs ðT Þ; sÞ þ w2 ru2 ðxs ðT Þ; sÞÞ þ w4 4 ; Dti oti oti DU oxs ðT Þ ¼ ðw1 ru1 ðxs ðT Þ; sÞ þ w2 ru2 ðxs ðT Þ; sÞÞ þ w3 : Dui;k oui;k
ð19Þ ð20Þ
Eqs. (19) and (20) along with the relations (11), (14)–(18) give all the necessary information to compute the gradient of the cost function U with respect to any element s of the space of schedules S. Therefore, standard gradient descent algorithms can be applied to find the local minima of U. We implemented the following scheme: Set variables. Fix the time horizon s; the maximum allowed value of the tumor mass during treatment Tmax; the number of vaccine administrations N; the bounds vmin, vmax; an initial value x0 of cells populations and an initial schedule s0 ¼ ðti ; ui;k Þ, i = 0, . . . , N 1 and k = 1, . . . , h; Step 0. Solve the system (1)–(3) with x0 as the initial value via the fourth-order Runge–Kutta integrator generating an approximation of the trajectory xs0 . At the same time, for every i = 0, . . . , N 1 and k = 1, . . . , h solve Eqs. (9) and (12) with initial conditions (10) and (13), respectively; Step 1. Compute the derivatives of U w.r.t. s0, using the Eqs. (19) and (20). to ti and h2 ouoUi;k to ui,k for some small positive Step 2. Update the schedule by adding h1 oU oti parameters h1,h2. GOTO Step 0. An important feature of the algorithm is that Eqs. (9) and (12) can be solved in parallel with the system (1)–(3). This results in a big computational advantage in comparison with a numerical evaluation of the gradient of U. Moreover, the exact computation of the gradient of U assures a greater accuracy. Optimizing the cost functions that are linear combinations of partial targets poses the question of how to calibrate the weights. On the one hand, it seems natural to think at the weights as a measure of the relative importance each target assumes in the optimization problem; on the other hand, this is usually false without further analysis, because the value of the functions uh may vary by orders of magnitude. For this reason, a proper scaling of the weights is necessary. In order to solve this problem we included in the algorithm a kind of ‘bootstrapping’ procedure consisting in running one optimization step with the only purpose of measuring a punctual value of all the partial costs u0i . Then, the weights wi are scaled according to wi ¼ wi =u0i .
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
9
In the following simulations, we set N = 5, s = 500 U, initial conditions x0 = (105, 105, 105), = 105 and d = 0.05. The initial time schedule was randomly generated and equal in all cases T to (106.8, 157.8, 164.9, 354.2, 475.3). The initial dosages were fixed to 103 for CTL therapy and to 101 for IL-2 therapy. In combined therapy, we set lower initial dosages for both treatments: 104 for CTL and 103 for IL-2 dosage. The weights were arbitrarily chosen as follows: w1 = w2 = w4 = 1 in all simulations; w3 = 101 for CTL therapy, w3 = 5Æ103 for IL-2 therapy and w3 = 5 · 104 for the combined therapy. The results of the simulations are reported in next section. max
5. Results We ran three groups of simulations corresponding to cellular immunotherapy, IL-2 immunotherapy and combined therapy, respectively. In view of applying the scaling procedure described in the previous section, we calculated u01 ¼ 1:70 106 ; u02 ¼ 1:69 106 ; u03 ¼ 5 103 ; u04 ¼ 1:27 for CTL therapy; u01 ¼ 2:56 104 ; u02 ¼ 1:63 106 ; u03 ¼ 5 101 ; u04 ¼ 1:27 for IL-2 therapy and u01 ¼ 3:24 105 ; u02 ¼ 1:77 106 ; u03 ¼ 5:5 103 ; u04 ¼ 1:27 for combined therapy. We then i ¼ wi =u0i , with wi specified for each case in the previous section. set w To attain the optimal solution, we ran 1000 optimization steps for CTL and IL-2 immunotherapy and 2000 steps for the combined therapy. No specific arrest condition was considered in the computer routine. During the run of the program, the target function was subjected to substantial reductions, proving the strong impact of the optimization. The total cost in function of the optimization step is displayed for the three cases (Figs. 2(b), 3(b), 4(b)). The resulting optimal objective function values were 0.30 for CTL treatment, 1.46 for IL-2 treatment and 0.025 for combined therapy. The simulations indicated that the combined therapy (Figs. 4(c)) is the most successful in controlling the tumor growth, keeping the tumor size below the threshold Tmax during the overall cure and determining a very low final value. As for the monotherapies, CTL treatment(Fig. 2(c)) was more effective than IL-2 treatment (Fig. 3(c)). It resulted in both a smaller integral cost and a smaller final mass in comparison with IL-2 therapy, but required lower amounts of drug. It is worth noting that the relative impact of the two monotherapies follows the specific choice of the parameters. For instance, the slower decay of the effectors in comparison to Interleukin-2 (i.e., l2 < l3), contributed to promote the greater efficacy of CTL treatment in our simulations. The vaccine cost dynamics during the optimization is also displayed (Figs. 2(b), 3(b), 4(b)). For each type of therapy, the resulting tumor dynamics (dark areas in Figs. 2(c), 3(c), 4(c), 3(c), 4(c)) was compared to the untreated case. The times of injection remained separate and were distributed at rather regular intervals within the treatment period (dark areas in Figs. 2(a), 3(a), 4(a), 3(a), 4(a)). Interestingly, the times of drug inoculation were synchronized with the peaks of the tumor dynamics. Observing that the tumor peaks tend to exceed the value Tmax and that the antitumor effect following an injection is instantaneous, such synchronization seems as a result of the minimization of /2, the integral cost. Intuitively, one would expect the highest doses in correspondence with the maximum peaks in the tumor dynamics. Of note, a comparison between the optimal protocols and the corresponding tumor dynamics displayed that no monotonic relation appeared between the intensity of the peaks
10
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
a
b
500
3 2.5
350
0.26 0.24
2
300
cost
injection time
400
0.28
total cost vaccine cost
450
250
0.22 0.2
1.5
0.18
200 1
150 100
0.16 0.14
0.5
0.12
50 0
0 0
200
400
600
800
0
1000
optimization step
optimization steps
c
0.004
10
CTL dosage
tumor size (adimensional)
0.1 200 400 600 800 1000
-4
0.003 0.002 0.001 0
0
250
500
time -5
5x10
max
T
10
-5
0
100
200
300
400
500
time units (adimensional)
Fig. 2. (a) The CTL schedule dynamics during the optimization. (b) The total cost function during the optimization (left scale) and the CTL dose cost (right scale). (c) The comparison between the tumor dynamics as a result of the optimal CTL treatment (filled) and the untreated case (unfilled). The horizontal line represents the threshold T max = 105. The inset plot shows the optimal CTL treatment.
of the tumor dynamics and the amount of drug delivered at their occurrence. This was rationalized by considering that the drug amount was optimal in relation to the total cost function, and not only in relation to the integral term. The inset plots in Fig. 2(c), 3(c), 4(c), show the optimal protocols for each therapy. In general, the simulations reveal a stable and efficient functioning of the optimization method.
6. Discussion and conclusions We present an optimization algorithm able to find optimal protocols consisting of repeated injections of one or more drugs. The protocols are optimal with respect to a set of medical targets
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
b
500
3
450
0.12
2.6
350
total cost
injection time
400
0.14
total cost vaccine cost
2.8
300 250 200 150
0.1
2.4
0.08
2.2 0.06
2
0.04
1.8
100
0.02
1.6
50
1.4
0 0
200
400
600
800
0
1000
vaccine cost
a
11
0 200 400 600 800 1000
optimization step
optimization steps
IL2 dosage
tumor size (adimensional)
c
10-4
3.5 3 2.5 2 1.5 1 0.5 0
0
250
500
time 5x10-5
T
10
max
-5
0
100
200
300
400
500
time units (adimensional)
Fig. 3. (a) The IL-2 schedule dynamics during the optimization. (b) The total cost function during the optimization (left scale) and the IL-2 dose cost (right scale). (c) The comparison between the tumor dynamics as a result of the optimal IL2 treatment (filled) and the untreated case (unfilled). The horizontal line represents the threshold Tmax = 105. The inset plot shows the optimal IL-2 treatment.
focused on controlling the tumor growth and limiting the impact of possible toxic side-effects. Unlike the standard optimal control approach that makes use of continuous control functions, our method approximates impulse-like drug administration. The algorithm relies on an exact computation of the gradient of the cost function through the use of variational equations that can be solved simultaneously with the system, thus granting both higher accuracy and lower computational cost in comparison to computing numerically the gradient. The numerical schema proved to be stable and robust even for the oscillatory dynamics of the PK-model. In all simulations, the optimal protocol achieved significant reduction of the total cost function. It is worth noting that the presented method results in a non-convex optimization problem. A systematic exploration of the landscape of the cost function may determine even more successful protocols.
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
a
b
500
3 2.5
350 300 250 200
0.016 0.014
2
total cost
injection time
400
0.018
total cost vaccine cost
450
0.012 0.01
1.5
0.008
1
150 100
0.006 0.004
0.5
0.002
50 0
0 0
500
1000
1500
0 500 1000 1500 2000
0
2000
vaccine cost
12
optimization step
optimization steps
c 0.005 CTL
tumor size (adimensional)
0.004
IL2
0.003 10
0.002
-4
0.001 0
0
250
500
time 5x10
-5
T
max
10-5 0
100
200
300
400
500
time units (adimensional)
Fig. 4. (a) The multitherapy schedule dynamics during the optimization. (b) The total cost function during the optimization (left scale) and the multitherapy dose cost (right scale). (c) The comparison between the tumor dynamics as a result of the optimal multitherapy (filled) and the untreated case (unfilled). The horizontal line represents the threshold Tmax = 105. The inset plot shows the optimal multitherapy.
The introduction of a time-separation cost resolved the problem of clustered injection times as found in previous works [1,2]. The optimization showed to keep the time of injections well separated within the time horizon of the treatment, a requirement that is consistent with the clinical practice. Moreover, the method was adapted to optimize multitherapies, that is, treatments composed by more drugs. We considered here a particular case of multitherapy, for which different factors (CTLs and IL-2 in our case) are delivered at the same time. While this restriction can be easily removed from a theoretical point of view, the implementation requires a greater level of sophistication. This fact is due to the possibility that injections corresponding to different drugs cross each other during the course of the optimization, leading to numerical difficulties.
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
13
In general, the applicability of the optimization to practical problems is tightly connected to the adherence of the model with reality; better models would imply more reliable results of the optimization and would allow comparisons to the established protocols. Further developments of this study include the extension of the algorithm to design feedback strategies, whose robustness is of primary importance when passing to clinical applications. Generalizations of the optimization method to handle the case of partial differential equations are an additional progress currently under study.
Acknowledgments This work was partially supported under the EC contract FP6-2004-IST-4, No. 028069 (ImmunoGrid). We thank the anonymous reviewers for useful comments.
References [1] F. Castiglione, B. Piccoli, Optimal control in a model of dendritic cell transfection cancer immunotherapy, Bull. Math. Biol. 68 (2006) 255. [2] B. Piccoli, F. Castiglione, Optimal vaccine scheduling in cancer immunotherapy, Physica A 370 (2) (2007) 672. [3] D. Kirschner, J.C. Panetta, Modeling immunotherapy of tumor-immune interaction, J. Math. Biol. 37 (1998) 235. [4] S. Ben-Efraim, Cancer immunotherapy: hopes and pitfalls: a review, Anticancer Res. 16 (5B) (1996) 3235. [5] R.B. Martin, Optimal control drug scheduling of cancer chemotherapy, Automatica 28 (1992) 1113. [6] T. Burden, J. Ernstberger, K.R. Fister, Optimal control applied to immunotherapy, Disc. Cont. Dyn. Sys. B. 4 (1) (2004) 135. [7] A. Swierniak, U. Ledzewicz, H. Schattler, Optimal control for a class of compartmental models in cancer chemotherapy, Int. J. Appl. Math. Comp. Sci. 13 (3) (2003) 357. [8] L. Preziosi (Ed.), Cancer Modeling and Simulation, Chapman and Hall/CRC Press (UK), London, 2003 (June 26). [9] N. Blumberg, C. Chuang-Stein, J.M. Heal, The relationship of blood transfusion, tumor staging, and cancer recurrence, Transfusion 30 (4) (1990) 291. [10] Yoshihiko Hirao, Eigoro Okajima, Seiichiro Ozono, Shoji Samma, Kenji Sasaki, Tadashi Hiramatsu, Katsuhiro Babaya, Shuji Watanabe, Yoshio Maruyama, A prospective randomized study of prophylaxis of tumor recurrence following transurethral resection of superficial bladder cancer-intravesical thio-TEPA versus oral UFT, Cancer Chemother Pharmacol 30 (1992) S26. [11] Richard A. Gatti, William A. Robinson, Amos S. Deinard, Mark Nesbit, Jeffrey J. McCullough, Mark Ballow, Robert A. Good, Cyclic leukocytosis in chronic myelogenous leukemia: New perspectives on pathogenesis and therapy, Blood 41 (6) (1973) 771. [12] B.J. Kennedy, Cyclic leukocyte oscillations in chronic mylegenous leukemia during hydroxyurea therapy, Blood 35 (6) (1970) 751. [13] Castiglione, Piccoli, Work in Progress.
Determination of the optimal therapeutic protocols in cancer immunotherapy Antonio Cappuccio *, Filippo Castiglione, Benedetto Piccoli Istituto Applicazioni Calcolo – CNR, V.le del Policlinico 137, 00161 Roma, Italy Received 27 April 2006; received in revised form 12 December 2006; accepted 23 February 2007 Available online 6 March 2007
Abstract Cancer immunotherapy aims at eliciting an immune system response against the tumor. However, it is often characterized by toxic side-effects. Limiting the tumor growth and, concurrently, avoiding the toxicity of a drug, is the problem of protocol design. We formulate this question as an optimization problem and derive an algorithm for its solution. Unlike the standard optimal control approach, the algorithm simulates impulse-like drug administrations. It relies on an exact computation of the gradient of the cost function with respect to any protocol by means of the variational equations, that can be solved in parallel with the system. In comparison with previous versions of this method [F. Castiglione, B. Piccoli, Optimal control in a model of dendritic cell transfection cancer immunotherapy, Bull. Math. Biol. 68 (2006) 255; B. Piccoli, F. Castiglione, Optimal vaccine scheduling in cancer immunotherapy, Physica A. 370 (2) (2007) 672], we optimize both the timing and the dosage of each administration and introduce a penalty term to avoid clustering of subsequent injections, a requirement consistent with the clinical practice. In addition, we implement the optimization scheme to simulate the case of multitherapies. The procedure works for any ODE system describing the pharmacokinetics and pharmacodynamics of an arbitrary number of therapeutic agents. In this work, it was tested for a well known model of the tumor– immune system interaction [D. Kirschner, J.C. Panetta, Modeling immunotherapy of tumor–immune interaction, J. Math. Biol. 37 (1998) 235]. Exploring three immunotherapeutic scenarios (CTL therapy, IL-2
*
Corresponding author. Tel.: +39 0688470241; fax: +39 064404306. E-mail addresses: [email protected] (A. Cappuccio), [email protected] (F. Castiglione), [email protected]. it (B. Piccoli). 0025-5564/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2007.02.009
2
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
therapy and combined therapy), we display the stability and efficacy of the optimization method, obtaining protocols that are successful compromises between various clinical requirements. 2007 Elsevier Inc. All rights reserved. Keywords: Optimization; Protocol design; Immunotherapy; Cancer; Drug holidays
1. Introduction Most cancer patients are treated with a combination of surgery, radiotherapy and chemotherapy. Surgery alone hardly suffices for a complete and long term eradication of cancerous cells. On the other hand, radiation and chemotherapy affect both malignant and healthy cells, thus causing severe side-effects. As a rule, the immune system of the patient scarcely supplements such antitumor treatments. Indeed, the host immune system eliminates antigens that recognizes as non-self, while generally ignores what it considers self. The tumor is composed of malignant self-cells and therefore it can elude the immune surveillance. Over the past years, cancer immunotherapy has gained interest as a way to enhance the anticancer immune system response. In clinical practice, this is achieved by using immunostimulants (e.g., interleukins), cell growth factors, vaccines, or transfecting monoclonal antibodies and T cells [4]. Planning an immunotherapeutic treatment is a complex matter. It depends on various aspects like the specific tumor line, its stage of development and the health condition of the patient. The choice of the therapeutic agents is the first step of the problem. In different contexts, a cure may involve the application of one drug only (i.e., monotherapy), or of more drugs (i.e., multitherapy). Once the proper factors are selected, the issue is designing the protocols of administration. The tumor usually follows a rapid growth trajectory in the absence of immune response, and, if inadequate quantities of immunostimulants are applied, the immune system fails to control its growth in the long run; to the contrary, intensive treatments are able in theory to challenge the tumor proliferation, but also induce toxic side-effects. Efficient cures require a balance between the contradictory targets of controlling the tumor growth and of limiting toxicities harmful to the patient. The realization of such a balance is of crucial importance to fully exploit the potential of a therapy. In clinical practice, side-effects are lowered by a proper calibration of the dosage and of the timing for each administration and, on top of that, by inserting temporary interruptions of the treatment (i.e., drug holidays). Given the complexity of the problem, criteria suggesting the optimal dosage and the convenient allocation of drug holidays are mainly driven by empirical methods. Thus, computational methods may provide a useful support to address the problem of selecting the best protocols. A convenient framework to formulate this class of problems is provided by the theory of optimal control. This, in few words, is a tool to sort out the optimal strategies to influence the dynamics of the considered system, i.e., minimizing or maximizing certain performance indexes usually called ‘cost functions’. Originally developed to handle problems of aerospace and cybernetics, optimal control is helping in designing novel applications also in the biomedical research [5–7].
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
3
Recently, an algorithm was introduced to find the optimal schedule consisting of a sequence of inoculations of an immunotherapeutic agent [1]. In the present approach, both the time and the dose of each ’drug’ injection (disregard the nature of the treatment) are viewed as variables to optimize. The cost function is chosen as a combination of four components to match a number of clinical requirements: (i) a final and a (ii) running cost on the tumor dynamics, (iii) a cost including the total amount of drug applied during the treatment and (iv) a cost related to the time-separation between two successive inoculations. The last term corresponds to the medical constraint of inserting drug holidays during the cure to lessen the toxicity due to the accumulation. Furthermore, the method was extended to perform the optimization in the case of multitherapies. The procedure works for any dynamical system describing the pharmacokinetics and pharmacodynamics of an arbitrary number of therapeutic agents. In this study, we applied the method to the Panetta–Kirschner (PK) model [3]. Such model captures some relevant features of the tumor– immune system interaction and, technically relevant, its oscillatory behavior allowed to test the numerical stability and robustness of the algorithm. The results of the simulations indicate that the optimal drug administration satisfies all the required clinical constraints. Moreover, the numerical stability of the optimization method proved to be satisfactory even for oscillatory dynamics exhibited in the PK-model. The paper is organized as follows: in Section 2, we present a brief description of the Panetta– Kirshner model; in Section 3, we introduce the optimization problem and its immunotherapeutic application in Section 4; in Section 5, the results are shown and their discussion is reported in Section 6.
2. The mathematical model and the problem of protocol design The model of Panetta and Kirschner [3] is a well known mathematical description of the tumor– immune system interaction [6,8]. Despite its simplicity, it exhibits rich dynamics that are in qualitative agreement with experimental findings, such as tumor recurrence [9,10] as well as short and long term tumor oscillations, [11,12]. The model represents the reciprocal interactions between the effector cells (E), the tumor cells (T) and the concentration of Interleukin-2 (IL). It consists of the following differential equations: p EI L E_ ¼ cT l2 E þ 1 þ s1 g1 þ I L aET T_ ¼ r2 T ð1 bT Þ ; g2 þ T p ET I_ L ¼ 2 l3 I L þ s2 : g3 þ T
ð1Þ ð2Þ ð3Þ
In brief, the effectors decay at rate l2 and are stimulated by the interaction with the tumor as well as by the presence of Interleukin-2; the tumor follows a logistic growth and is reduced by the effectors; Interleukin-2 is produced when the effectors interact with the tumor and decays at rate l3. The parameters s1 and s2 appearing in Eqs. (1) and (3) incorporate the therapeutic factors. The
4
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
model accounts for three treatment strategies: (I) immunotherapy by CTLs (s150, s2 = 0); (II) immunotherapy by IL-2 (s1 = 0, s250); (III) combined immunotherapy (s150, s250). A detailed stability analysis of the system (1)–(3) can be found in [3]. For completeness, we summarize the main results of that study. In the untreated case (s1 = s2 = 0), tumor eradication is impossible and the presence of oscillations is a typical evolution of the system (Fig. 1). The antigenicity of the tumor (parameter c), and the strength of the immune response (parameter a), determine the period of the oscillations. In case of treatment, significant correlations appear between the antigenicity of the tumor and the type and the intensity of treatment. For CTL therapy, if the intensity of the treatment and the tumor antigenicity are below a certain threshold dependent on the parameters, the tumor size sets to a large value. For larger values of c, the system shows bistability, that is, the system converges to either the large-tumor state or the tumor-free state. For IL-2 immunotherapy, setting s2 below a threshold is substantially equivalent to having no therapy. The tumor size reaches a large steady state, nearly its carrying capacity. In contrast, higher values of s2 are able to completely eradicate the tumor. For IL-2 therapy, the system admits a stable tumor-free steady state: if the treatment is strong enough it can clear the tumor regardless of the values of c, a. However, high amounts of drug are in most clinical cases accompanied by an increasing risk of side-effects that are very dangerous for the patient. Common side-effects of immunotherapeutic treatments include, for instance, nausea, vomit, decreased blood cell counts and fluid retention. These physiological complications pose severe constraints when determining realistic protocols. Reductions of the doses and interruptions of the cure are necessary remedies to avoid toxicities. Therefore, the question of protocol design is to incorporate into a therapy those constraints and concurrently control the tumor growth. 0.004
Tumor 0.00012 0.0001 8e-05 6e-05 4e-05 2e-05 0
tumor cells (arbitrary scale)
0.0035 0.003 0.0025
0
0.002
0.0009 0.0008 0.0007 0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0
CTL
250
500
IL-2
0
time (arbitrary)
250
500
time (arbitrary)
0.0015 0.001 0.0005 0 0
100
200
300
400
500
time (arbitrary)
Fig. 1. Typical tumor oscillations of the PK-model. Parameter values: c = 1.009, l2 = 0.0378, p1 = 0.044, g1 = 0.02, r2 = 0.123, b = 1, a = 0.018, g2 = 104, p2 = 0.9, g3 = 105, l2 = 1.8.
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
5
Techniques of continuous optimal control serves the purpose of finding the best therapeutic strategies and have been already applied to the Panetta–Kirshner model [6]. In clinical applications, the administrative procedures may not take place through continuous infusions (as in [6]), but rather by applying a finite number of relatively short inoculations. In the following, we describe an algorithm to optimize therapies making use of impulse-like drug administrations that is sensible to the above mentioned clinical requirements.
3. The optimization problem Before describing the method applied in this work, we recall a few concepts from the theory. We consider control systems of the following type x_ ¼ F ðx; uÞ ¼ f ðxÞ þ u;
ð4Þ
xð0Þ ¼ x0 :
ð5Þ
where x is the n-state, f(x) a field and u the control function. To fix the ideas, let us assume that the control acts on the first h components of the system, each representing the pharmacokinetics of a specific drug type. A multitherapeutic protocol consists of N administrations of h drug types at some times ti, i = 0, . . . , N 1. This can be represented by a control function uðtÞ ¼
N1 X h X i¼0
ui;k ðt ti Þv½ti ;ti þg ek :
ð6Þ
k¼1
where g is the duration of the administration and ek is the kth vector of the standard Rn basis. In previous works [1], we showed that if g is very short (in comparison with a characteristic time of the model), the control function u in (6) can be approximated by considering the limit g ! 0. In this limit, u can be formally written as uðtÞ ¼
N1 X h X i¼0
~ ui;k dti ðtÞ:
k¼1
where d(Æ) is the Dirac delta function and Z ti þg ~ ui;k ¼ ui;k ðt ti Þdt ti
Fixing a final time s, the space of schedules S is defined as the set S :¼ fti : t0 < t1 < < tN1 < s;
ui;k P 0; i ¼ 0; . . . ; N 1; k ¼ 1; . . . ; hg:
For every element s 2 S, we denote by xs(t) the corresponding trajectory with initial conditions as in (5). Furthermore, we introduce a final cost function U ¼ Uðxs ðsÞ; sÞ: We focus on the following optimization problem (P): minðUðxs ðsÞ; sÞÞ: s2S
6
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
Note that, in general, an analytic dependence of the cost function upon its variables is not available. Our goal is computing the derivatives of U w.r.t. its variables (ti, ui,k), i = 1, . . . , n; k = 1, . . . , h. Then, well known optimization schemes can be implemented to search for the minimum of U. Bearing in mind that the cost function depends on the schedule both explicitly and implicitly through the final point of the trajectory, we can write DUðxs ðsÞ; sÞ oxs ðsÞ oUðxs ðsÞ; sÞ ¼ rx Uðxs ðsÞ; sÞ þ ; Dti oti oti DUxs ðsÞ; sÞ oxs ðsÞ oUðxs ðsÞ; sÞ ¼ rx Uðxs ðsÞ; sÞ þ : Dui;k oui;k oui;k
ð7Þ ð8Þ
and oxous ðsÞ can be given [1,2,13]. Recently, it was proved that an exact computation of the terms oxots ðsÞ i i;k More precisely, the following propositions hold: Proposition 1. Consider an injection at time ti of h drug types with intensity ui,k, k = 1, . . . , h. If, for t > ti, the n-vector vi(t) solves the equation, v_ i ðtÞ ¼ rx f ðxs ðtÞÞ vi ðtÞ; vi ðti Þ ¼ f ðxs ðti ÞÞ f xs ðti Þ þ
h X
! ui;k ek
ð9Þ ð10Þ
k¼1
then vi ðsÞ ¼
oxs ðsÞ : oti
ð11Þ
Proposition 2. Consider an injection at time ti of drug type k with intensity ui,k. If, for t > ti and k = 1, . . . , h, the n-vectors zi,k(t) solve the equations ð12Þ z_ i;k ðtÞ ¼ rx f ðxs ðtÞÞ zi;k ðtÞ; zi;k ðti Þ ¼ ei;k ;
ð13Þ
then zi;k ðsÞ ¼
oxs ðsÞ : oui;k
ð14Þ
Remark 1. Propositions 1 and 2 hold even if the cost function includes an integral part. Indeed, if Z s U ¼ U1 ðxs ðsÞ; sÞ þ U2 ðxs ðtÞÞdt 0
it is sufficient to add a new variable xn+1 satisfying x_ nþ1 ¼ U2 ðxÞ with initial condition xn+1(0) = 0 and then consider the minimization of the extended final cost [2] min ðU1 ðxs ðsÞ; sÞ þ xnþ1 ðsÞÞ: s2S
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
7
4. Applications to cancer immunotherapy In our specific case, the variable x represents the vector (E, T, IL), the field f is given by the equations (1)–(3), while the control function u is given by the vector (s1, 0, 0) for CTL therapy, (0, 0, s2) for IL-2 therapy and (s1, 0, s2) for the combined therapy. The function U to optimize is chosen as a weighted sum of four terms: Uðxs ðsÞ; sÞ ¼
4 X
wh uh :
h¼1
with the following meaning (1) wh, h = 1,. . .,4 are constant weights of the combination; (2) u1 = T(x R s s(s)) is the final tumor mass; (3) u2 ¼ 0 ð½T ðxs ðtÞÞ T max þ Þ2 dt is an integral cost of the tumor mass exceeding a certain max threshold PN1TPh (see [2]), the symbol [Æ]+ indicating the positive part; (4) u3 ¼ i¼0 k¼1 ui;k is a cost associated to the total amount of drug delivered; (5) u4 = u4(t1 t0, . . . , tN1 tN2) is a term penalizing the interval between two following inoculations and is discussed below. The problem (P) turns then into ! 4 X min wh uk : u2S
h¼1
In order to apply the optimization algorithm, we need to compute the derivatives as in (7,8). Propositions 1 and 2 provide the derivatives of u1 and u2 with respect to every element s. P The function u3 ¼ N1 i¼0 ui has derivatives ou3 ¼ 1: oui;k
ð15Þ
for i = 0, . . . , N 1 and k = 0, . . . , h. The term u4 is written as a sum of functions of the time intervals nk = tk+1 tk between two successive administrations u4 ðt1 t0 ; . . . ; tN1 tN2 Þ ¼
N2 X
gðnk Þ:
k¼0
with g(Æ) a decreasing function. The partial derivatives of u4 w.r.t. variables ti are easily calculated. For i = 1, . . . , N 2, we get ou4 og og ¼ ð16Þ oni1 ni1 ¼ti ti1 oni ni ¼tiþ1 ti oti Moreover, we have
8
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
ou4 og ¼ ; on0 n0 ¼t1 t0 ot0 ou4 og ¼ on ot N1
ð17Þ ð18Þ
N2 nN2 ¼tN1 tN2
In the following, we choose gðnk Þ ¼ ednk : where d is a constant whose reciprocal value gives a measure of the drug holiday. In conclusion, we can write an analytic expression of (7,8) DU oxs ðT Þ ou ¼ ðw1 ru1 ðxs ðT Þ; sÞ þ w2 ru2 ðxs ðT Þ; sÞÞ þ w4 4 ; Dti oti oti DU oxs ðT Þ ¼ ðw1 ru1 ðxs ðT Þ; sÞ þ w2 ru2 ðxs ðT Þ; sÞÞ þ w3 : Dui;k oui;k
ð19Þ ð20Þ
Eqs. (19) and (20) along with the relations (11), (14)–(18) give all the necessary information to compute the gradient of the cost function U with respect to any element s of the space of schedules S. Therefore, standard gradient descent algorithms can be applied to find the local minima of U. We implemented the following scheme: Set variables. Fix the time horizon s; the maximum allowed value of the tumor mass during treatment Tmax; the number of vaccine administrations N; the bounds vmin, vmax; an initial value x0 of cells populations and an initial schedule s0 ¼ ðti ; ui;k Þ, i = 0, . . . , N 1 and k = 1, . . . , h; Step 0. Solve the system (1)–(3) with x0 as the initial value via the fourth-order Runge–Kutta integrator generating an approximation of the trajectory xs0 . At the same time, for every i = 0, . . . , N 1 and k = 1, . . . , h solve Eqs. (9) and (12) with initial conditions (10) and (13), respectively; Step 1. Compute the derivatives of U w.r.t. s0, using the Eqs. (19) and (20). to ti and h2 ouoUi;k to ui,k for some small positive Step 2. Update the schedule by adding h1 oU oti parameters h1,h2. GOTO Step 0. An important feature of the algorithm is that Eqs. (9) and (12) can be solved in parallel with the system (1)–(3). This results in a big computational advantage in comparison with a numerical evaluation of the gradient of U. Moreover, the exact computation of the gradient of U assures a greater accuracy. Optimizing the cost functions that are linear combinations of partial targets poses the question of how to calibrate the weights. On the one hand, it seems natural to think at the weights as a measure of the relative importance each target assumes in the optimization problem; on the other hand, this is usually false without further analysis, because the value of the functions uh may vary by orders of magnitude. For this reason, a proper scaling of the weights is necessary. In order to solve this problem we included in the algorithm a kind of ‘bootstrapping’ procedure consisting in running one optimization step with the only purpose of measuring a punctual value of all the partial costs u0i . Then, the weights wi are scaled according to wi ¼ wi =u0i .
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
9
In the following simulations, we set N = 5, s = 500 U, initial conditions x0 = (105, 105, 105), = 105 and d = 0.05. The initial time schedule was randomly generated and equal in all cases T to (106.8, 157.8, 164.9, 354.2, 475.3). The initial dosages were fixed to 103 for CTL therapy and to 101 for IL-2 therapy. In combined therapy, we set lower initial dosages for both treatments: 104 for CTL and 103 for IL-2 dosage. The weights were arbitrarily chosen as follows: w1 = w2 = w4 = 1 in all simulations; w3 = 101 for CTL therapy, w3 = 5Æ103 for IL-2 therapy and w3 = 5 · 104 for the combined therapy. The results of the simulations are reported in next section. max
5. Results We ran three groups of simulations corresponding to cellular immunotherapy, IL-2 immunotherapy and combined therapy, respectively. In view of applying the scaling procedure described in the previous section, we calculated u01 ¼ 1:70 106 ; u02 ¼ 1:69 106 ; u03 ¼ 5 103 ; u04 ¼ 1:27 for CTL therapy; u01 ¼ 2:56 104 ; u02 ¼ 1:63 106 ; u03 ¼ 5 101 ; u04 ¼ 1:27 for IL-2 therapy and u01 ¼ 3:24 105 ; u02 ¼ 1:77 106 ; u03 ¼ 5:5 103 ; u04 ¼ 1:27 for combined therapy. We then i ¼ wi =u0i , with wi specified for each case in the previous section. set w To attain the optimal solution, we ran 1000 optimization steps for CTL and IL-2 immunotherapy and 2000 steps for the combined therapy. No specific arrest condition was considered in the computer routine. During the run of the program, the target function was subjected to substantial reductions, proving the strong impact of the optimization. The total cost in function of the optimization step is displayed for the three cases (Figs. 2(b), 3(b), 4(b)). The resulting optimal objective function values were 0.30 for CTL treatment, 1.46 for IL-2 treatment and 0.025 for combined therapy. The simulations indicated that the combined therapy (Figs. 4(c)) is the most successful in controlling the tumor growth, keeping the tumor size below the threshold Tmax during the overall cure and determining a very low final value. As for the monotherapies, CTL treatment(Fig. 2(c)) was more effective than IL-2 treatment (Fig. 3(c)). It resulted in both a smaller integral cost and a smaller final mass in comparison with IL-2 therapy, but required lower amounts of drug. It is worth noting that the relative impact of the two monotherapies follows the specific choice of the parameters. For instance, the slower decay of the effectors in comparison to Interleukin-2 (i.e., l2 < l3), contributed to promote the greater efficacy of CTL treatment in our simulations. The vaccine cost dynamics during the optimization is also displayed (Figs. 2(b), 3(b), 4(b)). For each type of therapy, the resulting tumor dynamics (dark areas in Figs. 2(c), 3(c), 4(c), 3(c), 4(c)) was compared to the untreated case. The times of injection remained separate and were distributed at rather regular intervals within the treatment period (dark areas in Figs. 2(a), 3(a), 4(a), 3(a), 4(a)). Interestingly, the times of drug inoculation were synchronized with the peaks of the tumor dynamics. Observing that the tumor peaks tend to exceed the value Tmax and that the antitumor effect following an injection is instantaneous, such synchronization seems as a result of the minimization of /2, the integral cost. Intuitively, one would expect the highest doses in correspondence with the maximum peaks in the tumor dynamics. Of note, a comparison between the optimal protocols and the corresponding tumor dynamics displayed that no monotonic relation appeared between the intensity of the peaks
10
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
a
b
500
3 2.5
350
0.26 0.24
2
300
cost
injection time
400
0.28
total cost vaccine cost
450
250
0.22 0.2
1.5
0.18
200 1
150 100
0.16 0.14
0.5
0.12
50 0
0 0
200
400
600
800
0
1000
optimization step
optimization steps
c
0.004
10
CTL dosage
tumor size (adimensional)
0.1 200 400 600 800 1000
-4
0.003 0.002 0.001 0
0
250
500
time -5
5x10
max
T
10
-5
0
100
200
300
400
500
time units (adimensional)
Fig. 2. (a) The CTL schedule dynamics during the optimization. (b) The total cost function during the optimization (left scale) and the CTL dose cost (right scale). (c) The comparison between the tumor dynamics as a result of the optimal CTL treatment (filled) and the untreated case (unfilled). The horizontal line represents the threshold T max = 105. The inset plot shows the optimal CTL treatment.
of the tumor dynamics and the amount of drug delivered at their occurrence. This was rationalized by considering that the drug amount was optimal in relation to the total cost function, and not only in relation to the integral term. The inset plots in Fig. 2(c), 3(c), 4(c), show the optimal protocols for each therapy. In general, the simulations reveal a stable and efficient functioning of the optimization method.
6. Discussion and conclusions We present an optimization algorithm able to find optimal protocols consisting of repeated injections of one or more drugs. The protocols are optimal with respect to a set of medical targets
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
b
500
3
450
0.12
2.6
350
total cost
injection time
400
0.14
total cost vaccine cost
2.8
300 250 200 150
0.1
2.4
0.08
2.2 0.06
2
0.04
1.8
100
0.02
1.6
50
1.4
0 0
200
400
600
800
0
1000
vaccine cost
a
11
0 200 400 600 800 1000
optimization step
optimization steps
IL2 dosage
tumor size (adimensional)
c
10-4
3.5 3 2.5 2 1.5 1 0.5 0
0
250
500
time 5x10-5
T
10
max
-5
0
100
200
300
400
500
time units (adimensional)
Fig. 3. (a) The IL-2 schedule dynamics during the optimization. (b) The total cost function during the optimization (left scale) and the IL-2 dose cost (right scale). (c) The comparison between the tumor dynamics as a result of the optimal IL2 treatment (filled) and the untreated case (unfilled). The horizontal line represents the threshold Tmax = 105. The inset plot shows the optimal IL-2 treatment.
focused on controlling the tumor growth and limiting the impact of possible toxic side-effects. Unlike the standard optimal control approach that makes use of continuous control functions, our method approximates impulse-like drug administration. The algorithm relies on an exact computation of the gradient of the cost function through the use of variational equations that can be solved simultaneously with the system, thus granting both higher accuracy and lower computational cost in comparison to computing numerically the gradient. The numerical schema proved to be stable and robust even for the oscillatory dynamics of the PK-model. In all simulations, the optimal protocol achieved significant reduction of the total cost function. It is worth noting that the presented method results in a non-convex optimization problem. A systematic exploration of the landscape of the cost function may determine even more successful protocols.
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
a
b
500
3 2.5
350 300 250 200
0.016 0.014
2
total cost
injection time
400
0.018
total cost vaccine cost
450
0.012 0.01
1.5
0.008
1
150 100
0.006 0.004
0.5
0.002
50 0
0 0
500
1000
1500
0 500 1000 1500 2000
0
2000
vaccine cost
12
optimization step
optimization steps
c 0.005 CTL
tumor size (adimensional)
0.004
IL2
0.003 10
0.002
-4
0.001 0
0
250
500
time 5x10
-5
T
max
10-5 0
100
200
300
400
500
time units (adimensional)
Fig. 4. (a) The multitherapy schedule dynamics during the optimization. (b) The total cost function during the optimization (left scale) and the multitherapy dose cost (right scale). (c) The comparison between the tumor dynamics as a result of the optimal multitherapy (filled) and the untreated case (unfilled). The horizontal line represents the threshold Tmax = 105. The inset plot shows the optimal multitherapy.
The introduction of a time-separation cost resolved the problem of clustered injection times as found in previous works [1,2]. The optimization showed to keep the time of injections well separated within the time horizon of the treatment, a requirement that is consistent with the clinical practice. Moreover, the method was adapted to optimize multitherapies, that is, treatments composed by more drugs. We considered here a particular case of multitherapy, for which different factors (CTLs and IL-2 in our case) are delivered at the same time. While this restriction can be easily removed from a theoretical point of view, the implementation requires a greater level of sophistication. This fact is due to the possibility that injections corresponding to different drugs cross each other during the course of the optimization, leading to numerical difficulties.
A. Cappuccio et al. / Mathematical Biosciences 209 (2007) 1–13
13
In general, the applicability of the optimization to practical problems is tightly connected to the adherence of the model with reality; better models would imply more reliable results of the optimization and would allow comparisons to the established protocols. Further developments of this study include the extension of the algorithm to design feedback strategies, whose robustness is of primary importance when passing to clinical applications. Generalizations of the optimization method to handle the case of partial differential equations are an additional progress currently under study.
Acknowledgments This work was partially supported under the EC contract FP6-2004-IST-4, No. 028069 (ImmunoGrid). We thank the anonymous reviewers for useful comments.
References [1] F. Castiglione, B. Piccoli, Optimal control in a model of dendritic cell transfection cancer immunotherapy, Bull. Math. Biol. 68 (2006) 255. [2] B. Piccoli, F. Castiglione, Optimal vaccine scheduling in cancer immunotherapy, Physica A 370 (2) (2007) 672. [3] D. Kirschner, J.C. Panetta, Modeling immunotherapy of tumor-immune interaction, J. Math. Biol. 37 (1998) 235. [4] S. Ben-Efraim, Cancer immunotherapy: hopes and pitfalls: a review, Anticancer Res. 16 (5B) (1996) 3235. [5] R.B. Martin, Optimal control drug scheduling of cancer chemotherapy, Automatica 28 (1992) 1113. [6] T. Burden, J. Ernstberger, K.R. Fister, Optimal control applied to immunotherapy, Disc. Cont. Dyn. Sys. B. 4 (1) (2004) 135. [7] A. Swierniak, U. Ledzewicz, H. Schattler, Optimal control for a class of compartmental models in cancer chemotherapy, Int. J. Appl. Math. Comp. Sci. 13 (3) (2003) 357. [8] L. Preziosi (Ed.), Cancer Modeling and Simulation, Chapman and Hall/CRC Press (UK), London, 2003 (June 26). [9] N. Blumberg, C. Chuang-Stein, J.M. Heal, The relationship of blood transfusion, tumor staging, and cancer recurrence, Transfusion 30 (4) (1990) 291. [10] Yoshihiko Hirao, Eigoro Okajima, Seiichiro Ozono, Shoji Samma, Kenji Sasaki, Tadashi Hiramatsu, Katsuhiro Babaya, Shuji Watanabe, Yoshio Maruyama, A prospective randomized study of prophylaxis of tumor recurrence following transurethral resection of superficial bladder cancer-intravesical thio-TEPA versus oral UFT, Cancer Chemother Pharmacol 30 (1992) S26. [11] Richard A. Gatti, William A. Robinson, Amos S. Deinard, Mark Nesbit, Jeffrey J. McCullough, Mark Ballow, Robert A. Good, Cyclic leukocytosis in chronic myelogenous leukemia: New perspectives on pathogenesis and therapy, Blood 41 (6) (1973) 771. [12] B.J. Kennedy, Cyclic leukocyte oscillations in chronic mylegenous leukemia during hydroxyurea therapy, Blood 35 (6) (1970) 751. [13] Castiglione, Piccoli, Work in Progress.