Optimization of combination therapy for chronic myeloid leukemia with dosing constraints

Abstract

In this work, we demonstrate a mathematical technique for optimizing combination regimens with constraints. We apply the technique to a mathematical model for treatment of patients with chronic myeloid leukemia. The in-host model includes leukemic cell and immune system dynamics during treatment with tyrosine kinase inhibitors and immunomodulatory compounds. The model is minimal (semi-mechanistic) with just enough detail that all relevant therapeutic effects can be represented. The regimens are optimized to yield the highest possible reduction in disease burden, taking into account dosing constraints and side effect risks due to drug exposure. We compare the following three types of regimens: (1) regimens that are restricted to certain discrete dose levels, which can only change every three months; (2) optimal regimens determined using optimal control; and (3) regimens that are piecewise-constant like the first type of regimen, but are obtained as approximations to the optimal control regimens. All three types of regimens result in similar outcomes, but the last one is easy to compute in addition to being clinically feasible.

Appendix

1.1 Optimal constrained solutions with generalized Hill coefficients in the dynamics

In the model dynamics (1)–(3), we have used Michaelis-Menten terms to represent the drug effects; we have also have used analogous expressions to model some of the interactions between populations. More generally, the drug effects could be represented by a sigmoidal relationship of the form

$$\begin{aligned} \frac{u^n}{UC_{50}^n + u^n} \end{aligned}$$

(5)

with a Hill coefficient $n \ge 1$ (Prinz 2010). In principle, the coefficient n could be drug-dependent and could even depend on which particular drug interaction is being modeled leading to a large variety of situations. Here we give some optimal solutions using the same exponent n in each of the drug-related terms and also for the interaction terms among the populations. Figures 6, 7 and 8 show the optimal year-long constrained solutions for Hill coefficients $n=1.5$, $n=2$ and $n=2.7$ iterated over five years. We chose the value $n=2.7$ based on data for dasatinib (approximate mean $IC_{50}$ and $IC_{90}$ concentrations of CD34+ cells) in the paper by Ishida et al. (2016). The range of $IC_{50}$ values feasible for these data is large and skewed right (from 0.41 to 50 ng/mL with a median of 4.3). If we assume that a representative $IC_{50}$ value is 6.5, and plug in an $IC_{90}$ value of 14.3 for u in Eq. (5), we can numerically determine that $n = 2.7$ makes the expression in Eq. (5) close to 0.9, as desired (that is, the expression gives $90\%$ when we plug in an $IC_{90}$ value for u). The weights in the objective are the same ones that were used before (see the section Constrained Optimization of CML Therapy).

Fig. 6

$(\hbox {n} = 1.5)$ Solutions for the constrained optimization problem [Opt] with piecewise-constant controls in the class $\mathcal {U}_{pc}$, iterated year by year. The graphs show the solutions for years 1 (top), 3 (middle), and 5 (bottom). The dose levels are shown on the left and the resulting populations on the right. A Hill coefficient of $\mathbf {n=1.5}$ was used in the dynamics to describe all pharmacodynamic relations in Eqs. (1)–(3) in the form $\frac{u^n}{UC_{50}^n+u^n}$

Fig. 7

$(\hbox {n} = 2)$ Solutions for the constrained optimization problem [Opt] with piecewise-constant controls in the class $\mathcal {U}_{pc}$, iterated year by year. The graphs show the solutions for years 1 (top), 3 (middle), and 5 (bottom). The dose levels are shown on the left and the resulting populations on the right. A Hill coefficient of $\mathbf {n=2}$ was used in the dynamics to describe all pharmacodynamic relations in Eqs. (1)–(3) in the form $\frac{u^n}{UC_{50}^n+u^n}$

Fig. 8

$(\hbox {n} = 2.7)$ Solutions for the constrained optimization problem [Opt] with piecewise-constant controls in the class $\mathcal {U}_{pc}$, iterated year by year. The graphs show the solutions for years 1 (top), 3 (middle), and 5 (bottom). The dose levels are shown on the left and the resulting populations on the right. A Hill coefficient of $\mathbf {n=2.7}$ was used in the dynamics to describe all pharmacodynamic relations in Eqs. (1)–(3) in the form $\frac{u^n}{UC_{50}^n+u^n}$

Table 3 Cell population values (in cells/$\mu$L) for the constrained optimal solutions for years 1, 3 and 5 for various values of Hill coefficient n

