# Optimal Feedback Control of Cancer Chemotherapy Using Hamilton–Jacobi–Bellman Equation

Authors: Yong Dam Jeong; Kwang Su Kim; Yunil Roh; Sooyoun Choi; Shingo Iwami; Il Hyo Jung
Journal: Complexity
 (2022)

---

## Abstract

Cancer chemotherapy has been the most common cancer treatment. However, it has side effects that kill both tumor cells and immune cells, which can ravage the patient’s immune system. Chemotherapy should be administered depending on the patient’s immunity as well as the level of cancer cells. Thus, we need to design an efficient treatment protocol. In this work, we study a feedback control problem of tumor-immune system to design an optimal chemotherapy strategy. For this, we first propose a mathematical model of tumor-immune interactions and conduct stability analysis of two equilibria. Next, the feedback control is found by solving the Hamilton–Jacobi–Bellman (HJB) equation. Here, we use an upwind finite-difference method for a numerical approximate solution of the HJB equation. Numerical simulations show that the feedback control can help determine the treatment protocol of chemotherapy for tumor and immune cells depending on the side effects.


---

## Body

## 1. Introduction

Cancer is known as one of leading causes of death worldwide. While many studies are under way to understand the tumor-immune interactions, it is the challenging problem to find an efficient cancer treatment in modern medicine. Chemotherapy has been one of the most common cancer treatment options, but because it kills both tumor cells and immune cells, it has side effects that destroy the patient’s immune system. Immunotherapy compensates for the deficiencies of chemotherapy by raising the effectiveness of immunity against tumor cells, so it is one of the most prominent cancer treatments, especially in a multifaceted approach [1]. However, immunotherapy may make immune cells over-activated, potentially causing autoimmune diseases. Furthermore, it is limited to cancer types and has expensive costs. Thus, chemotherapy is still being used as a common cancer treatment.There exist many research papers on mathematical modeling about tumor-immune system interactions [2–13]. De Boer et al. [14] presented a mathematical model based on ordinary differential equations (ODEs), which is for the macrophage T lymphocyte interactions generating an immune response against tumors. Furthermore, there are many studies on analyzing various cancer treatments through mathematical modeling [15–23]. De Pillis et al. [18] explored a computationally optimal strategy combining chemotherapy and immunotherapy treatments to minimize tumors and toxicity to patients while maximizing the immunity. Rihan et al. [15] studied an optimal control problem using a delay differential model (DDE) describing the tumor-immune interactions in presence of immuno-chemotherapy.Most of the papers mentioned above are related to the results of open-loop control which is a function of time. By solving the optimality system from Pontryagin’s maximum principle (PMP) [24], the control can be obtained. However, the principle is only a necessary condition for the control to be optimal and the change of states is not reflected in the open-loop control, so there are the difficulties in applying these results to biological phenomena practically.On the contrary, in the case of feedback control, since the perturbation of the state is considered, the shortcomings of open-loop control can be compensated. Feedback control can be found by solving a nonlinear partial differential equation called Hamilton–Jacobi–Bellman (HJB) equation which arises from the idea of dynamic programming [25]. Because the value function, which is the solution of HJB equation, was not generally guaranteed to be smooth, the concept of viscosity solution came up as an alternative [26], but this also could not be obtained analytically. Thus, numerical approximation method is needed to solve the problem.In this paper, we study an optimal feedback control problem of a tumor-immune model. To minimize tumor cells and cost associated with chemotherapy while reducing immune cells destroyed by chemotherapy, we find an optimal chemotherapy strategy that is dependent on tumor and immune cells. Furthermore, we investigate the change of optimal chemotherapy strategy according to the level of side effects and cost. For this, we derive the HJB equation through Bellman’s principle. As a numerical approximation method for the viscosity solution, an approximated solution for the HJB equation, an upwind finite-difference method [27], is used.The organization of this study is as follows. In Section2, we propose a mathematical model of the tumor-immune system and explain the parameters used in the model. In addition, we study the positiveness and boundedness of solutions of the model. In Section 3, we conduct the stability analysis for two equilibria. In Section 4, the optimal control problems are discussed. To design the feedback control problem, we define a value function and then derive the HJB equation. For solving this problem numerically, the upwind finite-difference method is described. Also, an open-loop control problem is considered for comparison with the feedback control. Numerical simulation results of the optimal chemotherapy strategy are presented in Section 5. Finally, we provide discussion and conclusions in Section 6.


## 2. Mathematical Model

### 2.1. Tumor-Immune Dynamics

We consider a mathematical model for tumor-immune dynamics, as represented by the scheme in Figure1. The model is given by the following differential equations:(1)dEdt=σ+cT−δE,dTdt=αT1−βT−ωET.Figure 1
Schematic diagram of model (1) for interaction of tumor cells and immune cells (commonly called effector cells).Tumor cellsT are assumed to grow in a logistic fashion. α is a growth rate of tumor cells and β is the inverse of carrying capacity for tumor cells. Effector cells E kill tumor cells at a net rate ω [7]. We assume that the effector cells are produced at a constant rate σ [28]. δ denotes the rate of effector cell turnover. Furthermore, it is assumed that, due to the direct presence of tumor cells, the effector cells grow at a rate c. The explanation of the parameters used in model (1) is in Table 1.Table 1
Variables and parameters in model (1).
SymbolDescriptionUnitσThe rate of effector cell proliferationcells day−1cThe antigenicity of tumor cellday−1δThe rate of effector cell turnoverday−1αThe growth rate of tumor cellday−1βThe inverse of tumor cell’s carrying capacitycells−1ωThe death rate of tumor cell by the interaction with an effector cellcells−1 day−1


### 2.2. Positiveness and Boundedness

