# Optimal Finite Cancer Treatment Duration by Using Mixed Vaccine Therapy and Chemotherapy: State Dependent Riccati Equation Control

Authors: Ali Ghaffari; Mostafa Nazari; Farhad Arab
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

The main objective of this paper is to propose an optimal finite duration treatment method for cancer. A mathematical model is proposed to show the interactions between healthy and cancerous cells in the human body. To extend the existing models, the effect of vaccine therapy and chemotherapy are also added to the model. The equilibrium points and the related local stability are derived and discussed. It is shown that the dynamics of the cancer model must be changed and modified for finite treatment duration. Therefore, the vaccine therapy is used to change the parameters of the system and the chemotherapy is applied for pushing the system to the domain of attraction of the healthy state. For optimal chemotherapy, an optimal control is used based on state dependent Riccati equation (SDRE). It is shown that, in spite of eliminating the treatment, the system approaches the healthy state conditions. The results show that the development of optimal vaccine-chemotherapy protocols for removing tumor cells would be an appropriate strategy in cancer treatment. Also, the present study states that a proper treatment method not only reduces the population of the cancer cells but also changes the dynamics of the cancer.


---

## Body

## 1. Introduction

Modeling and treatment of cancer are the main focus of many researchers worldwide from clinicians, biologists, mathematicians, and control engineers. Cancer mathematical models create an appropriate insight into the behavior of cells in the presence of cancer cells and their interaction with drugs. On the other hand, preparations of such drugs and medical examination have high risk and cost. These illustrate the importance of mathematical and suitable control modeling for cancer treatment. The cancer treatment models, in addition, will enable researchers to forecast and adjust the behavior of the cancerous tumor [1]. The modeling approaches to study disease dynamics include but are not limited to the following: optimization, compartmental, and dynamical system approaches [2]. In this study, we use a dynamical system approach which shows the interaction among cells and drugs. In [3], a review of nonspatial tumor-immune models is presented.In order to avoid the adverse side effects of such drugs and preserve the level of drug dosage, drugs should be used based on a regular program. Different control methods have been used for solving this problem. Using these methods along with optimizing the amount of drugs used yields to effective diminishing of cancer cells [4].Theory of optimal control has been used in modeling of chemotherapy treatment problems. In this problem, the optimal controller is gained by solving a series of differential equations [5]. Currently, many researchers presented mathematical models to simulate the behavior of the drug and its effects on the body [6]. Chemotherapy treatment program was introduced as an optimal control problem by Swan and Vincent in 1977 [7]. In 1990, Swan studied application of optimal control theory in cancer chemotherapy and described great variation among these models [1]. In 2000, Clare and his colleagues introduced several models in the field of application of chemotherapy in the treatment of breast cancer [8]. In 2001, Parker and Doyle performed a thorough review of articles that build the mathematical models of drug delivery and allocated small parts of the cancer optimal chemotherapy [9]. In 2005, Harrold and, in 2009, Harrold and Parker recognized deficiencies and weaknesses in the treatment of chemotherapy in the clinical programs [10, 11]. In 2007, Nanda et al. applied an optimal control model of two-drug chemotherapy for leukemia [12]. In 2011, Shi et al. presented a summary of the optimization models in chemotherapy treatment programs [6]. In 2013, Moradi et al. designed an optimal robust control of drug delivery in cancer chemotherapy [4]. However, the recent studies assumed that the dynamics of the cancer during treatment is time invariant. In other words, the authors considered the effects of therapeutic inputs only on the system states. However, the dynamics of cancer changes during its progression [13]. As an example, wrecking inputs such as external stresses can disable the DNA repair genes [14]. These inputs are able to change the functions of growth-inhibiting signals (TGF-b), regulatory growth signals (TGF-a), and apoptosis (TP53) [13]. Therefore, an effective treatment method should correct these destructive changes in the dynamic behavior of the system.In this study, a system of ordinary differential equations (ODE) is considered to present the interaction between healthy and cancerous cells. The present study extends the existing mathematical model of [15]. Those studies have investigated the effects of therapeutic inputs on the system states. However, the important shortcoming is that the cancer relapses after elimination of the therapy. In this paper, a method for finite duration treatment is proposed such that at the end of treatment the system approaches its healthy equilibrium point. Moreover, the presented model is analyzed by adding vaccine and chemotherapy treatment terms. The vaccine has an effect on some parameters of the system [16], while chemotherapy has an effect on the cells populations. For optimal chemotherapy, the state dependent Riccati equation (SDRE) based optimal control technique is applied to the nonlinear model.The organization of the paper is as follows. In the next section, the nonlinear mathematical model is analyzed. We extend this model by adding vaccine and chemotherapy treatment terms. Then, the SDRE based optimal control is applied to the nonlinear cancer dynamics in Section3. In this work, the amount of chemotherapy drug is considered the control input to the system. The aim of the mixed vaccine and chemotherapy treatments is to present an optimal finite time duration treatment such that the cancer is not able to relapse. In the last section, the simulation results are discussed.The main highlights of the present study can be summarized as follows:(i)
changing the dynamics of the cancer model during the treatment,(ii)
inserting the effects of the vaccine therapy in the cancer model,(iii)
applying SDRE optimal control to the nonlinear cancer dynamics.


