Leukemia

We study the optimal control of leukemic cell population presented in [2], namely the optimal use of a combination of cytostatic and cytotoxic drugs in order to reduce the total cancerous cells. The model from [1] and [3] considers two sub-populations of leukemic cells: the inactive resting cells, and the proliferating cells. Resting cells go into the proliferating phase at constant rate $\beta$ . Proliferating cells either die by apoptosis with rate $\gamma$ or divide during mitosis into two daughter cells. Mitosis has a duration of $2\tau$ and the daughter cells start in the resting phase. We introduce the age variable a for the time spent in the proliferating phase, and note $R(t)$ the resting population and $p(a,t)$ the proliferating population. We consider two drugs: a cytotoxic that induces a death rate $u(t)$ in the second half of the proliferating phase, and a cytostatic $v$ that results in a fraction $k(t)$ of inhibited resting cells that are blocked from entering the proliferating phase. This age-structured model is represented on following figure.

Fig1: Age-structured model

The system follows the dynamics:

$\frac{dR}{dt}(t) = -(1-k(t)) \beta R(t) + 2p(t,2\tau)\\ \frac{dk}{dt}(t) = v(t) (1-k(t)) - \alpha k(t)\\ \frac{\partial p}{\partial t}(t,a) + \frac{\partial p}{\partial a}(t,a) = -(\gamma + \chi_{(\tau,2\tau)}(a)u(t)) p(t,a),\ 2<a <2\tau\\ p(t,0) = (1-k(t)) \beta R(t)$

In order to handle the transport equation for $p$ , we note the total and second half proliferating populations

$P(t) = \int_0^{2\tau} p(t,a)da, \quad P_2(t) = \int_\tau^{2\tau} p(t,a)da$

As detailed in [1] and [2], by the method of characteristics the dynamics can be rewritten as a system of delay differential equations corresponding to a controlled version of Mackey’s model [3].

$\begin{array}{ll} \frac{dR}{dt}(t) = & -(1-k(t)) \beta R(t) + 2(1-k(t-2\tau)\beta R(t-2\tau)e^{-(\gamma 2 \tau + \int_{t-\tau}^\tau u(s) ds)}\\ \frac{dP}{dt}(t) = & -(\gamma P(t) + u(t)P_2(t)) + (1-k(t))\beta R(t) \\ & - (1-k(t-2\tau))\beta R(t-2\tau) e^{-(\gamma 2 \tau + \int_{t-\tau}^\tau u(s) ds)}\\ \frac{dP_2}{dt}(t) = & -(\gamma + u(t)P_2(t)) + (1-k(t-\tau))\beta R(t-\tau)e^{-\gamma \tau}\\ & - (1-k(t-2\tau))\beta R(t-2\tau) e^{-(\gamma 2 \tau + \int_{t-\tau}^\tau u(s) ds)}\\ \frac{dk}{dt}(t) = & v(t) (1-k(t)) - \alpha k(t) \end{array}$

We refer the reader to [2] for the formal analysis of this system, namely the existence of solutions and the introduction of age-weighted populations $\tilde{P}, \tilde{P_2}$ to counteract the horizon effect that would cause to always stop all treatment at the end of the time interval. Taking into account constraints on the maximal cumulative dose of drugs, the resulting optimal control problem writes as

$\left \lbrace \begin{array}{l} Min \ R(T) + \tilde{P}(T)\\ \frac{dR}{dt}(t) = -(1-k(t)) \beta R(t) + \tilde{p}(t,2\tau)\\ \frac{d\tilde{P}}{dt}(t) = \lambda \tilde{P}(t) - u(t) \tilde{P}_2(t) + \tilde{p}(t,0) - \tilde{p}(t,2\tau)\\ \frac{d\tilde{P}_2}{dt}(t) = (\lambda - u(t)) \tilde{P}_2(t) + \tilde{p}(t,\tau) - \tilde{p}(t,2\tau)\\ \frac{dk}{dt}(t) = v(t) (1-k(t)) - \alpha k(t)\\ 0 \le u(t) \le \bar{u}, \quad 0 \le v(t) \le \bar{v}\\ \int_0^T u(t) dt \le \bar{U} \end{array} \right .$

where \lambda is defined by the equation $\lambda+\beta = 2 \beta e^{-(\lambda+\gamma)2\tau}$ and

$\begin{array}{l} \tilde{p}(t,0) = (1-k(t))(\lambda+\beta)R(t)\\ \tilde{p}(t,\tau) = (1-k(t-\tau))(\lambda+\beta)R(t-\tau)e^{\lambda \tau}\\ \tilde{p}(t,2\tau) = (1-k(t-2\tau))(\lambda+\beta)R(t-2\tau)e^{\lambda 2\tau - y(t)}\\ y(t) = \int_{t-\tau}^t u(s)ds \end{array}$

with the border cases

$\begin{array}{ll} \tilde{p}(t,\tau) = p_0(\tau-t)2 e^{-\gamma t - (\lambda+\gamma)\tau} \ if\ t <\tau\\ \tilde{p}(t,2\tau) = p_0(2\tau-t)2 e^{-\gamma t - y(t)} \ if\ t <2\tau \\ y(t) = \int_0^t u(s)ds \ if\ t<\tau\\ \end{array}$

From [2] we take the initial conditions

$R(0)=R_0 = 4\ 10^5,\ k(0)=0,\ y(0)=0,\\ \tilde{P}(0)=\int_0^{2\tau}e^{(\lambda+\gamma)(a-2\tau)}p_0(a)da, \quad \tilde{P}_2(0)=\int_\tau^{2\tau}e^{(\lambda+\gamma)(a-2\tau)}p_0(a)da.\\$

and the parameters

$T=5,\ \tau=1,\ p_0(a)=0.5\ 10^5,\ \alpha=1,\ \beta=2,\ \bar{v}=2,\ \gamma=0.05.$

Numerical simulations are carried out for three cases: $\bar{u}=0.2,\ \bar{U}=2\bar{u}$ (Fig.2); $\bar{u}=1$ , no limit $\bar{U}$ (Fig.3); and $\bar{u}=1,\ \bar{U}=2\bar{u}$ (Fig.4).

Fig2: Low cytotoxic dose $\bar{u}=0.2$ , with cumulative limit. Optimal solution is bang-bang with cytotoxic administered when sub-population $P_2$ is high.

Fig3: High cytotoxic dose $\bar{u}=1$ , without cumulative limit. Optimal cytostatic is bang-bang.

Fig4: High cytotoxic dose $\bar{u}=1$ , with cumulative limit. Both optimal cytostatic and cytotoxic are bang-bang.

References:
[1] M. Adimy, F. Crauste, and A. El Abdllaoui. Discrete maturity-structured model of cell differentiation with applications to acute myelogenous leukemia. J. Biol. Systems, 16(3):395–424, 2008
[2] X. Dupuis. Optimal control of leukemic cell population dynamics. Math. Model. Nat. Phenom., 9(1):4–26, 2014.
[3] M.C. Mackey. Unified hypothesis for the origin of aplastic anemia and periodic hematopoiesis. Blood, 51(5):941–956, 1978.

Links

Contact