To validate the mathematical model biologically, we should prove that solutions to system (1) are positive and bounded. For instance, it is not biologically feasible that a cell is negative. Furthermore, because the human body is made up of a finite number of cells, the cells should remain finite. These are summarized in the following theorems.Theorem 1.
LetEt,Tt be a solution of model (1). If the initial conditions satisfy E0>0 and T0>0, then the solution will remain positive, for all t>0.Proof.
Assume that there existst∗=inft>0|Et⋅Tt=0. First, we consider the case of Et∗=0. Then, we know that Et≥0 and Tt>0, for all t∈0,t∗. It follows that, for all t∈0,t∗,(2)dEdt=σ+cT−δE>−δE.
By Gronwall’s inequality, we obtain the inequality as follows. For allt∈0,t∗,(3)Et>E0eδt.
This is a contradiction with the fact thatEt∗=0.
Now, suppose thatTt∗=0 and note that Et>0, for all t∈0,t∗. It is followed from the extreme value theorem that E has positive maximum value EM. Then, it is obtained that, for all t∈0,t∗,(4)dTdt=αT1−βT−ωET≥−αβT2−ωEMT=−αβTT+ωEMαβ.
By the partial fraction decomposition, we obtain the inequality as follows; for allt∈0,t∗,(5)Tt≥BωEM/αβe−ωEMt1−Be−ωEMt>0,where B=T0/T0+ωEM/αβ. This contradicts the assumption of Tt∗=0. Thus, Et,Tt will remain positive.Theorem 2.
LetEt,Tt be a solution of model (1). If E0≥0 and T0≥0, then Et and Tt are bounded.Proof.
SupposeTt is unbounded. First, we consider the case of T0>1/β. Let t1=inft>0|Tt=1/β. Then, it follows that dT/dt|t=t1<0. However, since Tt is unbounded, there exists t2>t1 such that Tt2=1/β, which also implies that dT/dt|t=t2<0. By the recursive procedures, we obtain that Tt≤1/β for t≥t1. This contradicts the assumption. So, Tt≤T0, for all t>0. Now, consider T0≤1/β. Let t∗=inft>0|Tt>1/β. Then, there exists t¯∈0,t∗ such that Tt¯=1/β. It follows that dT/dt|t=t¯<0. The remaining proof is similar to the first case. Thus, Tt≤1/β, for all t>0. Therefore, Tt is bounded. Then, we can obtain the following inequality:(6)dEdt=σ+cT−δE<σ+cTM−δE,where TM=maxT0,1/β. It is followed that(7)Et<σ+cTMδ+E0−σ+cTMδe−δt.
Taking thelimsup of both sides, we obtain(8)limsupt⟶∞Et≤σ+cTMδ.
Therefore,Et is also bounded.


## 2.1. Tumor-Immune Dynamics

We consider a mathematical model for tumor-immune dynamics, as represented by the scheme in Figure1. The model is given by the following differential equations:(1)dEdt=σ+cT−δE,dTdt=αT1−βT−ωET.Figure 1
Schematic diagram of model (1) for interaction of tumor cells and immune cells (commonly called effector cells).Tumor cellsT are assumed to grow in a logistic fashion. α is a growth rate of tumor cells and β is the inverse of carrying capacity for tumor cells. Effector cells E kill tumor cells at a net rate ω [7]. We assume that the effector cells are produced at a constant rate σ [28]. δ denotes the rate of effector cell turnover. Furthermore, it is assumed that, due to the direct presence of tumor cells, the effector cells grow at a rate c. The explanation of the parameters used in model (1) is in Table 1.Table 1
Variables and parameters in model (1).
SymbolDescriptionUnitσThe rate of effector cell proliferationcells day−1cThe antigenicity of tumor cellday−1δThe rate of effector cell turnoverday−1αThe growth rate of tumor cellday−1βThe inverse of tumor cell’s carrying capacitycells−1ωThe death rate of tumor cell by the interaction with an effector cellcells−1 day−1


## 2.2. Positiveness and Boundedness

To validate the mathematical model biologically, we should prove that solutions to system (1) are positive and bounded. For instance, it is not biologically feasible that a cell is negative. Furthermore, because the human body is made up of a finite number of cells, the cells should remain finite. These are summarized in the following theorems.Theorem 1.
LetEt,Tt be a solution of model (1). If the initial conditions satisfy E0>0 and T0>0, then the solution will remain positive, for all t>0.Proof.
Assume that there existst∗=inft>0|Et⋅Tt=0. First, we consider the case of Et∗=0. Then, we know that Et≥0 and Tt>0, for all t∈0,t∗. It follows that, for all t∈0,t∗,(2)dEdt=σ+cT−δE>−δE.
By Gronwall’s inequality, we obtain the inequality as follows. For allt∈0,t∗,(3)Et>E0eδt.
This is a contradiction with the fact thatEt∗=0.
Now, suppose thatTt∗=0 and note that Et>0, for all t∈0,t∗. It is followed from the extreme value theorem that E has positive maximum value EM. Then, it is obtained that, for all t∈0,t∗,(4)dTdt=αT1−βT−ωET≥−αβT2−ωEMT=−αβTT+ωEMαβ.
By the partial fraction decomposition, we obtain the inequality as follows; for allt∈0,t∗,(5)Tt≥BωEM/αβe−ωEMt1−Be−ωEMt>0,where B=T0/T0+ωEM/αβ. This contradicts the assumption of Tt∗=0. Thus, Et,Tt will remain positive.Theorem 2.
LetEt,Tt be a solution of model (1). If E0≥0 and T0≥0, then Et and Tt are bounded.Proof.
SupposeTt is unbounded. First, we consider the case of T0>1/β. Let t1=inft>0|Tt=1/β. Then, it follows that dT/dt|t=t1<0. However, since Tt is unbounded, there exists t2>t1 such that Tt2=1/β, which also implies that dT/dt|t=t2<0. By the recursive procedures, we obtain that Tt≤1/β for t≥t1. This contradicts the assumption. So, Tt≤T0, for all t>0. Now, consider T0≤1/β. Let t∗=inft>0|Tt>1/β. Then, there exists t¯∈0,t∗ such that Tt¯=1/β. It follows that dT/dt|t=t¯<0. The remaining proof is similar to the first case. Thus, Tt≤1/β, for all t>0. Therefore, Tt is bounded. Then, we can obtain the following inequality:(6)dEdt=σ+cT−δE<σ+cTM−δE,where TM=maxT0,1/β. It is followed that(7)Et<σ+cTMδ+E0−σ+cTMδe−δt.
Taking thelimsup of both sides, we obtain(8)limsupt⟶∞Et≤σ+cTMδ.
Therefore,Et is also bounded.


## 3. Stability Analysis

### 3.1. Local Stability

From model (1), two equilibria can be obtained. One of these is the tumor-free equilibrium, that is,(9)P0=E0,T0=σδ,0.To check the local stability, we compute the Jacobian atP0 as follows:(10)JP0=−δc0α−σωδ.If the trace is negative and the determinant is positive, thenP0 is locally asymptotically stable. Thus, the condition αδ−σω>0 should be satisfied. We define the criterion K as follows:(11)K=αδσω.The local stability ofP0 is summarized in the following theorem.Theorem 3.
IfK<1, model (1) has an tumor-free equilibrium P0 which is locally asymptotically stable.
The other equilibrium isP1=E1,T1, where(12)E1=σ+cT1δ,T1=αδ−σωαβδ+cω.IfK>1, the positiveness of P1 is guaranteed. Note that if K<1, the equilibrium P1 does not exist and P0=P1 when K=1. The local stability of the equilibrium P1 is summarized in the next theorem.Theorem 4.
IfK>1, model (1) has an unstable equilibrium P0 and a locally asymptotically stable positive equilibrium P1.Proof.
Computing the Jacobian atP1 gives(13)JP1=−δcωT1α−2αβT1−ωE1.
To show thatP1 is locally asymptotically stable, we compute the following conditions:(14)TrJP1=−δ+α−2αβT1−ωE1=−δ2+αδ−2αβδT1−cωT1−σωδ=−δ−αβK−1αβδ+cω,DetJP1=−αδ+2αβδT1+ωδE1+cωT1=αδ−σω+2αβδT1+2cωT1=αδ−σω.=σωK−1.
SinceK>1, the trace is negative and the determinant is positive. Therefore, the equilibrium P1 is locally asymptotically stable. In this case, the tumor-free equilibrium P0 is unstable, which is a saddle.