## 2. The Mathematical Model

The presented population model originates from [15]. The dynamic behavior of the body organ which is affected by the cancer is given by the following equations:
(1)ddt[xy]=[f1(x,y)f2(x,y)],
where
(2)x˙=a1x(1-xK1)-(d1+c)x-b1xy,y˙=a2y(1-yK2)-d2y+cx-(b2xy+g)y.
In (2), x and y are the healthy and cancer cells concentrations, respectively. The state variables should be physiologically possible; therefore, their values are nonnegative; that is, x≥0 and y≥0. The coefficients a1 and a2 represent the growth rate of the healthy and cancer cells, respectively. The growth rate of healthy and tumor tissues decelerates as the concentrations of both the healthy and tumor tissues approach the carrying capacities K1 and K2, respectively [17]. The effect of the immune system is to kill the mutated and cancer cells at proportional rates d1 and d2. The immune system agents force the cancer cells to suicide through apoptosis [18]. The coefficient c represents the portion of the healthy cells, whose genome is disordered by the external stresses. These cells start the neoplastic transformation and are added to the tumor cells [19]. The tumor competes with healthy tissue for resources, such as blood, nutrients, and space [20]. Moreover, all the cancer cells compete with each other. The competition coefficients between different cells are b1, b2, and g.The aim of this paper is the total recovery of the patient after a finite duration treatment such that the cancer is not able to relapse. In other words, the population of the cancer cells must go to the healthy state after elimination of treatment. However, based on the parameters presented in Table1, the trajectory of the system is shown in Figure 1(a). The equilibrium point 2 in the figure is the only equilibrium point of the system in the first orthant. Based on the stability analysis which will be discussed in Section 2.1, the dynamics of the system represents the cancer state. So, the treatment must be applied during the entire life of the patient. Otherwise, after elimination of the input, the system comes back to the equilibrium point 2. So, for finite duration treatment, the dynamics of the system must be changed. Since the vaccine therapy impacts some parameters of the system, mixed vaccine therapy and chemotherapy treatments are used. The duty of the vaccine therapy is to change the dynamics of the system and the duty of the chemotherapy is to push the system toward the domain of attraction of the healthy equilibrium point.Table 1
Estimated parameters.
Parameters
Estimated value [16, 17, 21, 22]
a
1
3.5
a
2
6.7
K
1
80000
K
2
90000
d
1
0.05
d
2
0.03
c
0.05
b
1
0.2
b
2
0.11
g
0.1
K
x
0.6
K
y
0.6
K
a
1
5
K
b
2
0.4
μ
10The trajectory pattern of the free system: (a) before the vaccination and (b) after the vaccination.
(a)
(b)A complete effect of the vaccine therapy in time duration is needed; therefore, saturation dynamic is applied in the present model to change the parameters of the system. The effect of vaccine therapy is considered on parametersa1 and b2 [16], which is included in the mathematical model by the term vv(t)≥0. The rate of changing these parameters is assumed to be proportional to the input magnitude vv(t). The values of μa1 and μb2 are dependent on the dynamics of a1 and b2, respectively. The biotransformation coefficients saturate at finite limits Ka1and Kb2, which are related to the biological limits of body organs and the accumulation of external effect. The trajectory pattern of the free system after the vaccine therapy treatment is shown in Figure 1(b). As shown in this figure, the healthy state is entered in the first orthant.Also, the effect of chemotherapy is included by the termM(t) for which vM(t)≥0 is the amount of chemotherapy agent injected per day per liter of blood. Some chemotherapeutic drugs, such as doxorubicin, are only effective during certain phases of the cell cycle, and pharmacokinetics also indicates that the effectiveness of chemotherapy is bounded [16]. Therefore, we use a saturation term 1.2M/(0.8+M) to represent the chemotherapy fractional cell kill. Note that the kill rate is almost linear at low concentrations of drug, while it plateaus at a higher drug concentration. Kx is the fractional healthy cell kill by chemotherapy and Kx is the fractional tumor cell kill by chemotherapy [16]. The chemotherapy drug concentration decays in the body and γ is a constant in Day-1 related to the drug elimination rate [16]. So, the modified equations of the system with treatment are as follows:
(3)x˙=a1x(1-xK1)-(d1+c)x-b1xy-Kx(1.2M0.8+M)x,y˙=a2y(1-yK2)-d2y+cx-(b2xy+g)y-Ky(1.2M0.8+M)y,da1dt=μa1vv(t)(1-a1Ka1),db2dt=μb2vv(t)(1-b2Kb2),dMdt=-μM+vM(t).
Note that the system is an autonomous system with differentiable functions, which satisfies existence and uniqueness of initial value problems.Proposition 1.
The subspaceX={(x(t),y(t))∣x≥0,y≥0} is invariant under (2).Proof.
(i) It can easily be concluded that if(x(0),y(0))=(0,0), then (x(t),y(t))=(0,0) for all t≥0.
(ii) If(x(t),y(t)) approaches the vertical axis from subspace X,x˙|x=0=0,y˙|x=0=[a2(1-(y/K2))-(d2+g)]y, and x does not decrease, hence, its value remains on the subspace X.
(iii) If(x(t),y(t)) approaches the horizontal axis from subspace X,y˙|y=0=cx≥0,x˙|y=0=[a1(1-(x/K1))-(d1+c)]x, and y does not decrease, hence, its value remains positive.
Therefore, the subspaceX={(x(t),y(t))∣x≥0,y≥0} is invariant under (2).


