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.

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 Section 2, 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 Figure 1. The model is given by the following differential equations:

Figure 1 

Schematic diagram of model (1) for interaction of tumor cells and immune cells (commonly called effector cells).

Tumor cells 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 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 . The explanation of the parameters used in model (1) is in Table 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. Let be a solution of model (1). If the initial conditions satisfy and , then the solution will remain positive, for all .

Proof. Assume that there exists . First, we consider the case of . Then, we know that and , for all . It follows that, for all ,By Gronwall’s inequality, we obtain the inequality as follows. For all ,This is a contradiction with the fact that .
Now, suppose that and note that , for all . It is followed from the extreme value theorem that has positive maximum value . Then, it is obtained that, for all ,By the partial fraction decomposition, we obtain the inequality as follows; for all ,where . This contradicts the assumption of . Thus, will remain positive.

Theorem 2. Let be a solution of model (1). If and , then and are bounded.

Proof. Suppose is unbounded. First, we consider the case of . Let . Then, it follows that . However, since is unbounded, there exists such that , which also implies that . By the recursive procedures, we obtain that for . This contradicts the assumption. So, , for all . Now, consider . Let . Then, there exists such that . It follows that . The remaining proof is similar to the first case. Thus, , for all . Therefore, is bounded. Then, we can obtain the following inequality:where . It is followed thatTaking the of both sides, we obtainTherefore, 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,

To check the local stability, we compute the Jacobian at as follows:

If the trace is negative and the determinant is positive, then is locally asymptotically stable. Thus, the condition should be satisfied. We define the criterion as follows:

The local stability of is summarized in the following theorem.

Theorem 3. If , model (1) has an tumor-free equilibrium which is locally asymptotically stable.
The other equilibrium is , where

If , the positiveness of is guaranteed. Note that if , the equilibrium does not exist and when . The local stability of the equilibrium is summarized in the next theorem.

Theorem 4. If , model (1) has an unstable equilibrium and a locally asymptotically stable positive equilibrium .

Proof. Computing the Jacobian at givesTo show that is locally asymptotically stable, we compute the following conditions:Since , the trace is negative and the determinant is positive. Therefore, the equilibrium is locally asymptotically stable. In this case, the tumor-free equilibrium is unstable, which is a saddle.

3.2. Global Stability

Furthermore, we investigate the global stability of the equilibria and . These are summarized in the following theorems.

Theorem 5. If , the tumor-free equilibrium of model (1) is globally asymptotically stable.

Proof. Consider the following Lyapunov function on :It is easy to see that at the equilibrium and , for all . Thus, is positively definite and is also clearly radially unbounded. Taking the derivative of with respect to , we obtain thatBy the arithmetic-geometric mean and the assumption , we obtain the following inequality:Hence, we have , for all . Therefore, by Lyapunov’s theorem, the tumor-free equilibrium is globally asymptotically stable.

Theorem 6. Assume ; the equilibrium of model (1) is globally asymptotically stable whenever .

Proof. Let and set the functions:Define the function on . Then, is continuously differentiable on . Furthermore, we obtain thatApplying Dulac–Bendixson criterion, system (1) has no periodic orbits or graphics on . Since Theorem 2 guaranteed the boundedness of solutions of system (1), is positively invariant and contains an omega limit set of every initial conditions. Thus, Poincaré–Bendixson theorem can be applied. When , if , then the solutions will remain on the -axis in -plane and converge to . If , that is, , then we want to show that the equilibrium does not belong to the omega limit set of the point . Suppose the omega limit set contains . By Theorem 4, is unstable saddle, and so it has stable manifold which is given by the -axis. Because each solution starting from remains on the -axis and converges to , thus, the stable manifold of is not in . Hence, the omega limit set should contain another equilibrium . But, since is locally asymptotically stable, the solutions starting close to , approach to it. So, the equilibrium is not in the omega limit set. Therefore, the omega limit set of contains only . □

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:where and . Then, the new dependent variables and are dimensionless quantities. Hence, the system for them becomeswhere , , and .

4.1. Feedback Control Problem

4.1.1. Hamilton–Jacobi–Bellman Equation

Let and be an initial time and a final time, respectively, and let be a measurable function from to . The control means the treatment effort in time , which represents the intensity of chemotherapy. The class of all admissible controls is denoted by

Then, the system with the control is as follows:where and are the chemotherapy-induced death coefficients of effector cell and tumor cell, respectively. Especially, the coefficient 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:where the Lagrangian and is an initial condition. In the Lagrangian , 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:

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:where and the initial condition . The value function in characterizing the optimal feedback law is described in the following two propositions [29].

Proposition 1. Let be the value function. Then, if there exists a control such thatthen is an optimal control, where is the state corresponding to .

The optimal control is denoted asfor almost all .

Proposition 2. Let be the value function. Then, is an optimal control-trajectory pair in feedback form if and only iffor almost all .
By Propositions 1 and 2, the feedback law via the value function can be constructed. This is summarized in the following theorem.

Theorem 7. Let be the value function. Suppose satisfiesThen, is the feedback law of the optimal control problem, where is the solution of the following control system:

4.1.2. Upwind Finite-Difference Method

Solutions of HJB equation are continuous, but nonsmooth in . 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 , where is the viscosity. Then, the solution of this equation becomes smooth and converges to the viscosity solution of original HJB equation as [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 and , it is followed from Theorem 2 that , for all . Therefore, we restrict the feedback control problem to the domain , where . Then, we can reformulate equation (26) as follows:

We need to be determine a suitable boundary condition, for . Since the only possible information is the value function at the final time , we pose the Dirichlet boundary condition, that is, [32].

Here, we apply the upwind finite-difference method to (32) with [33, 34]. Let , and . To discretize (32), we assume that, for , , , with , and , , with , are the partitions along the and axes, respectively. Here, we assume the meshes for state variables are uniform, that is, . Then, two stencils can be considered. One consists of the six vertices , , , , , and . The other consists of the six vertices , , , , , and . On the cuboid consisting these stencils, (32) can be approximated by the difference equations as follows; for and ,where