### 3.2. Global Stability

Furthermore, we investigate the global stability of the equilibriaP0 and P1. These are summarized in the following theorems.Theorem 5.
IfK≤1, the tumor-free equilibrium P0 of model (1) is globally asymptotically stable.Proof.
Consider the following Lyapunov function onℝ+2:(15)L=E−E0−E0lnEE0+cαT.
It is easy to see thatL=0 at the equilibrium P0 and L>0, for all E,T≠σ/δ,0. Thus, L is positively definite and is also clearly radially unbounded. Taking the derivative of L with respect to t, we obtain that(16)dLdt=1−E0EE′+cαT′=1−E0Eσ+cT−δE+cT1−βT−cωαET=−cβT2+σ+δE0−δE−σE0E+2cT−cE0TE−cωαET=−cβT2+σ2−δEσ−σδE+cT2−σδE−δEσK.
By the arithmetic-geometric mean and the assumptionK≤1, we obtain the following inequality:(17)dLdt≤−cβT2<0.
Hence, we haveL′<0, for all E,T≠σ/δ,0. Therefore, by Lyapunov’s theorem, the tumor-free equilibrium P0 is globally asymptotically stable.Theorem 6.
AssumeK>1; the equilibrium P1 of model (1) is globally asymptotically stable whenever T0>0.Proof.
LetX=ℝ+2 and set the functions:(18)F1E,T=σ+cT−δE,F2E,T=αT1−βT−ωET.
Define the functionDE,T=1/T on X. Then, D is continuously differentiable on X. Furthermore, we obtain that(19)∂DF1∂E+∂DF2∂T=−δT−αβ<0.
Applying Dulac–Bendixson criterion, system (1) has no periodic orbits or graphics on X. Since Theorem 2 guaranteed the boundedness of solutions of system (1), X is positively invariant and contains an omega limit set of every initial conditions. Thus, Poincaré–Bendixson theorem can be applied. When K>1, if T0=0, then the solutions will remain on the E-axis in E,T-plane and converge to P0. If T0>0, that is, E0,T0∈X, then we want to show that the equilibrium P0 does not belong to the omega limit set of the point E0,T0. Suppose the omega limit set contains P0. By Theorem 4, P0 is unstable saddle, and so it has stable manifold which is given by the E-axis. Because each solution starting from T0=0 remains on the E-axis and converges to P0, thus, the stable manifold of P0 is not in X. Hence, the omega limit set should contain another equilibrium P1. But, since P1 is locally asymptotically stable, the solutions starting close to P1, approach to it. So, the equilibrium P0 is not in the omega limit set. Therefore, the omega limit set of E0,T0 contains only P1. □


## 3.1. Local Stability

From model (1), two equilibria can be obtained. One of these is the tumor-free equilibrium, that is,(9)P0=E0,T0=σδ,0.To check the local stability, we compute the Jacobian atP0 as follows:(10)JP0=−δc0α−σωδ.If the trace is negative and the determinant is positive, thenP0 is locally asymptotically stable. Thus, the condition αδ−σω>0 should be satisfied. We define the criterion K as follows:(11)K=αδσω.The local stability ofP0 is summarized in the following theorem.Theorem 3.
IfK<1, model (1) has an tumor-free equilibrium P0 which is locally asymptotically stable.
The other equilibrium isP1=E1,T1, where(12)E1=σ+cT1δ,T1=αδ−σωαβδ+cω.IfK>1, the positiveness of P1 is guaranteed. Note that if K<1, the equilibrium P1 does not exist and P0=P1 when K=1. The local stability of the equilibrium P1 is summarized in the next theorem.Theorem 4.
IfK>1, model (1) has an unstable equilibrium P0 and a locally asymptotically stable positive equilibrium P1.Proof.
Computing the Jacobian atP1 gives(13)JP1=−δcωT1α−2αβT1−ωE1.
To show thatP1 is locally asymptotically stable, we compute the following conditions:(14)TrJP1=−δ+α−2αβT1−ωE1=−δ2+αδ−2αβδT1−cωT1−σωδ=−δ−αβK−1αβδ+cω,DetJP1=−αδ+2αβδT1+ωδE1+cωT1=αδ−σω+2αβδT1+2cωT1=αδ−σω.=σωK−1.
SinceK>1, the trace is negative and the determinant is positive. Therefore, the equilibrium P1 is locally asymptotically stable. In this case, the tumor-free equilibrium P0 is unstable, which is a saddle.


## 3.2. Global Stability

Furthermore, we investigate the global stability of the equilibriaP0 and P1. These are summarized in the following theorems.Theorem 5.
IfK≤1, the tumor-free equilibrium P0 of model (1) is globally asymptotically stable.Proof.
Consider the following Lyapunov function onℝ+2:(15)L=E−E0−E0lnEE0+cαT.
It is easy to see thatL=0 at the equilibrium P0 and L>0, for all E,T≠σ/δ,0. Thus, L is positively definite and is also clearly radially unbounded. Taking the derivative of L with respect to t, we obtain that(16)dLdt=1−E0EE′+cαT′=1−E0Eσ+cT−δE+cT1−βT−cωαET=−cβT2+σ+δE0−δE−σE0E+2cT−cE0TE−cωαET=−cβT2+σ2−δEσ−σδE+cT2−σδE−δEσK.
By the arithmetic-geometric mean and the assumptionK≤1, we obtain the following inequality:(17)dLdt≤−cβT2<0.
Hence, we haveL′<0, for all E,T≠σ/δ,0. Therefore, by Lyapunov’s theorem, the tumor-free equilibrium P0 is globally asymptotically stable.Theorem 6.
AssumeK>1; the equilibrium P1 of model (1) is globally asymptotically stable whenever T0>0.Proof.
LetX=ℝ+2 and set the functions:(18)F1E,T=σ+cT−δE,F2E,T=αT1−βT−ωET.
Define the functionDE,T=1/T on X. Then, D is continuously differentiable on X. Furthermore, we obtain that(19)∂DF1∂E+∂DF2∂T=−δT−αβ<0.
Applying Dulac–Bendixson criterion, system (1) has no periodic orbits or graphics on X. Since Theorem 2 guaranteed the boundedness of solutions of system (1), X is positively invariant and contains an omega limit set of every initial conditions. Thus, Poincaré–Bendixson theorem can be applied. When K>1, if T0=0, then the solutions will remain on the E-axis in E,T-plane and converge to P0. If T0>0, that is, E0,T0∈X, then we want to show that the equilibrium P0 does not belong to the omega limit set of the point E0,T0. Suppose the omega limit set contains P0. By Theorem 4, P0 is unstable saddle, and so it has stable manifold which is given by the E-axis. Because each solution starting from T0=0 remains on the E-axis and converges to P0, thus, the stable manifold of P0 is not in X. Hence, the omega limit set should contain another equilibrium P1. But, since P1 is locally asymptotically stable, the solutions starting close to P1, approach to it. So, the equilibrium P0 is not in the omega limit set. Therefore, the omega limit set of E0,T0 contains only P1. □