### 2.1. Equilibrium Points

Four equilibrium points of (2) are calculated as follows:
(4)(1)x=0,y=0,(5)(2)x=0,y=K2(a2-d2-g)a2=-β2α2,(6)(3)x=α6+α7,y=α1(α6+α7)+β1,(7)(4)x=α6-α7,y=α1(α6-α7)+β1,
where α1=-a1/K1b1, β1=(a1-d1-c)/b1, α2=-a2/K2, β2=a2-d2-g, α3=α1(α1α2-b2), α4=2α1α2β1-b2β1+α1β2+c, α5=β1(α2β1+β2),α6=-α4/2α3, and α7=α42-4α3α5/2α3.The linearization of (2) around the arbitrary equilibrium point (x-,y-) is given by
(8)Z˙*=AZ*+BH.O.T+C,
where
(9)Z*=[xy],A=[-2a1K1x-+b1β1-b1y--b1x-c-b2y--2a2K2y-+a2-d2-g-b2x-],BH.O.T=-[a1K1x2+b1xya2K2y2+b2xy],C=[-a1K1x-2+(a1-d1-c)x--b1x-y--a2K2y-2+(a2-d2-g)y--b2x-y-+cx-].
The higher order terms are neglected around the origin (x,y)=(0,0) and the last term C is equal to zero at the equilibrium points.


### 2.2. Stability Analysis

