# American Institute of Mathematical Sciences

May  2019, 24(5): 2383-2405. doi: 10.3934/dcdsb.2019100

## Optimal control for cancer chemotherapy under tumor heterogeneity with Michealis-Menten pharmacodynamics

 1 Department of Mechanical and Aerospace Engineering, University of Texas at Arlington, Arlington, TX, 76010, USA 2 Dept. of Electrical and Systems Engineering, Washington University, St. Louis, Mo, 63130, USA

* Corresponding author: Heinz Schättler

Received  January 2018 Revised  January 2019 Published  March 2019

We analyze mathematical models for cancer chemotherapy under tumor heterogeneity as optimal control problems. Tumor heterogeneity is incorporated into the model through a potentially large number of states which represent sub-populations of tumor cells with different chemotherapeutic sensitivities. In this paper, a Michaelis-Menten type functional form is used to model the pharmacodynamic effects of the drug concentrations. In the objective, a weighted average of the tumor volume and the total amounts of drugs administered (taken as an indirect measurement for the side effects) is minimized. Mathematically, incorporating a Michaelis-Menten form creates a nonlinear structure in the controls with partial convexity properties in the Hamiltonian function for the optimal control problem. As a result, optimal controls are continuous and this fact can be utilized to set up an efficient numerical computation of extremals. Second order Jacobi type necessary and sufficient conditions for optimality are formulated that allow to verify the strong local optimality of numerically computed extremal controlled trajectories. Examples which illustrate the cases of both strong locally optimal and non-optimal extremals are given.

Citation: Shuo Wang, Heinz Schättler. Optimal control for cancer chemotherapy under tumor heterogeneity with Michealis-Menten pharmacodynamics. Discrete & Continuous Dynamical Systems - B, 2019, 24 (5) : 2383-2405. doi: 10.3934/dcdsb.2019100
##### References:

show all references

##### References:
A numerically computed extremal control (top row, left) for the effectiveness function $\varphi(x) = \frac{2}{1+x^2}$ (top row, right), weights $\alpha = (1, \ldots, 1)$, $\beta = 500\alpha$ and $\gamma = 5000$, and therapy horizon $T = 10$. The corresponding evolution of the individual cell populations is shown in the middle row (left) along with the evolution of the average of the populations (right). A comparison of the initial density $n(0, x)\equiv \frac{200}{21}$ and the terminal density $n(10, x)$ shown as red curve is given in the bottom row on the left. This control is a strong local minimum as the entries of the solution $S$ to the Riccati differential equation (13) along the extremal controlled trajectory (shown in the bottom row, right) remain finite over the full interval $[0, T]$
A numerically computed extremal control (top row, left) for the effectiveness function $\varphi_1$ and $\varphi_2$ given in (22), (top row, right), $\alpha = (1, \ldots, 1)$, $\beta = 50 \alpha$, $\gamma = (600, 1400)$ and time horizon $T = 10$. The corresponding evolution of the individual cell populations is shown in the middle row (left) along with the evolution of the average of the populations (right). A comparison of the initial density $n(0, x)\equiv \frac{200}{21}$ and the terminal density $n(10, x)$ shown as red curve is given in the bottom row on the left. This control is a strong local minimum as the entries of the solution $S$ to the Riccati differential equation (13) along the extremal controlled trajectory (shown in the bottom row, right) remain finite over the full interval $[0, T]$
A numerically computed extremal control (top row, left) for the effectiveness function $\varphi_1$ and $\varphi_2$ given in (23), (top row, right), $\alpha = (1, \ldots, 1)$, $\beta = 50 \alpha$, $\gamma = (600, 1400)$ and time horizon $T = 10$. The corresponding evolution of the individual cell populations is shown in the middle row (left) along with the evolution of the average of the populations (right). A comparison of the initial density $n(0, x)\equiv \frac{200}{21}$ and the terminal density $n(10, x)$ shown as red curve is given in the bottom row on the left. This control is not a strong local minimum as the entries of the solution $S$ to the Riccati differential equation (13) along the extremal controlled trajectory (shown in the bottom row, right) blow up in the interval
Examples of the system response, for the same parameters as in Fig. 3, for the constant controls $u_1 = 1$ and $u_2 = 2.5$
Illustration of the functions $\Upsilon$ (top, left), $\eta$ (top, right) and the full Hamiltonian $H$ (bottom) as function of the controls. These graphs depict the actual values for the solution shown in Fig. 6 below at the time $t = 19$
A numerically computed locally optimal control (top row, left) for $\alpha = (1, \ldots, 1)$, $\beta = 15 \alpha$, $\gamma = (90, 30)$ with therapy horizon $T = 20$ and effectiveness functions $\varphi_1 = \frac{2}{1+x^2}$ and $\varphi_2 = \frac{2}{1+3x^6}$ (top row, right). The corresponding evolution of the individual cell populations is shown in the middle row (left) along with the evolution of the average of the populations (right). A comparison of the initial density $n(0, x)\equiv \frac{200}{21}$ and the terminal density $n(20, x)$ shown as red curve is given in the bottom row on the left. This control is a strong local minimum as the entries of the solution $S$ to the Riccati differential equation (30) along the extremal controlled trajectory (shown in the bottom row, right) remain finite in the full interval $[0, 20]$
A numerically computed locally optimal control (top row, left) for α = (1, ...., 1), β = 15α, γ = (90, 30) with therapy horizon T = 20 and effectiveness functions ${\varphi _1}{\rm{ = }}\frac{2}{{1 + 5{x^2}}}$ and ${\varphi _2}{\rm{ = }}\frac{2}{{1 + 60{{\left( {1 - x} \right)}^6}}}$ (top row, right). The corresponding evolution of the individual cell populations is shown in the middle row (left) along with the evolution of the average of the populations (right). A comparison of the initial density n(0, x) ≡ $\frac{{200}}{{21}}$ and the terminal density n(20, x) shown as red curve is given in the bottom row on the left. This control is a strong local minimum as the entries of the solution S to the Riccati differential equation (30) along the extremal controlled trajectory (shown in the bottom row, right) blow up in the interval.
A numerically computed extremal control (top row, left) for α = (1, ..., 1), β = 15α, γ = (90, 30) with shortened therapy horizon T = 5 and effectiveness function ${\varphi _1}{\rm{ = }}\frac{2}{{1 + 5{x^2}}}$ and ${\varphi _2}{\rm{ = }}\frac{2}{{1 + 60{{\left( {1 - x} \right)}^6}}}$ (top row, right). The corresponding evolution of the individual cell populations is shown in the middle row (left) along with the evolution of the average of the populations (right). A comparison of the initial density n(0, x) ≡ $\frac{{200}}{{21}}$ and the terminal density n(5, x) shown as red curve is given in the bottom row on the left. This control still is a strong local minimum as the entries of the solution S to the Riccati differential equation (30) along the extremal controlled trajectory (shown in the bottom row, right) remain finite in the full interval [0, T].
Shooting Method
 Require: $N(t)>0$, $\forall t \in [0, T]$ and pick a threshold $\epsilon >0$. Ensure: optimal control   start with an initial guess $\lambda(0)$ and the stepsize $s=0.01$. while $\|\lambda(T)-\alpha\|>\epsilon$ do  solve $N, \; \lambda$.   $\lambda(0)\Leftarrow \lambda(0)-s\Big( \frac{\delta \lambda(T)}{\delta \lambda(0)}\Big) ^{-1} (\lambda(T)-\alpha)$ end while
 Require: $N(t)>0$, $\forall t \in [0, T]$ and pick a threshold $\epsilon >0$. Ensure: optimal control   start with an initial guess $\lambda(0)$ and the stepsize $s=0.01$. while $\|\lambda(T)-\alpha\|>\epsilon$ do  solve $N, \; \lambda$.   $\lambda(0)\Leftarrow \lambda(0)-s\Big( \frac{\delta \lambda(T)}{\delta \lambda(0)}\Big) ^{-1} (\lambda(T)-\alpha)$ end while