## 4. Optimal Control Problems

In this section, we consider optimal control problems of model (1) to investigate the optimal chemotherapy strategy depending on the level of effector cells and tumor cells. To begin with, we design a feedback control problem. Additionally, an open-loop control problem of chemotherapy is considered for comparison with the feedback control.First, we simplify model (1) into the nondimensional form through the following process:(20)x=x1,x2=EE∗,TT∗,where E∗=σ+c/β/δ and T∗=1/β. Then, the new dependent variables x1t and x2t are dimensionless quantities. Hence, the system for them becomes(21)dx1dt=f1t,x=σ¯+c¯x2−δx1,dx2dt=f2t,x=αx21−x2−ω¯x1x2,where σ¯=σ/E∗, c¯=cT∗/E∗, and ω¯=ωE∗.


### 4.1. Feedback Control Problem

#### 4.1.1. Hamilton–Jacobi–Bellman Equation

Letτ and tf be an initial time and a final time, respectively, and let u be a measurable function from τ,tf to 0,1. The control ut means the treatment effort in time t, which represents the intensity of chemotherapy. The class of all admissible controls is denoted by(22)Uf=L∞τ,tf;0,1.Then, the system with the control is as follows:(23)dx1dt=f1t,x,u=σ¯+c¯x2−δx1−qEux1,dx2dt=f2t,x,u=αx21−x2−ω¯x1x2−qTux2,where qE and qT are the chemotherapy-induced death coefficients of effector cell and tumor cell, respectively. Especially, the coefficient qE is related to the side effect of chemotherapy. Thus, to reduce tumor cells and death of effector cells via chemotherapy while minimizing the cost to chemotherapy, we consider the following objective functional:(24)Jτ,a,b,u=∫τtfLt,xt,utdt,where the Lagrangian Lt,xt,ut=qEutx1t+1−qTutx2t+B/2ut2 and xτ=x1τ,x2τ=a,b is an initial condition. In the Lagrangian L, first and second terms mean the death of effector cells due to side effects and the tumor cells remaining after chemotherapy, respectively. The last term represents the cost of chemotherapy.Define the value function:(25)vτ,xτ=minu∈UfJτ,xτ,u.The optimal control problem can be formulated as the HJB equation via Bellman’s principle of the dynamic programming. Thus, we can obtain the HJB equation as follows:(26)vt+minu∈Uf∇xv⋅ft,x,u+Lt,x,u=0,∀t,x∈0,tf×ℝ2,where f=f1,f2 and the initial condition vtf,x=0. The value function in characterizing the optimal feedback law is described in the following two propositions [29].Proposition 1.
Letv∈C10,tf×ℝn be the value function. Then, if there exists a control u∗⋅∈Uf such that(27)∇xvt,x∗t⋅ft,x∗t,u∗t+Lt,x∗t,u∗t=minu∈Uf∇xvt,x∗t⋅ft,x∗t,u+Lt,x∗t,u.then u∗⋅ is an optimal control, where x∗=x1∗,x2∗ is the state corresponding to u∗.The optimal control is denoted as(28)u∗t∈argminu∈Uf∇xvt,x∗t⋅ft,x∗t,u+Lt,x∗t,u,for almost all t∈0,T.Proposition 2.
Letv∈C10,tf×ℝ2 be the value function. Then, u∗t,x∗t is an optimal control-trajectory pair in feedback form if and only if(29)vtt,x∗t+∇xvt,x∗t⋅ft,x∗t,u∗t+Lt,x∗t,u∗t=0,for almost all t∈0,T.
By Propositions1 and 2, the feedback law via the value function can be constructed. This is summarized in the following theorem.Theorem 7.
Letv∈C10,tf×ℝ2 be the value function. Suppose ut,x satisfies(30)∇xvt,x⋅ft,x,ut,x+Lt,x,ut,x=minu∈Uf∇xvt,x⋅ft,x,ut,x+Lt,x,ut,x.
Then,uz∗t=ut,xzt is the feedback law of the optimal control problem, where xzt is the solution of the following control system:(31)xz′t=ft,xzt,ut,xzt,∀xz0=z,t∈0,tf.


#### 4.1.2. Upwind Finite-Difference Method

Solutions of HJB equation are continuous, but nonsmooth inℝn. In order to complement the nonsmooth, the idea of viscosity solutions is used by Lions [26]. The idea is to consider the convection-diffusion equation by adding the diffusion term ϵ∇2v, where ϵ>0 is the viscosity. Then, the solution of this equation becomes smooth and converges to the viscosity solution of original HJB equation as ϵ⟶0 [30]. It can be shown that the existence and uniqueness of viscosity solution to HJB equation are guaranteed [30, 31]. However, this equation generally is not analytically solvable and so we use the upwind finite-difference method as a numerical approximation method.Before addressing this method, we settle for computational problems that arise from using the finite-difference method. Since the variables of model (1) have wide ranges of numbers, solving a feedback control problem requires a high amount of numerical computations. So, by scaling the variables to the appropriate level, we transformed model (1) to dimensionless model (21). Assuming that E0≤E∗ and T0≤T∗, it is followed from Theorem 2 that xt=x1t,x2t∈0,1×0,1, for all t>0. Therefore, we restrict the feedback control problem to the domain 0,tf×Ω, where Ω=0,1×0,1. Then, we can reformulate equation (26) as follows:(32)ε∇2v+vt+minu∈Uf∇xvt,x⋅ft,x,u+Lt,x,u=0,∀t,x∈0,tf×Ω,vt,x=gt,x,∀t,x∈0,tf×∂Ω,vtf,x=0,∀x∈Ω.We need to be determine a suitable boundary condition, fort,x∈0,tf×∂Ω. Since the only possible information is the value function at the final time tf, we pose the Dirichlet boundary condition, that is, gt,x=vtf,x [32].Here, we apply the upwind finite-difference method to (32) with ϵ≪1 [33, 34]. Let i=i1,i2,i1+=i1+1,i2,i1−=i1−1,i2,i2+=i1,i2+1, and i2−=i1,i2−1. To discretize (32), we assume that, for k=1,2, xk,ik=ikΔxk, ik=0,1,…,Mk, with Δxk=1/Mk, and tj=j−1Δt, j=1,2,…,Mt+1, with Δt=tf/Mt, are the partitions along the xk and t axes, respectively. Here, we assume the meshes for state variables are uniform, that is, M1=M2=M. Then, two stencils can be considered. One consists of the six vertices x1,i1−1,tj, x1,i1,tj, x1,i1+1,tj, x1,i1−1,tj+1, x1,i1,tj+1, and x1,i1+1,tj+1. The other consists of the six vertices x2,i2−1,tj, x2,i2,tj, x2,i2+1,tj, x2,i2−1,tj+1, x2,i2,tj+1, and x2,i2+1,tj+1. On the cuboid consisting these stencils, (32) can be approximated by the difference equations as follows; for i1,i2=j,…,M−j and j=1,2,…,Mt,(33)vij+1−vijΔt=∑k=121+sgnfk,ijuij2fk,ijuijvik+j−vijΔxk+1−sgnfk,ijuij2fk,ijuijvij−vik−jΔxk+Lijuij,uij=argsupu∈Uf−∑k=121+sgnfk,iju2fk,ijuvik+j−vijΔxk+1−sgnfk,iju2fk,ijuvij−vik−jΔxk−Liju,where(34)f1,ijuij=f1tj,x1,i1,x2,i2,uij,f2,ijuij=f2tj,x1,i1,x2,i2,uij,Lijuij=Ltj,x1,i1,x2,i2,uij.Based on the upwind finite-difference method, the algorithm for solving HJB equation is summarized in Algorithm1.Algorithm 1: Algorithm for solving HJB equation
Set the initial conditionvi1ik=0Mk and the feedback control ui1ik=0Mk for k=1,2forj=1,2,…,Mtdofori1=j,j+1,…,M1−jdofori2=j,j+1,…,M2−jdoWiju=−∑k=121+sgnfk,iju/2fk,ijuvik+j−vij/Δxk+1−sgnfk,iju/2fk,ijuvij−vik−j/Δxk−Lijuuij=argsupu∈UfWijuvij+1=vij−ΔtWijuijendendend