We study the stability of the equilibrium points of (2) in this section. The results of the analysis are stated as follows.(i)
If the equilibrium point 3 is located in the first orthant and the equilibrium point 2 is not, then the state variables of (2) will converge to the equilibrium point 3 (healthy state).(ii)
If the equilibrium points 2 and 3 are located in the first orthant, then the state variables of (2) will converge to one of these two points.(iii)
If the equilibrium point 2 is located in the first orthant and the equilibrium point 3 is not, then the state variables of (2) will converge to the equilibrium point 2 (cancer state).Proof.
(1) At the equilibrium point 1,x=0,y=0, the eigenvalues of A are
(10)λ1=a1-(d1+c),λ2=a2-(d2+g).
Based on Table 1, the value of a1 is larger than the sum of d1 and c. Also, a2 is larger than the sum of d2 and g. Thus, the eigenvalues of the equilibrium point 1 are always positive and the origin is an unstable node.
(2) At the equilibrium point 2,x=0,y=-β2/α2, the eigenvalues of A are
(11)λ1=b1(β1+β2α2),λ2=-(a2-d2-g).
In (5), y is positive; then β2/α2<0. Also, from Table 1, we notice that the value of -β2/α2 is larger than β1; then the equilibrium point 2, if it exists, is a stable node.
(3) At the equilibrium point 3,x=α6+α7,y=α1(α6+α7)+β1. Noticing that x and y in (6) are in the first orthant, then α6+α7>0 and α1(α6+α7)+β1>0. The principal minors of -A at this equilibrium point are as follows:
(12)Δ1=a1K1(α6+α7)>0,Δ2=(2a12a2K1K2b1[1-1K1(α6+α7)]+a1K1(2b2(α6+α7)+d2+g-2a2(d1+c)K2b1-a2)+b1(c-b2a1-d1-cb1))(α6+α7)>0.
Then, the coefficient matrix A is negative definite at this equilibrium point and it is a stable node.
(4) At the equilibrium point 4,x=α6-α7,y=α1(α6-α7)+β1. Noticing that x and y in (7) are in the first orthant, then α6-α7>0 and α1(α6-α7)+β1>0. The principal minors of -A at this equilibrium point are as follows:
(13)Δ1=a1K1(α6-α7)>0,Δ2=(2a12a2K1K2b1[1-1K1(α6-α7)]+a1K1(2b2(α6-α7)+d2+g-2a2(d1+c)K2b1-a2)+b1(c-b2a1-d1-cb1))(α6-α7)<0.
Therefore, if this equilibrium point exists, then it is an unstable saddle point.


## 2.1. Equilibrium Points

Four equilibrium points of (2) are calculated as follows:
(4)(1)x=0,y=0,(5)(2)x=0,y=K2(a2-d2-g)a2=-β2α2,(6)(3)x=α6+α7,y=α1(α6+α7)+β1,(7)(4)x=α6-α7,y=α1(α6-α7)+β1,
where α1=-a1/K1b1, β1=(a1-d1-c)/b1, α2=-a2/K2, β2=a2-d2-g, α3=α1(α1α2-b2), α4=2α1α2β1-b2β1+α1β2+c, α5=β1(α2β1+β2),α6=-α4/2α3, and α7=α42-4α3α5/2α3.The linearization of (2) around the arbitrary equilibrium point (x-,y-) is given by
(8)Z˙*=AZ*+BH.O.T+C,
where
(9)Z*=[xy],A=[-2a1K1x-+b1β1-b1y--b1x-c-b2y--2a2K2y-+a2-d2-g-b2x-],BH.O.T=-[a1K1x2+b1xya2K2y2+b2xy],C=[-a1K1x-2+(a1-d1-c)x--b1x-y--a2K2y-2+(a2-d2-g)y--b2x-y-+cx-].
The higher order terms are neglected around the origin (x,y)=(0,0) and the last term C is equal to zero at the equilibrium points.


## 2.2. Stability Analysis