The graphs in Figs. 6, 7, and 8 show the changes in the optimal constrained controls for various values of the exponent n. In the first year, all controls are at full dose irrespective of the exponent, but as the values for P decrease, changes in the controls $u_1$ and $u_2$ arise. Generally, the doses for $u_1$ and $u_2$ are lower in subsequent years and on similar levels as for the case $n=1$ shown in Fig. 2, but the timing of the switching varies. The use of lower doses towards the end of the year and the associated increase in tumor volume then leads to higher doses at the beginning of the subsequent years. The control $u_3$ generally is at full dose throughout all five years, and only for $n=2.7$ does the dose drop, with $u_3=0$ in the last quarter of year 5. Table 3 lists the cell population values of the states after years 1, 3, and 5. Overall, for these scenarios the Michaelis-Menten functional relation with $n=1$ corresponds to a conservative estimate for the effectiveness of therapy as far as the terminal values for the states Q and E are concerned. The final values for P depend more on the overall amounts of drugs used and here there is no clear relation with the exponent n. This mostly is caused by the fact that the effects of one year’s regimen propagate to the next year which almost leads to periodic behavior: if the initial tumor volumes are low, low drug doses are administered, which allows P to increase, which then leads to the administration of higher doses and lowering of P in the following year.

1.2 Optimal solutions for different weights $\gamma _i$ in the objective functional

The choice of weights in the objective is an essential aspect of the optimization and generally one does not know the “correct” values to use. Although data are available on the severity of adverse events for the three classes of therapies included in this work (Imatinib 2001, Dasatinib 2006, Nivolumab 2014), and one can form an opinion of a relative ranking for the three compounds, we do not make any such claims in this work. We simply include hypothetical illustrations for how optimal solutions react to different assessments of the side effects of the drugs. We explore the structure of optimal solutions for another example here, setting

$$\begin{aligned} \gamma _1 = \frac{3}{u_1^{\max }}, \quad \gamma _2 = \frac{5}{u_2^{\max }} \quad \text{ and } \quad \gamma _3 =\frac{3}{u_3^{\max }} \, . \end{aligned}$$

(6)

In this case, we have assigned a higher toxicity penalty to $u_2$ than $u_1$, and we have increased the toxicity penalty for $u_3$. Otherwise, the same initial conditions, dose level limits and parameters as before were used. Figure 9 compares unconstrained optimal controls with optimal constrained regimens in $\mathcal{U}_{pc}$. With these weights, the disease burden in year 1 is still so high that full-dose therapies are optimal. However, once the values of P decrease, because of the higher weights on the controls (i.e., assumed-higher side effects), the optimal dose levels decrease more than in the case considered previously. The dose levels for $u_3$ are now so small that they are all approximated with 0. Indeed, the optimal solutions in the class $\mathcal{U}_{pc}$ often take the value zero during the last quarter of the later years.

Fig. 9

Comparison of the unconstrained optimal control solutions in $\mathcal{U}$ (left) with constrained optimal regimens in $\mathcal{U}_{pc}$ for years 1, 3, and 5 (top to bottom) and weights $\gamma _1 = \frac{3}{u_1^{\max }}$, $\gamma _2 = \frac{5}{u_2^{\max }}$ and $\gamma _3 =\frac{3}{u_3^{\max }}$

We also include several examples for a five-year treatment period. As before, the weights $\gamma _i$ are scaled with their maximum values in the form $\gamma _i = \tilde{\gamma }_i/u_i^{\max }$ for $i = 1, 2, 3$. We have fixed $\tilde{\gamma }_3=20$, but have used different weights on the controls $u_1$ and $u_2$ to give one or the other a higher weight for toxicity. The numerical values of the weights for our runs are given in Table 4. In runs A and B, shown as the top and middle rows of Fig. 10, respectively, the left-hand column shows control $u_1$ with a higher assigned toxicity than $u_2$, with the highest weight on $u_1$ in run A. The optimal controls and their piecewise-constant approximations incorporate these adjustments in the weights and result in lower drug dose levels.

Table 4 Toxicity penalty weights on the controls for the five-year runs shown in Fig. 10. The weight $\tilde{\gamma }_3$ is set to 20 for each of the runs

Fig. 10

Optimal controls in $\mathcal{U}$ (left) and the resulting populations (right) over a five-year treatment period with different weights $\tilde{\gamma }$: $\tilde{\gamma }_1=250$ and $\tilde{\gamma }_2=100$ for run A (top row), $\tilde{\gamma }_1=150$ and $\tilde{\gamma }_2=50$ for run B (middle row), and $\tilde{\gamma }_1=50$ and $\tilde{\gamma }_2=250$ for run C (bottom row). The weight $\tilde{\gamma }_3$ is fixed to 20 for all of these runs