### 4.2. Open-Loop Control Problem

In open-loop control problem,tf and u are the same as meanings in the feedback control problem. Then, an objective functional for open-loop control problem is given by(35)Gu=∫0tfLt,xt,utdt,where x0=x10,x20 is an appropriate initial condition. Similarly, the open-loop control problem is to find an optimal control u∗ such that(36)Gu∗=minGu,u∈Uo.subject to system (23), where the control set is defined as(37)Uo=u|utis Lebesgue measurable,0≤ut≤1∀t∈0,tf.First, an existence of optimal control for problem (36) and the optimality system is derived.Theorem 8.
Given the objective functionalGu and the control set Uo, there exist an optimal control u∗ such that Gu∗=minGu,u∈Uo.Proof.
From [35], we prove the existence of an optimal control. By Theorem 1, the state variables are nonnegative. Note that the control is also nonnegative. The set of all controls uinUo is convex and closed. For the minimizing problem, the convexity of G in u is satisfied and the optimality system is bounded. So, the compactness necessary for the existence of optimal control is guaranteed. Furthermore, the Lagrangian L is convex on the control set Uo. Since the state variables are bounded, we can easily confirm that there exist a constant κ>1 and numbers ρ1,ρ2 such that(38)Gu≥ρ1u2κ/2−ρ2=ρ1uκ−ρ2.
If follows the existence of an optimal control.
In order to find a minimal value of the Lagrangian, we define a Hamiltonian:(39)H=x2+qEux1−qTux2+B2u2+λ1σ¯+c¯x2−δx1−qEux1+λ2αx21−x2−ω¯x1x2−qTux2,where λ1 and λ2 are the adjoint variables. To derive a necessary conditions for the optimal solution of optimal control problem (23) and (35), we use Pontryagin’s Maximum Principle [24]. The adjoint system and control characterization are described in the following theorem.Theorem 9.
Given an optimal controlu∗ and a solution x∗=x1∗,x2∗, there exist adjoint variables λ1 and λ2 satisfying(40)λ1′t=−∂H∂x1=qEut+λ1tδ+qEut+λ2tω¯x2t,λ2′t=−∂H∂x2=−1+qTut−c¯λ1t−λ2tα−2αx2t−ω¯x1t−qTut.under transversality conditions,(41)λitf=0,i=1,2.
Moreover, the control functionu∗ is given by(42)u∗=min1,max0,R,where R=λ1−1qEx1+λ2+1qTx2/B.Proof.
From Pontryagin’s Maximum Principle [24], the adjoint system can be obtained:(43)λ1′t=−∂H∂x1,λ2′t=−∂H∂x2.with λitf=0 for i=1,2. By solving the equation(44)∂H∂u=0and using properties of the control set Uo, we can characterize the optimal control (42).


## 4.1. Feedback Control Problem

### 4.1.1. Hamilton–Jacobi–Bellman Equation

Letτ and tf be an initial time and a final time, respectively, and let u be a measurable function from τ,tf to 0,1. The control ut means the treatment effort in time t, which represents the intensity of chemotherapy. The class of all admissible controls is denoted by(22)Uf=L∞τ,tf;0,1.Then, the system with the control is as follows:(23)dx1dt=f1t,x,u=σ¯+c¯x2−δx1−qEux1,dx2dt=f2t,x,u=αx21−x2−ω¯x1x2−qTux2,where qE and qT are the chemotherapy-induced death coefficients of effector cell and tumor cell, respectively. Especially, the coefficient qE is related to the side effect of chemotherapy. Thus, to reduce tumor cells and death of effector cells via chemotherapy while minimizing the cost to chemotherapy, we consider the following objective functional:(24)Jτ,a,b,u=∫τtfLt,xt,utdt,where the Lagrangian Lt,xt,ut=qEutx1t+1−qTutx2t+B/2ut2 and xτ=x1τ,x2τ=a,b is an initial condition. In the Lagrangian L, first and second terms mean the death of effector cells due to side effects and the tumor cells remaining after chemotherapy, respectively. The last term represents the cost of chemotherapy.Define the value function:(25)vτ,xτ=minu∈UfJτ,xτ,u.The optimal control problem can be formulated as the HJB equation via Bellman’s principle of the dynamic programming. Thus, we can obtain the HJB equation as follows:(26)vt+minu∈Uf∇xv⋅ft,x,u+Lt,x,u=0,∀t,x∈0,tf×ℝ2,where f=f1,f2 and the initial condition vtf,x=0. The value function in characterizing the optimal feedback law is described in the following two propositions [29].Proposition 1.
Letv∈C10,tf×ℝn be the value function. Then, if there exists a control u∗⋅∈Uf such that(27)∇xvt,x∗t⋅ft,x∗t,u∗t+Lt,x∗t,u∗t=minu∈Uf∇xvt,x∗t⋅ft,x∗t,u+Lt,x∗t,u.then u∗⋅ is an optimal control, where x∗=x1∗,x2∗ is the state corresponding to u∗.The optimal control is denoted as(28)u∗t∈argminu∈Uf∇xvt,x∗t⋅ft,x∗t,u+Lt,x∗t,u,for almost all t∈0,T.Proposition 2.
Letv∈C10,tf×ℝ2 be the value function. Then, u∗t,x∗t is an optimal control-trajectory pair in feedback form if and only if(29)vtt,x∗t+∇xvt,x∗t⋅ft,x∗t,u∗t+Lt,x∗t,u∗t=0,for almost all t∈0,T.
By Propositions1 and 2, the feedback law via the value function can be constructed. This is summarized in the following theorem.Theorem 7.
Letv∈C10,tf×ℝ2 be the value function. Suppose ut,x satisfies(30)∇xvt,x⋅ft,x,ut,x+Lt,x,ut,x=minu∈Uf∇xvt,x⋅ft,x,ut,x+Lt,x,ut,x.
Then,uz∗t=ut,xzt is the feedback law of the optimal control problem, where xzt is the solution of the following control system:(31)xz′t=ft,xzt,ut,xzt,∀xz0=z,t∈0,tf.