We study the stability of the equilibrium points of (2) in this section. The results of the analysis are stated as follows.(i)
If the equilibrium point 3 is located in the first orthant and the equilibrium point 2 is not, then the state variables of (2) will converge to the equilibrium point 3 (healthy state).(ii)
If the equilibrium points 2 and 3 are located in the first orthant, then the state variables of (2) will converge to one of these two points.(iii)
If the equilibrium point 2 is located in the first orthant and the equilibrium point 3 is not, then the state variables of (2) will converge to the equilibrium point 2 (cancer state).Proof.
(1) At the equilibrium point 1,x=0,y=0, the eigenvalues of A are
(10)λ1=a1-(d1+c),λ2=a2-(d2+g).
Based on Table 1, the value of a1 is larger than the sum of d1 and c. Also, a2 is larger than the sum of d2 and g. Thus, the eigenvalues of the equilibrium point 1 are always positive and the origin is an unstable node.
(2) At the equilibrium point 2,x=0,y=-β2/α2, the eigenvalues of A are
(11)λ1=b1(β1+β2α2),λ2=-(a2-d2-g).
In (5), y is positive; then β2/α2<0. Also, from Table 1, we notice that the value of -β2/α2 is larger than β1; then the equilibrium point 2, if it exists, is a stable node.
(3) At the equilibrium point 3,x=α6+α7,y=α1(α6+α7)+β1. Noticing that x and y in (6) are in the first orthant, then α6+α7>0 and α1(α6+α7)+β1>0. The principal minors of -A at this equilibrium point are as follows:
(12)Δ1=a1K1(α6+α7)>0,Δ2=(2a12a2K1K2b1[1-1K1(α6+α7)]+a1K1(2b2(α6+α7)+d2+g-2a2(d1+c)K2b1-a2)+b1(c-b2a1-d1-cb1))(α6+α7)>0.
Then, the coefficient matrix A is negative definite at this equilibrium point and it is a stable node.
(4) At the equilibrium point 4,x=α6-α7,y=α1(α6-α7)+β1. Noticing that x and y in (7) are in the first orthant, then α6-α7>0 and α1(α6-α7)+β1>0. The principal minors of -A at this equilibrium point are as follows:
(13)Δ1=a1K1(α6-α7)>0,Δ2=(2a12a2K1K2b1[1-1K1(α6-α7)]+a1K1(2b2(α6-α7)+d2+g-2a2(d1+c)K2b1-a2)+b1(c-b2a1-d1-cb1))(α6-α7)<0.
Therefore, if this equilibrium point exists, then it is an unstable saddle point.


## 3. Optimal Control for Mixed Drug Administration

A novel optimal control is the SDRE based optimal control. The theoretical background of the “SDRE based optimal control” has not been completely analyzed. This method becomes in the attention of many control engineers due to its computational simplicity and satisfactory experimental results [5]. In this paper, we apply SDRE based optimal control to the nonlinear cancer dynamics. The control input to the model is the amount of chemotherapy drug.


### 3.1. SDRE Optimal Control Theory