Table 5 Parameter values used in the calculations for the dynamical system shown in Fig. 11. For all parameters not shown here, we used values from Table 1

Fig. 11

Optimal controls in $\mathcal{U}$ (left) and the resulting cell populations (right) over a three-month treatment period with MMR. The weights used are: (top) $\tilde{\gamma }_1=10$, $\tilde{\gamma }_2=3$, and $\tilde{\gamma }_3=1$; and (bottom) $\tilde{\gamma }_1=25$, $\tilde{\gamma }_2=40$, and $\tilde{\gamma }_3=3$

1.3 Maintenance programs and elimination therapies

Model parameters are patient-specific and responses range from long-term survival to almost no effect from the therapy. So far we have considered situations that lie in between these two extremes. Our numerical examples illustrate scenarios that correspond to what could be considered management of the disease over a longer treatment period in a maintenance program. In this section, we include examples in which almost complete eradication of the disease is reached within a short period. Specifically, we include examples in which a reduction of the original BCR-ABL1 levels to less than $0.1 \%$ of the original disease burden [known as major molecular response, or MMR (Baccarani et al. 2009)] is accomplished in a three-month period. Such behavior is only seen for certain sets of parameter values. The parameter values shown in Table 5 were used in the dynamics for the numerical examples in this section. Parameters that are not listed in Table 5 are the same as in Table 1 in the paper.

We used the following weights in the objective:

$$\begin{aligned} \alpha _1 = 5000,&\qquad \beta _1 = 5* \alpha _1, \\ \alpha _2 = 1.67,&\qquad \beta _2 = 0.5* \alpha _2. \end{aligned}$$

We used various weights in the objective for the toxicity on the controls $u_1$, $u_2$, and $u_3$. In the run shown in the top row of Fig. 11 we used $\tilde{\gamma }_1=10$, $\tilde{\gamma }_2=3$, and $\tilde{\gamma }_3=1$. For these weights, the drugs $u_2$ and $u_3$ end up at full dose levels throughout the three-month therapy interval . The drug assigned a higher toxicity penalty (in this case, $u_1$) starts at full dose, but very quickly is reduced to lower dose levels. The values for the leukemic populations after three months are given by $P(3)=0.895$ and $Q(3)=0.028$. This reflects a reduction of P to less than $0.1 \%$ of its initial value (MMR) and a large reduction in the value of Q, as well. In the second run, shown in the bottom row of Fig. 11, the control $u_2$ is assigned a higher toxicity penalty. The weights used are $\tilde{\gamma }_1=25$, $\tilde{\gamma }_2=40$, and $\tilde{\gamma }_3=3$. Although $u_2$ is assigned higher toxicity in this run, overall it is the more effective TKI drug in our model and thus it is only reduced towards the end of therapy. Control $u_1$ follows the same pattern as in the other run, since its assigned toxicity penalty is still quite high. As a result, the leukemic cell reductions are smaller and now the values are $P(3)=1.950$ and $Q(3)=0.029$. The higher P level illustrates the role of the second drug $u_2$ in reducing the leukemic cell populations in our model.

1.4 Summary

In this work, we have computed optimal regimens for combination therapy of CML patients, with constraints on dose levels and the times when they can be changed. Using optimal control, we also computed regimens that are optimal and allowed to change continuously at any time. A comparison showed that the two approaches lead to similar outcomes. The first class of solutions, $\mathcal{U}_{pc}$, is motivated by practical concerns and considers only certain dose levels along with a schedule that only allows the dose to change every three months. The associated optimization problem only has a finite number of choices, but it quickly becomes intractable due to computational complexity. For the second class, $\mathcal{U}$, only the ranges of the doses have limits, and arbitrary changes are allowed within these limits. So this becomes a finite-dimensional continuous-time optimal control problem (Bonnard and Chyba 2003; Bressan and Piccoli 2007; Pontryagin et al. 1964; Schättler and Ledzewicz 2012), which takes only a few minutes to run. We examined a third class of solutions, obtained by approximating optimal solutions $u_{*} \in \mathcal{U}$ by functions in $\mathcal{U}_{pc}$ which are closest to the average of $u_{*}$ on each quarter year. This third class has the advantage of being tractable, as the optimal control problem is numerically straightforward, as is the calculation of the average values of the controls over three-month intervals. The third class is also clinically feasible, as it only has certain fixed dose levels, which can only change every three months. Examples in this work show that these approximations have almost identical cell population results, which is why the third class of solutions is the one we propose should be used for clinical applications.