### 4.1.2. Upwind Finite-Difference Method

Solutions of HJB equation are continuous, but nonsmooth inℝn. In order to complement the nonsmooth, the idea of viscosity solutions is used by Lions [26]. The idea is to consider the convection-diffusion equation by adding the diffusion term ϵ∇2v, where ϵ>0 is the viscosity. Then, the solution of this equation becomes smooth and converges to the viscosity solution of original HJB equation as ϵ⟶0 [30]. It can be shown that the existence and uniqueness of viscosity solution to HJB equation are guaranteed [30, 31]. However, this equation generally is not analytically solvable and so we use the upwind finite-difference method as a numerical approximation method.Before addressing this method, we settle for computational problems that arise from using the finite-difference method. Since the variables of model (1) have wide ranges of numbers, solving a feedback control problem requires a high amount of numerical computations. So, by scaling the variables to the appropriate level, we transformed model (1) to dimensionless model (21). Assuming that E0≤E∗ and T0≤T∗, it is followed from Theorem 2 that xt=x1t,x2t∈0,1×0,1, for all t>0. Therefore, we restrict the feedback control problem to the domain 0,tf×Ω, where Ω=0,1×0,1. Then, we can reformulate equation (26) as follows:(32)ε∇2v+vt+minu∈Uf∇xvt,x⋅ft,x,u+Lt,x,u=0,∀t,x∈0,tf×Ω,vt,x=gt,x,∀t,x∈0,tf×∂Ω,vtf,x=0,∀x∈Ω.We need to be determine a suitable boundary condition, fort,x∈0,tf×∂Ω. Since the only possible information is the value function at the final time tf, we pose the Dirichlet boundary condition, that is, gt,x=vtf,x [32].Here, we apply the upwind finite-difference method to (32) with ϵ≪1 [33, 34]. Let i=i1,i2,i1+=i1+1,i2,i1−=i1−1,i2,i2+=i1,i2+1, and i2−=i1,i2−1. To discretize (32), we assume that, for k=1,2, xk,ik=ikΔxk, ik=0,1,…,Mk, with Δxk=1/Mk, and tj=j−1Δt, j=1,2,…,Mt+1, with Δt=tf/Mt, are the partitions along the xk and t axes, respectively. Here, we assume the meshes for state variables are uniform, that is, M1=M2=M. Then, two stencils can be considered. One consists of the six vertices x1,i1−1,tj, x1,i1,tj, x1,i1+1,tj, x1,i1−1,tj+1, x1,i1,tj+1, and x1,i1+1,tj+1. The other consists of the six vertices x2,i2−1,tj, x2,i2,tj, x2,i2+1,tj, x2,i2−1,tj+1, x2,i2,tj+1, and x2,i2+1,tj+1. On the cuboid consisting these stencils, (32) can be approximated by the difference equations as follows; for i1,i2=j,…,M−j and j=1,2,…,Mt,(33)vij+1−vijΔt=∑k=121+sgnfk,ijuij2fk,ijuijvik+j−vijΔxk+1−sgnfk,ijuij2fk,ijuijvij−vik−jΔxk+Lijuij,uij=argsupu∈Uf−∑k=121+sgnfk,iju2fk,ijuvik+j−vijΔxk+1−sgnfk,iju2fk,ijuvij−vik−jΔxk−Liju,where(34)f1,ijuij=f1tj,x1,i1,x2,i2,uij,f2,ijuij=f2tj,x1,i1,x2,i2,uij,Lijuij=Ltj,x1,i1,x2,i2,uij.Based on the upwind finite-difference method, the algorithm for solving HJB equation is summarized in Algorithm1.Algorithm 1: Algorithm for solving HJB equation
Set the initial conditionvi1ik=0Mk and the feedback control ui1ik=0Mk for k=1,2forj=1,2,…,Mtdofori1=j,j+1,…,M1−jdofori2=j,j+1,…,M2−jdoWiju=−∑k=121+sgnfk,iju/2fk,ijuvik+j−vij/Δxk+1−sgnfk,iju/2fk,ijuvij−vik−j/Δxk−Lijuuij=argsupu∈UfWijuvij+1=vij−ΔtWijuijendendend


## 4.1.1. Hamilton–Jacobi–Bellman Equation

Letτ and tf be an initial time and a final time, respectively, and let u be a measurable function from τ,tf to 0,1. The control ut means the treatment effort in time t, which represents the intensity of chemotherapy. The class of all admissible controls is denoted by(22)Uf=L∞τ,tf;0,1.Then, the system with the control is as follows:(23)dx1dt=f1t,x,u=σ¯+c¯x2−δx1−qEux1,dx2dt=f2t,x,u=αx21−x2−ω¯x1x2−qTux2,where qE and qT are the chemotherapy-induced death coefficients of effector cell and tumor cell, respectively. Especially, the coefficient qE is related to the side effect of chemotherapy. Thus, to reduce tumor cells and death of effector cells via chemotherapy while minimizing the cost to chemotherapy, we consider the following objective functional:(24)Jτ,a,b,u=∫τtfLt,xt,utdt,where the Lagrangian Lt,xt,ut=qEutx1t+1−qTutx2t+B/2ut2 and xτ=x1τ,x2τ=a,b is an initial condition. In the Lagrangian L, first and second terms mean the death of effector cells due to side effects and the tumor cells remaining after chemotherapy, respectively. The last term represents the cost of chemotherapy.Define the value function:(25)vτ,xτ=minu∈UfJτ,xτ,u.The optimal control problem can be formulated as the HJB equation via Bellman’s principle of the dynamic programming. Thus, we can obtain the HJB equation as follows:(26)vt+minu∈Uf∇xv⋅ft,x,u+Lt,x,u=0,∀t,x∈0,tf×ℝ2,where f=f1,f2 and the initial condition vtf,x=0. The value function in characterizing the optimal feedback law is described in the following two propositions [29].Proposition 1.
Letv∈C10,tf×ℝn be the value function. Then, if there exists a control u∗⋅∈Uf such that(27)∇xvt,x∗t⋅ft,x∗t,u∗t+Lt,x∗t,u∗t=minu∈Uf∇xvt,x∗t⋅ft,x∗t,u+Lt,x∗t,u.then u∗⋅ is an optimal control, where x∗=x1∗,x2∗ is the state corresponding to u∗.The optimal control is denoted as(28)u∗t∈argminu∈Uf∇xvt,x∗t⋅ft,x∗t,u+Lt,x∗t,u,for almost all t∈0,T.Proposition 2.
Letv∈C10,tf×ℝ2 be the value function. Then, u∗t,x∗t is an optimal control-trajectory pair in feedback form if and only if(29)vtt,x∗t+∇xvt,x∗t⋅ft,x∗t,u∗t+Lt,x∗t,u∗t=0,for almost all t∈0,T.
By Propositions1 and 2, the feedback law via the value function can be constructed. This is summarized in the following theorem.Theorem 7.
Letv∈C10,tf×ℝ2 be the value function. Suppose ut,x satisfies(30)∇xvt,x⋅ft,x,ut,x+Lt,x,ut,x=minu∈Uf∇xvt,x⋅ft,x,ut,x+Lt,x,ut,x.
Then,uz∗t=ut,xzt is the feedback law of the optimal control problem, where xzt is the solution of the following control system:(31)xz′t=ft,xzt,ut,xzt,∀xz0=z,t∈0,tf.