Consider the deterministic, infinite horizon nonlinear optimal regulation (stabilization) problem, such that it is full state observable, time invariant, and affine in the input, represented in the following form:(14)x=f(x)+B(x)u(t),x(0)=x0,
where x∈Rn is the state vector, u∈Rm is the input vector, t∈[0,∞) with C1(Rn) functionsf:Rn→Rn and B:Rn→Rn×m, and B(x)≠0 for all x. Without loss of generality, the origin x=0 is assumed to be an equilibrium point. The minimization of the infinite time performance index,
(15)J(x0,u(·))=12∫0∞{xT(t)Q(t)x(t)+uT(t)R(x)u(t)}dt,
is considered, which is nonquadratic in x but quadratic in u. The state and input weighting matrices are assumed to be state dependent such thatQ:Rn→Rn and R:Rn→Rm×m. It is assumed that Q and R are symmetric and R is positive definite:
(16)Q(x)≥0,R(x)>0.
Sincef(0)=0 andf(·)∈C1(Rn), the system (14) can be written in pseudo-linear form:
(17)x˙=A(x)x+B(x)u,
where f(x)=A(x)x. In (17), A(x)∈Rn×n and B(x)∈Rn×m are state dependent coefficient (SDC) matrices which bring the nonlinear system described by (14) into a linear-like representation. These matrices are not unique. However, the recommended selection of the matrices A(x) and B(x) is that they are controllable. The state dependent controllability matrix is as follows:
(18)M(x)=[B(x)A(x)B(x)⋯An-1(x)B(x)].In order to control the nonlinear system, the above matrix must have full rank for the domain for which the nonlinear system is controlled.Some optimal control problems need constraints that must be applied on state variables or the control input. Choice of weight matricesQ(x) andR(x) plays an important role in satisfying these optimal control problems constraints.Hamiltonian matrix for the optimal control problem is as follows:(19)H(x,u,λ)=12(xT(t)Q(t)x(t)+uT(t)R(x)u(t))+λT(A(x)x+B(x)u)-w-T(u-umin⁡)-w^T(umax⁡-u).w-T and w^T are m dimensional nonnegative vectors and are presented to apply constraint to the control input and they must satisfy the following conditions:
(20)w-T(u-umin⁡)=w^T(umax⁡-u)=0.
From the Hamiltonian, the necessary conditions for optimality are
(21)x˙=∂H∂λ=A(x)x+B(x)u,λ˙=-∂H∂x=-Q(x)-[dA(x)xdx]Tλ-[dB(x)udx]Tλ,0=∂H∂u=R(x)u+BT(x)λ-w-+w^.The last equation of (21) gives the optimal control of the following form:
(22)u(x)=-R-1(x)(BT(x)λ-w-+w^).By applying the theory of LQR, the adjoint state vector has the form given by(23)λ=P(x).Finally, the control input is obtained in the following form:(24)u(x)=-R-1(x)BT(x)P(x)x,
which is the control input with the following feedback gain:
(25)K(x)=-R-1(x)BT(x)P(x).P(x) is a symmetric state dependent and positive definite matrix which is given by the solution of algebraic Riccati equations:
(26)AT(x)P(x)+P(x)A(x)-P(x)B(x)R-1(x)BT(x)P(x)+Q(x)=0.Dynamics of the closed loop system is obtained according to the following equation:(27)x˙=(A(x)-B(x)K(x))x.


### 3.2. SDRE Optimal Control Design

The equilibrium point 3 (healthy state) is entered to the first orthant by using vaccine therapy. Now, by using chemotherapy, the trajectory of the system must be pushed toward the domain of attraction of the healthy state. Afterwards, the trajectory of the system approaches the equilibrium point 3 even after eliminating the treatment.In order to design the SDRE based optimal control, we must rewrite the system in the form of (17) by shifting the healthy state to the origin. New state variables are defined as follows:
(28)x1=x-77.89,x2=y-0.158,x3=M.In this case, the system of equations is as follows:(29)x˙1=a1(x1+77.89)(1-x1+77.89K1)-(d1+c)(x1+77.89)-b1(x1+77.89)(x2+0.158)-Kxx3(x1+77.89)1.2+x3,x˙2=a2(x2+0.158)(1-x2+0.158K2)-d2(x2+0.158)+c(x1+77.89)-(b2(x1+77.89)+g)(x2+0.158)-Kyx3(x2+0.158)1.2+x3,x˙3=-μx3+vM(t).
To use the SDRE method, the above equations must be represented in the form of pseudo-linear given by (17). The matrices A(x) and B(x) are
(30)A=[A11-b1(x1+77.89)-Kxx1+77.891.2+x3-b2(x1+0.158)+cA22-Kyx2+0.1581.2+x300-μ],A11=a1(1-x1+77.89K1)-77.89a1K1-(d1+c)-0.158b1,A22=a2(1-x2+0.158K2)-0.158a2K2-d2-g-77.89b2,B=[001]T.The first step to create an optimal control problem is deriving a desirable cost function. We use the following matrix for the cost function:(31)Q(x)=[00001000001],R=4×1012.


## 3.1. SDRE Optimal Control Theory