Iterative Method
 Require: $N(t)>0$, $\forall t \in [0, T]$ and pick a threshold $\epsilon >0$. Ensure: initial guess   Start with an initial trajectory $N^{(0)}$ and $k=1$. while $\|N^{(k)}-N^{(k-1)}\|>\epsilon$ do   $A^{(k-1)} \Leftarrow$ function of $N^{(k-1)}$   $O^{(k-1)} \Leftarrow$ function of $N^{(k-1)}$   Solve $S^{(k)}, \; \nu^{(k)}$ backward in time by assuming $\lambda^{(k)} = S^{(k)}N^{(k)}+\nu^{(k)}$.  Update $N^{(k)}$ forward in time:  $\qquad \dot{N}^{(k)} =[A^{(k-1)}+O^{(k-1)}S^{(k)}]N^{(k)}+O^{(k-1)}\nu^{(k)}$)  $u^{(k)} \Leftarrow$ function of $N^{(k)}$, $\lambda^{(k)}$  $k+1 \Leftarrow k$ end while
 Require: $N(t)>0$, $\forall t \in [0, T]$ and pick a threshold $\epsilon >0$. Ensure: initial guess   Start with an initial trajectory $N^{(0)}$ and $k=1$. while $\|N^{(k)}-N^{(k-1)}\|>\epsilon$ do   $A^{(k-1)} \Leftarrow$ function of $N^{(k-1)}$   $O^{(k-1)} \Leftarrow$ function of $N^{(k-1)}$   Solve $S^{(k)}, \; \nu^{(k)}$ backward in time by assuming $\lambda^{(k)} = S^{(k)}N^{(k)}+\nu^{(k)}$.  Update $N^{(k)}$ forward in time:  $\qquad \dot{N}^{(k)} =[A^{(k-1)}+O^{(k-1)}S^{(k)}]N^{(k)}+O^{(k-1)}\nu^{(k)}$)  $u^{(k)} \Leftarrow$ function of $N^{(k)}$, $\lambda^{(k)}$  $k+1 \Leftarrow k$ end while
 [1] Hong Niu, Zhijiang Feng, Qijin Xiao, Yajun Zhang. A PID control method based on optimal control strategy. Numerical Algebra, Control & Optimization, 2021, 11 (1) : 117-126. doi: 10.3934/naco.2020019 [2] Lars Grüne, Matthias A. Müller, Christopher M. Kellett, Steven R. Weller. Strict dissipativity for discrete time discounted optimal control problems. Mathematical Control & Related Fields, 2020  doi: 10.3934/mcrf.2020046 [3] Hui Lv, Xing'an Wang. Dissipative control for uncertain singular markovian jump systems via hybrid impulsive control. Numerical Algebra, Control & Optimization, 2021, 11 (1) : 127-142. doi: 10.3934/naco.2020020 [4] Awais Younus, Zoubia Dastgeer, Nudrat Ishaq, Abdul Ghaffar, Kottakkaran Sooppy Nisar, Devendra Kumar. On the observability of conformable linear time-invariant control systems. Discrete & Continuous Dynamical Systems - S, 2020  doi: 10.3934/dcdss.2020444 [5] Youming Guo, Tingting Li. Optimal control strategies for an online game addiction model with low and high risk exposure. Discrete & Continuous Dynamical Systems - B, 2020  doi: 10.3934/dcdsb.2020347 [6] Bernard Bonnard, Jérémy Rouot. Geometric optimal techniques to control the muscular force response to functional electrical stimulation using a non-isometric force-fatigue model. Journal of Geometric Mechanics, 2020  doi: 10.3934/jgm.2020032 [7] Xuefeng Zhang, Yingbo Zhang. Fault-tolerant control against actuator failures for uncertain singular fractional order systems. Numerical Algebra, Control & Optimization, 2021, 11 (1) : 1-12. doi: 10.3934/naco.2020011 [8] Zuliang Lu, Fei Huang, Xiankui Wu, Lin Li, Shang Liu. Convergence and quasi-optimality of $L^2-$norms based an adaptive finite element method for nonlinear optimal control problems. Electronic Research Archive, 2020, 28 (4) : 1459-1486. doi: 10.3934/era.2020077 [9] Laurence Cherfils, Stefania Gatti, Alain Miranville, Rémy Guillevin. Analysis of a model for tumor growth and lactate exchanges in a glioma. Discrete & Continuous Dynamical Systems - S, 2020  doi: 10.3934/dcdss.2020457 [10] Jianquan Li, Xin Xie, Dian Zhang, Jia Li, Xiaolin Lin. Qualitative analysis of a simple tumor-immune system with time delay of tumor action. Discrete & Continuous Dynamical Systems - B, 2020  doi: 10.3934/dcdsb.2020341 [11] Håkon Hoel, Gaukhar Shaimerdenova, Raúl Tempone. Multilevel Ensemble Kalman Filtering based on a sample average of independent EnKF estimators. Foundations of Data Science, 2020  doi: 10.3934/fods.2020017 [12] Zhiyan Ding, Qin Li, Jianfeng Lu. Ensemble Kalman Inversion for nonlinear problems: Weights, consistency, and variance bounds. Foundations of Data Science, 2020  doi: 10.3934/fods.2020018 [13] José Madrid, João P. G. Ramos. On optimal autocorrelation inequalities on the real line. Communications on Pure & Applied Analysis, 2021, 20 (1) : 369-388. doi: 10.3934/cpaa.2020271 [14] Tommi Brander, Joonas Ilmavirta, Petteri Piiroinen, Teemu Tyni. Optimal recovery of a radiating source with multiple frequencies along one line. Inverse Problems & Imaging, 2020, 14 (6) : 967-983. doi: 10.3934/ipi.2020044 [15] Lingju Kong, Roger Nichols. On principal eigenvalues of biharmonic systems. Communications on Pure & Applied Analysis, 2021, 20 (1) : 1-15. doi: 10.3934/cpaa.2020254 [16] Peizhao Yu, Guoshan Zhang, Yi Zhang. Decoupling of cubic polynomial matrix systems. Numerical Algebra, Control & Optimization, 2021, 11 (1) : 13-26. doi: 10.3934/naco.2020012 [17] Yongge Tian, Pengyang Xie. Simultaneous optimal predictions under two seemingly unrelated linear random-effects models. Journal of Industrial & Management Optimization, 2020  doi: 10.3934/jimo.2020168 [18] M. S. Lee, H. G. Harno, B. S. Goh, K. H. Lim. On the bang-bang control approach via a component-wise line search strategy for unconstrained optimization. Numerical Algebra, Control & Optimization, 2021, 11 (1) : 45-61. doi: 10.3934/naco.2020014 [19] Ilyasse Lamrani, Imad El Harraki, Ali Boutoulout, Fatima-Zahrae El Alaoui. Feedback stabilization of bilinear coupled hyperbolic systems. Discrete & Continuous Dynamical Systems - S, 2020  doi: 10.3934/dcdss.2020434 [20] Felix Finster, Jürg Fröhlich, Marco Oppio, Claudio F. Paganini. Causal fermion systems and the ETH approach to quantum theory. Discrete & Continuous Dynamical Systems - S, 2020  doi: 10.3934/dcdss.2020451

2019 Impact Factor: 1.27