## 4.1.2. Upwind Finite-Difference Method

Solutions of HJB equation are continuous, but nonsmooth inℝn. In order to complement the nonsmooth, the idea of viscosity solutions is used by Lions [26]. The idea is to consider the convection-diffusion equation by adding the diffusion term ϵ∇2v, where ϵ>0 is the viscosity. Then, the solution of this equation becomes smooth and converges to the viscosity solution of original HJB equation as ϵ⟶0 [30]. It can be shown that the existence and uniqueness of viscosity solution to HJB equation are guaranteed [30, 31]. However, this equation generally is not analytically solvable and so we use the upwind finite-difference method as a numerical approximation method.Before addressing this method, we settle for computational problems that arise from using the finite-difference method. Since the variables of model (1) have wide ranges of numbers, solving a feedback control problem requires a high amount of numerical computations. So, by scaling the variables to the appropriate level, we transformed model (1) to dimensionless model (21). Assuming that E0≤E∗ and T0≤T∗, it is followed from Theorem 2 that xt=x1t,x2t∈0,1×0,1, for all t>0. Therefore, we restrict the feedback control problem to the domain 0,tf×Ω, where Ω=0,1×0,1. Then, we can reformulate equation (26) as follows:(32)ε∇2v+vt+minu∈Uf∇xvt,x⋅ft,x,u+Lt,x,u=0,∀t,x∈0,tf×Ω,vt,x=gt,x,∀t,x∈0,tf×∂Ω,vtf,x=0,∀x∈Ω.We need to be determine a suitable boundary condition, fort,x∈0,tf×∂Ω. Since the only possible information is the value function at the final time tf, we pose the Dirichlet boundary condition, that is, gt,x=vtf,x [32].Here, we apply the upwind finite-difference method to (32) with ϵ≪1 [33, 34]. Let i=i1,i2,i1+=i1+1,i2,i1−=i1−1,i2,i2+=i1,i2+1, and i2−=i1,i2−1. To discretize (32), we assume that, for k=1,2, xk,ik=ikΔxk, ik=0,1,…,Mk, with Δxk=1/Mk, and tj=j−1Δt, j=1,2,…,Mt+1, with Δt=tf/Mt, are the partitions along the xk and t axes, respectively. Here, we assume the meshes for state variables are uniform, that is, M1=M2=M. Then, two stencils can be considered. One consists of the six vertices x1,i1−1,tj, x1,i1,tj, x1,i1+1,tj, x1,i1−1,tj+1, x1,i1,tj+1, and x1,i1+1,tj+1. The other consists of the six vertices x2,i2−1,tj, x2,i2,tj, x2,i2+1,tj, x2,i2−1,tj+1, x2,i2,tj+1, and x2,i2+1,tj+1. On the cuboid consisting these stencils, (32) can be approximated by the difference equations as follows; for i1,i2=j,…,M−j and j=1,2,…,Mt,(33)vij+1−vijΔt=∑k=121+sgnfk,ijuij2fk,ijuijvik+j−vijΔxk+1−sgnfk,ijuij2fk,ijuijvij−vik−jΔxk+Lijuij,uij=argsupu∈Uf−∑k=121+sgnfk,iju2fk,ijuvik+j−vijΔxk+1−sgnfk,iju2fk,ijuvij−vik−jΔxk−Liju,where(34)f1,ijuij=f1tj,x1,i1,x2,i2,uij,f2,ijuij=f2tj,x1,i1,x2,i2,uij,Lijuij=Ltj,x1,i1,x2,i2,uij.Based on the upwind finite-difference method, the algorithm for solving HJB equation is summarized in Algorithm1.Algorithm 1: Algorithm for solving HJB equation
Set the initial conditionvi1ik=0Mk and the feedback control ui1ik=0Mk for k=1,2forj=1,2,…,Mtdofori1=j,j+1,…,M1−jdofori2=j,j+1,…,M2−jdoWiju=−∑k=121+sgnfk,iju/2fk,ijuvik+j−vij/Δxk+1−sgnfk,iju/2fk,ijuvij−vik−j/Δxk−Lijuuij=argsupu∈UfWijuvij+1=vij−ΔtWijuijendendend


## 4.2. Open-Loop Control Problem

In open-loop control problem,tf and u are the same as meanings in the feedback control problem. Then, an objective functional for open-loop control problem is given by(35)Gu=∫0tfLt,xt,utdt,where x0=x10,x20 is an appropriate initial condition. Similarly, the open-loop control problem is to find an optimal control u∗ such that(36)Gu∗=minGu,u∈Uo.subject to system (23), where the control set is defined as(37)Uo=u|utis Lebesgue measurable,0≤ut≤1∀t∈0,tf.First, an existence of optimal control for problem (36) and the optimality system is derived.Theorem 8.
Given the objective functionalGu and the control set Uo, there exist an optimal control u∗ such that Gu∗=minGu,u∈Uo.Proof.
From [35], we prove the existence of an optimal control. By Theorem 1, the state variables are nonnegative. Note that the control is also nonnegative. The set of all controls uinUo is convex and closed. For the minimizing problem, the convexity of G in u is satisfied and the optimality system is bounded. So, the compactness necessary for the existence of optimal control is guaranteed. Furthermore, the Lagrangian L is convex on the control set Uo. Since the state variables are bounded, we can easily confirm that there exist a constant κ>1 and numbers ρ1,ρ2 such that(38)Gu≥ρ1u2κ/2−ρ2=ρ1uκ−ρ2.
If follows the existence of an optimal control.
In order to find a minimal value of the Lagrangian, we define a Hamiltonian:(39)H=x2+qEux1−qTux2+B2u2+λ1σ¯+c¯x2−δx1−qEux1+λ2αx21−x2−ω¯x1x2−qTux2,where λ1 and λ2 are the adjoint variables. To derive a necessary conditions for the optimal solution of optimal control problem (23) and (35), we use Pontryagin’s Maximum Principle [24]. The adjoint system and control characterization are described in the following theorem.Theorem 9.
Given an optimal controlu∗ and a solution x∗=x1∗,x2∗, there exist adjoint variables λ1 and λ2 satisfying(40)λ1′t=−∂H∂x1=qEut+λ1tδ+qEut+λ2tω¯x2t,λ2′t=−∂H∂x2=−1+qTut−c¯λ1t−λ2tα−2αx2t−ω¯x1t−qTut.under transversality conditions,(41)λitf=0,i=1,2.
Moreover, the control functionu∗ is given by(42)u∗=min1,max0,R,where R=λ1−1qEx1+λ2+1qTx2/B.Proof.
From Pontryagin’s Maximum Principle [24], the adjoint system can be obtained:(43)λ1′t=−∂H∂x1,λ2′t=−∂H∂x2.with λitf=0 for i=1,2. By solving the equation(44)∂H∂u=0and using properties of the control set Uo, we can characterize the optimal control (42).