Consider the deterministic, infinite horizon nonlinear optimal regulation (stabilization) problem, such that it is full state observable, time invariant, and affine in the input, represented in the following form:(14)x=f(x)+B(x)u(t),x(0)=x0,
where x∈Rn is the state vector, u∈Rm is the input vector, t∈[0,∞) with C1(Rn) functionsf:Rn→Rn and B:Rn→Rn×m, and B(x)≠0 for all x. Without loss of generality, the origin x=0 is assumed to be an equilibrium point. The minimization of the infinite time performance index,
(15)J(x0,u(·))=12∫0∞{xT(t)Q(t)x(t)+uT(t)R(x)u(t)}dt,
is considered, which is nonquadratic in x but quadratic in u. The state and input weighting matrices are assumed to be state dependent such thatQ:Rn→Rn and R:Rn→Rm×m. It is assumed that Q and R are symmetric and R is positive definite:
(16)Q(x)≥0,R(x)>0.
Sincef(0)=0 andf(·)∈C1(Rn), the system (14) can be written in pseudo-linear form:
(17)x˙=A(x)x+B(x)u,
where f(x)=A(x)x. In (17), A(x)∈Rn×n and B(x)∈Rn×m are state dependent coefficient (SDC) matrices which bring the nonlinear system described by (14) into a linear-like representation. These matrices are not unique. However, the recommended selection of the matrices A(x) and B(x) is that they are controllable. The state dependent controllability matrix is as follows:
(18)M(x)=[B(x)A(x)B(x)⋯An-1(x)B(x)].In order to control the nonlinear system, the above matrix must have full rank for the domain for which the nonlinear system is controlled.Some optimal control problems need constraints that must be applied on state variables or the control input. Choice of weight matricesQ(x) andR(x) plays an important role in satisfying these optimal control problems constraints.Hamiltonian matrix for the optimal control problem is as follows:(19)H(x,u,λ)=12(xT(t)Q(t)x(t)+uT(t)R(x)u(t))+λT(A(x)x+B(x)u)-w-T(u-umin⁡)-w^T(umax⁡-u).w-T and w^T are m dimensional nonnegative vectors and are presented to apply constraint to the control input and they must satisfy the following conditions:
(20)w-T(u-umin⁡)=w^T(umax⁡-u)=0.
From the Hamiltonian, the necessary conditions for optimality are
(21)x˙=∂H∂λ=A(x)x+B(x)u,λ˙=-∂H∂x=-Q(x)-[dA(x)xdx]Tλ-[dB(x)udx]Tλ,0=∂H∂u=R(x)u+BT(x)λ-w-+w^.The last equation of (21) gives the optimal control of the following form:
(22)u(x)=-R-1(x)(BT(x)λ-w-+w^).By applying the theory of LQR, the adjoint state vector has the form given by(23)λ=P(x).Finally, the control input is obtained in the following form:(24)u(x)=-R-1(x)BT(x)P(x)x,
which is the control input with the following feedback gain:
(25)K(x)=-R-1(x)BT(x)P(x).P(x) is a symmetric state dependent and positive definite matrix which is given by the solution of algebraic Riccati equations:
(26)AT(x)P(x)+P(x)A(x)-P(x)B(x)R-1(x)BT(x)P(x)+Q(x)=0.Dynamics of the closed loop system is obtained according to the following equation:(27)x˙=(A(x)-B(x)K(x))x.


## 3.2. SDRE Optimal Control Design

The equilibrium point 3 (healthy state) is entered to the first orthant by using vaccine therapy. Now, by using chemotherapy, the trajectory of the system must be pushed toward the domain of attraction of the healthy state. Afterwards, the trajectory of the system approaches the equilibrium point 3 even after eliminating the treatment.In order to design the SDRE based optimal control, we must rewrite the system in the form of (17) by shifting the healthy state to the origin. New state variables are defined as follows:
(28)x1=x-77.89,x2=y-0.158,x3=M.In this case, the system of equations is as follows:(29)x˙1=a1(x1+77.89)(1-x1+77.89K1)-(d1+c)(x1+77.89)-b1(x1+77.89)(x2+0.158)-Kxx3(x1+77.89)1.2+x3,x˙2=a2(x2+0.158)(1-x2+0.158K2)-d2(x2+0.158)+c(x1+77.89)-(b2(x1+77.89)+g)(x2+0.158)-Kyx3(x2+0.158)1.2+x3,x˙3=-μx3+vM(t).
To use the SDRE method, the above equations must be represented in the form of pseudo-linear given by (17). The matrices A(x) and B(x) are
(30)A=[A11-b1(x1+77.89)-Kxx1+77.891.2+x3-b2(x1+0.158)+cA22-Kyx2+0.1581.2+x300-μ],A11=a1(1-x1+77.89K1)-77.89a1K1-(d1+c)-0.158b1,A22=a2(1-x2+0.158K2)-0.158a2K2-d2-g-77.89b2,B=[001]T.The first step to create an optimal control problem is deriving a desirable cost function. We use the following matrix for the cost function:(31)Q(x)=[00001000001],R=4×1012.


## 4. Numerical Simulations

In this section, we simulate the behavior of the model by considering the combined treatment. We assume that the initial state of the system is in the domain of attraction of the equilibrium point 2. As shown in Figure2, the chemotherapy treatment pushes the trajectory of the system to the domain of attraction of equilibrium point 3 in an optimal manner. It is created by using vaccine therapy and changing the trajectory pattern of the system in the first orthant. Then, the chemotherapy is stopped and the system approaches the healthy state without any treatment. In other words, at the end of the treatment, the system is placed in a self-destruction cycle by changing the dynamics of the system.The behavior of the system during chemotherapy: (a) the time evolution of the system; (b) trajectory pattern of the system and system response after vaccine therapy.
(a)
(b)The simulation results show that the combined vaccine therapy and chemotherapy treatment is effective for finite duration treatment. In other words, changing the dynamics of the cancer in order to have finite duration treatment is essential. In other words, if the vaccination is not used, the system goes back to its cancer state after the elimination of the treatment (Figure4). In Figure 4, the cancer cells are reduced by chemotherapy, but, after elimination of the treatment, the system approaches the only equilibrium point existing in the first orthant.The behavior of the system during finite duration chemotherapy without vaccine therapy: (a) the time evolution of the system; (b) trajectory pattern of the system and system response without vaccine therapy.
(a)
(b)In [23], the authors proposed on-off regimens for minimizing the number of tumor cells and preserving the healthy cells in an admissible level. However, these suggested regimens are not able to complete elimination of tumor cells. Three different types of quadratic and linear cost functions are considered in [24, 25]. The optimal chemotherapy regimens which are calculated based on these cost functions are able to eradicate the tumor, but they do not consider finite duration treatment method. In [16], de Pillis et al. proposed a mixed immunotherapy and chemotherapy protocol for cancer treatment. The main shortcoming of this protocol is that, after elimination of the treatment, the cancer relapses due to lack of a stable tumor free equilibrium point. In addition, this proposed method is open loop which has many deficiencies such as unrobustness in dealing with parameter variation. In [5, 26, 27], the authors presented SDRE method for cancer treatment. In the model used in those papers, the tumor free equilibrium point is stable. Therefore, changes in the dynamics of the system do not need and the chemotherapy treatment is sufficient for finite duration treatment. However, in the extended version of this model presented in [16], the authors showed that the tumor free equilibrium point is unstable.In the model presented in this paper, it is shown in Figure4 that the chemotherapy alone is not an adequate approach for finite cancer treatment duration. This could be interpreted due to lack of a stable healthy state at the beginning of the treatment (Figure 1). To overcome this problem, the mixed vaccine-chemotherapy is used. In addition, the chemotherapy terms are exerted in a saturation manner, which is in accordance with the physical observations [16].We may conclude that if there is no stable healthy equilibrium point in the model, the dynamics of the system must be changed to reinforce the immune system. This is in accordance with the physical observation. Many evidences exist which show that in some cases the immune system is capable of diminishing the tumor cells without the assistance of external treatments [28].It has to be noted that, in the proposed treatment method, there are also side effects of chemotherapy.


## 5. Conclusion

In this paper, we have modified the existing mathematical models of cancer by mixed vaccine therapy and chemotherapy. We showed that, to obtain the finite duration treatment, a change in the dynamics of the system is necessary. In other words, a suitable cancer treatment method is a method that reduces the population of the tumor cells and also changes the dynamics of the cancer. To change the dynamics of the system, the vaccine therapy is used for changing the parameters of the system and the chemotherapy is also employed for pushing the system to the domain of attraction of the healthy state in an optimal manner. The SDRE optimal control is used for chemotherapy. It is shown that this method has fast and easy derivation for suboptimal control for the chemotherapy problem.



---

*Source: 363109-2014-03-13.xml*