## 5. Numerical Results

In numerical simulations of feedback control and open-loop control problems, we set the case where the criterionK is greater than 1, i.e., tumor survive persistently. The parameter values are summarized in Table 2. Because chemotherapy usually proceeds only for the initial about 5 days due to toxicity, we assume that model (23) is defined over a time interval [0, 5]. Also, we set weight constant B=1. In particular, for the feedback control problem, we solve the HJB equation numerically by using mesh size for state variables and time with Δx1=Δx2=0.02 and Δt=0.01, respectively.Table 2
Parameter values for model (1).
SymbolValueReferencesσ1.3×104[7]c0.103[12]δ0.0412[7]α0.18[7]β2.0×10−9[7]ω1.101×10−7[7]qE—qT0.9[19]We consider two scenarios forqE. One is qE=0.2, which is a case of lower side effects. Upper panels in Figure 2 represent values of feedback control ut,x1,x2 on the domain Ω at t= 0 and 2 with qE=0.2. At t=0, there are higher control values if effector cells are low and tumor cells are high. The controls are not required in the opposite case where tumor cells are high and effector cells are low. When the level of effector cells is fixed, the controls have higher values which decrease as the number of tumor cells increases. On the contrary, the control values become lower as tumor cells go higher. At t=2, the controls are needed highest when the level of tumor cells is more than about 0.4, which suggests that, in this case, the chemotherapy should depend on the level of tumor cells rather than the effector cells. The other is the case where qE=0.6, which is represented in the lower panels in Figure 2. Similar to the case of qE=0.2, the higher degree of chemotherapy is required when the effector cells are low and the tumor cells are high. Compared to the case of lower side effects, the chemotherapy should be conducted when the immunity is much higher.Figure 2
Feedback controls for initial time valuesτ=0 and 2 with qE=0.2 (upper panels) and with qE=0.6 (lower panels).Furthermore, in order to investigate the utility of chemotherapy (UOC), we consider a measure, UOC, as the function ofqE which is defined by(45)UOCqE=∬ΩuqE0,x1,x2dx1dx2−∬Ωu10,x1,x2dx1dx2∬Ωu00,x1,x2dx1dx2−∬Ωu10,x1,x2dx1dx2,where u0 and u1 represent feedback controls in the case of qE=0 and qE=1, respectively. The numerator describes the difference between the level of chemotherapy adapted when the degree of side effects is qE and in the case of 100% side effects. The denominator indicates the difference between the level of chemotherapy adapted in the absence of side effects and the level of chemotherapy adapted when there are fully side effects. That is, UOC represents the ratio of efficacy of chemotherapy according to the degree of side effects. Figure 3 shows UOC as the chemotherapy-induced death coefficient of effector cell, qE, varies. UOC decreases from 1 to 0 as qE increases and falls below 50% if qE is higher than 0.4. A case with qE=0.2 has UOC as 70%, but it is about 30% when qE is 0.6.Figure 3
Utility of chemotherapy (UOC) for the chemotherapy-induced death coefficient of effector cellqE.In addition, we examine the efficacy of feedback control for chemotherapy as compared to the open-loop control. This is shown in Figure4. In the case of qE=0.2, the open-loop control initially reduces tumor cells slightly more than the feedback control, but, at the same time, the effector cells decrease more. Through the feedback control, it can be confirmed that effector cells survived more and the tumor cells were reduced more at the final time. However, in this case, there was no big difference between feedback control and open-loop control. Especially, considering the higher side effects, the difference became higher.Figure 4
Numerical simulations of tumor-immune model (1) with chemotherapy. The blue line means the feedback control and the green line represents the open-loop control. The upper panels indicate the number of effector cells for qE=0.2 (a) and qE=0.6 (b), respectively. The lower panels indicate the number of tumor cells for qE=0.2 (c) and qE=0.6 (d), respectively. The initial values are assumed to be E0,T0=5.0×106,5.0×108.
(a)(b)(c)(d)


## 6. Discussion and Conclusion

Chemotherapy is known as the common cancer treatment, which can directly kill tumor cells. Since there are side effects of killing effector cells and alternative treatment such as immunotherapy have been introduced, the use of chemotherapy has decreased. However, due to the advantages of time and cost, it is still the most widely used. For the effective chemotherapy, we should take into account for not only the size of tumor cells but also the level of effector cells. Thus, in this study, we set up the mathematical model for tumor-immune system and considered feedback control problem using HJB equation to show the effective chemotherapy with simulation. The upwind finite-difference method was used to solve the problem numerically. As a result, the optimal chemotherapy was designed depending on the level of effector cells as well as the size of tumor cells. Also, we calculated the degree of chemotherapy according to the extent of side effects and then suggested the utility of chemotherapy (UOC).Furthermore, considering the open-loop control for chemotherapy through PMP, we compared the results of feedback control with that. In the case of feedback control, the effector cells were reduced less and the tumor cells were decreased more at the same extent of side effects. Consequently, the feedback control could suggest the effective chemotherapy compared to the open-loop control. PMP is widely used as a conventional optimization method, but it is only a necessary condition for an optimum, which is just a local solution for optimal control problems. On the contrary, the dynamic programming can guarantee necessary and sufficient conditions for an optimum. However, it has also drawbacks of complicated process and a lot of computation times. We were able to reduce the complexity and computation times by using the simple model for tumor-immune system, but this is one of the limitations in our study. Therefore, in this regard, we plan two points as future works. First is to investigate more sophisticated biological mechanisms for the tumor-immune system and then establish a mathematical model for those. Secondly, we will improve the existing numerical method to effectively reduce the computation time of dynamic programming, which is expected to make it easier to apply mathematical models to feedback control problems.



---

*Source: 2158052-2022-04-13.xml*