Optimal control of effector-tumor-normal cells dynamics in presence of adoptive immunotherapy

1.   Introduction

Cancer is still considered as one of the dreaded diseases among the human population. The majority of types of cancer precipitate through the growth of malignant tumor cells. So, special care is needed to control the growth of malignant tumors ^[2]. Effective treatments for cancer patients include surgery, chemotherapy, immunotherapy, and radiation therapy. There are a variety of immunosuppressive therapies including cell transfer, dendritic cells, vaccines, monoclonal antibodies, cytokines, and adjuvant therapy to name a few ^[3,4,5]. In cancer, an increase in the number of growing tumor cells kills the tissues of our body ^[6]. The response of a tumor to treatment depends on several factors including the severity of the tumor, the patient's response to treatment, and the location of the tumor, etc. In some cases, immunotherapy is very effective and such a treatment stimulates the immune system to fight the tumor cells and to have fewer side effects than other therapies ^[3,4,7,8]. Therefore, our goal is not only to stop cancer but also to explore effective treatment strategies to eradicate cancer without putting the patients' health condition to further risk.

Adoptive immunotherapy targets the growth and development of antibodies by increasing the number of effectors or immune cells that help our immune system to fight cancer and other diseases. It is made up of white blood cells, tissues, and the lymphatic system. Adoptive immunotherapy, which acquires implants in the immune system to fight cancer, enhances the body's natural ability to fight tumor cells. This treatment consists of tumor antigens on cancer cells, whose molecules are identified on the face and bound with antibody proteins. Tumor antigens are usually proteins or other macronutrients (e.g., carbohydrates) and common antibodies that bind to foreign viruses, but anti-immunotherapy antibodies bind to antigen markers and target cancer cells so that the immune system blocks or kills tumors ^[9]. The mathematical model of tumor-immune interactions is discussed by Qomlaqi et al. ^[5] which can be used to develop a systematic approach to immunotherapy treatment. Interested readers are referred to the papers ^{[8,10,11,12,13,14,15]} to have more information about the use of adoptive immunotherapy in cancer treatment.

Mathematical modeling and optimal control theory have played an important role in answering such important research questions which are found generally through experimentation, but the performance of the required experiments becomes too costly or otherwise difficult to perform. Kuznetsov et al. ^[16] described the cytotoxic $T$-lymphocyte response to tumor development. Kuznetsov and Knott ^[7] discussed the methods of tumor growth, its suppression, and regeneration. De Pillis and Raduskaya ^[1] presented the competition between normal cells and tumor cells in a model considering the role of chemotherapy into account. De Pillis et al. ^[17] introduced a mathematical model in tumor development using combined immunotherapy and chemotherapy. Nagata et al. ^[18] discussed the dynamical behaviours of $T$ cell tumor response under the support of dendritic cells. Many researchers have worked on tumor growth models with optimal control strategies. Optimal control strategies are very helpful in finding effective treatment strategies so that the dose of the drug can be controlled with time to lessen different side effects which otherwise may put the patients' health under threat. Interested readers are referred to the papers ^{[2,5,6,9,19,20,21]}.

In modern times, a variety of optimal control methods have been introduced to find the best treatment strategy which can cure cancer, inflicting minimum side effects ^{[6,22,23,24,25]}. Khajanchi and Ghosh ^[6] demonstrated the method of combined optimal control and developed strategies that could increase the number of effector cells and minimize the overall dose of the drug that could eradicate tumor cell load. Nastitie and Khusnul Arif ^[22] discussed the cancer treatment process through a combination of radiotherapy and anti-angiogenic therapy aimed at reducing cancer size by adding control doses of radiotherapy and anti-angiogenic. Glick and Mastroberradino ^[23] used anti-angiogenesis doses at the end of two weeks of treatment to determine the optimal dose to reduce tumor size and reduce toxicity to the patient. Bukkuri ^[24] carried out optimal control analysis of combined chemotherapy-immunotherapy in the case of a cancer model. Irana Ira et al. ^[25] discussed the dynamics of tumor cell growth and its optimal control. It is important to note the fact that the application of appropriate control theory requires a variety of populations which ultimately increases the number of equations in the mathematical model under consideration making it a difficult problem to handle with. To have an idea in this context, interested readers are referred to the papers ^{[26,27,28,29]}.

In the present paper, we have discussed about a mathematical model having a cancerous cell in the form of tumors and have adopted the method of adoptive immunotherapy to get rid of the cancerous tumor cells. To make the treatment strategy a better one in which the patient's health condition does not come under threat due to different side effects, we have converted it to an optimal control problem and solved it.

The rest of this paper is as follows: In Section 2, we illustrate the model formulation. Positiveness and boundedness of the solutions of the model are verified in Section 3. The existence of equilibrium points and the local stability of those are discussed in Sections 4 and 5, respectively. In Section 6, the Lyapunov function is constructed to investigate the global stability of a locally stable healthy equilibrium point. In Section 7, we discuss the problem of optimal control concerning the model using Pontryagin's minimum principle and the Hamiltonian function. Numerical simulations and corresponding discussions are presented in Section 8. Finally, the paper ends with a detailed review in Section 9.

2.   Model formulation

Our considered model in this paper is a slight modification of the model proposed by De Pillis et al. ^[1]. Modifications are done in the first and second equations of the model due to the incorporation of adoptive immunotherapy. In the considered model, $E\left(t\right)$, $T\left(t\right)$, and $H\left(t\right)$ denote the total density of effector or immune cells, tumor cells, and healthy or normal cells at any time $t>0$, respectively. The system of differential equations describing the growth, death, and interactions of this population is given by

with initial conditions $E\left(0\right) = {E}_{0}$, $T\left(0\right) = {T}_{0}$, and $H\left(0\right) = {H}_{0}$ where each initial value is positive.

In the first equation of model (2.1), the effector cells have a constant source rate $\mu$, while the density of effector cells is proportional to death using the term ${d}_{1}E$. Response of tumor-specific effector cells is regulated by the Michaelis-Menten form $\rho ET/\left(\sigma +T\right)$ which gives a saturation effect. The third term ${c}_{1}ET$ represents the decomposing rate of effector cells by tumor cells. The last term ${v}_{1}$ represents the external source of adoptive cellular immunotherapy. This last term has not been applied in the model proposed by De Pillis et al. ^[1].

In the second equation of model (2.1), the tumor cell population is logistically increasing, as justified in ^[6,7,8,9], and is killed by the immune cells and normal cells through a mass-action dynamics which are represented by the terms ${c}_{2}ET$ and ${c}_{3}TH$. Tumor cell growth rate is ${r}_{1}$ and ${b}_{1}$ is the maximum carrying capacity of tumor cells. The final term ${v}_{2}T$ refers to tumor cells killed by an external injection of adoptive immunotherapy. This last term has not been used in the model introduced by De Pillis et al. ^[1].

In the third equation of model (2.1), the normal cell population is also growing logistically, ${r}_{2}$ is representing the growth rates of normal cells, and the maximum carrying capacities of normal cells is one ^[1]. The second term ${c}_{4}TH$ refers to the rate of normal cell death owing to tumor cells.

3.   Positive invariance and boundedness

Before we proceed with the mathematical analysis, we must show that the model is biologically feasible and that the values of all parameters are non-negative. According to the standard comparison theory, it follows

Integration of the above leads to

Again

Proceeding as above, we have

and similarly, we find

Thus, the feasible region for model (2.1) to be realistic is defined as $\psi = \left\{\right(E, T, H)\in {R}_{+}^{3}\}$.

We further assume that the initial values are positive i.e. $E\left(0\right)\ge 0$, $T\left(0\right)\ge 0$, and $H\left(0\right)\ge 0$, then $E\left(t\right)\ge 0$, $T\left(t\right)\ge 0$, and $H\left(t\right)\ge 0$ for all $t\ge 0$. The trajectories evolve in the attracting regions

The domain $\psi$ is positive invariant for model (2.1) and thus biologically meaningful for the cell densities. We will analyze the model quantitative behaviour in the domain $\psi$ ^[30]. This verifies that model (2.1) is biologically feasible.

4.   Existing of equilibrium points

We consider a case where the external source of adoptive cellular immunotherapy and tumor cells killed by an external injection of adoptive immunotherapy per day is a constant parameter ^[6].

Equilibrium points are found by

The simplification of the above results shows that

(ⅰ) ${S}_{1}\left({E}_{1}, {T}_{1}, {H}_{1}\right)$ is a healthy equilibrium (where there is no existence of tumor) point where ${E}_{1} = \left(\mu +{v}_{1}\right)/{d}_{1}$, ${T}_{1} = 0$, and ${H}_{1} = 1.$

(ⅱ) ${S}_{2}({E}_{2}, {T}_{2}{, H}_{2})$is a co-existing or unhealthy equilibrium point (where there is the existence of tumor) where

and ${T}_{2}$ can be found from the solution of the equation

or

where

Note: Equilibrium points are valid if ${r}_{2}>{c}_{4}{T}_{2}$ and $({d}_{1}+{c}_{1}{T}_{2})(\sigma +{T}_{2})>\rho {T}_{2}$ which can be seen easily from the above equations.

5.   Local stability analysis

Without treatment case $({\mathit{\boldsymbol{v}}_1} = {\mathit{\boldsymbol{v}}_2} = {\bf{0}})$: In this section, we study the nature of stability of the equilibrium points ${S}_{1}$ and ${S}_{2}$ by considering ${v}_{1}=0$ and ${v}_{2}=0$; these equilibrium points are

and

The Jacobian matrix of the system (2.1) is

where

(ⅰ) At the first equilibrium point, the eigenvalues of the Jacobian matrix ${J}_{S*}$ are

By applying the condition of stability, the necessary condition for asymptotic stability of equilibrium point ${S}_{1}^{*}$ is found to be

and it will be unstable when ${r}_{1}>\frac{{c}_{2}\mu }{{d}_{1}}+{c}_{3}$.

(ⅱ) The eigenvalues associated with co-existing equilibrium point ${S}_{2}^{*}\left({E}_{2}^{*}, {T}_{2}^{*}, {H}_{2}^{*}\right)$ are derived from the Jacobian matrix ${J}_{S*}$. The characteristic equation is

where (in (5.1))

and

From here, we get

and

By Routh-Hurwitz stability criteria, if ${A}_{11}>0$ and ${A}_{11}{A}_{12}-{A}_{13}>0$, then ${S}_{2}^{*}$ is stable and becomes unstable when conditions are not satisfied. Validity of (5.2) and (5.3) are verified by putting the parameter values from Table 1.

With treatment case (${\mathit{\boldsymbol{v}}_1} \ne {\bf{0}}\; {\bf{and}}\; {\mathit{\boldsymbol{v}}_2} \ne {\bf{0}}$): In this section, we study the nature of stability of the equilibrium points ${S}_{1}$ and ${S}_{2}$ of the system and here ${v}_{1}\ne 0$ and ${v}_{2}\ne 0$.

The Jacobian matrix of the system (2.1) is

where

(ⅰ) At the first equilibrium point, the eigenvalues of the Jacobian matrix ${J}_{S}$ are

By applying the condition of stability, the necessary condition for asymptotic stability of equilibrium point ${S}_{1}$ is found to be

and it will be unstable when ${r}_{1}>\frac{{c}_{2}(\mu +{v}_{1})}{{d}_{1}}+{c}_{3}+{v}_{2}$.

(ⅱ) The eigenvalues associated with co-existing equilibrium point ${S}_{2}\left({E}_{2}, {T}_{2}, {H}_{2}\right)$ are derived from the Jacobian matrix ${J}_{S}$. The characteristic equation is

where (in (5.4))

and

From here, we get

and

By Routh-Hurwitz stability criteria, if ${B}_{11}>0$ and ${B}_{11}{B}_{12}-{B}_{13}>0$, then ${S}_{2}$ is stable and becomes unstable when conditions are not satisfied. Validity of (5.5) and (5.6) are verified by putting the parameter values from Table 1.

Now, we conclude that ${S}_{1}$ and ${S}_{2}$ as equilibrium points are biologically feasible and stable under the condition of biological existence.

6.   Global stability analysis of the healthy equilibrium point ${\mathit{S}}_{1}$ in treatment case

For the behaviour of the system (2.1) far away from the equilibrium point ${S}_{1}$, we analyze the global stability of ${S}_{1}$ in this section. Let's define the Lyapunov function of the model (2.1) as

Now, we differentiate w.r.t. time to obtain

where

By noting the second component of the vector $V$ in (6.1), we must have:

where such a condition, namely (6.2), results in ${V}^{T}Y>0.$Furthermore, by considering the values of parameters from Table 1, if $E = \left(\mu +{v}_{1}\right)/{d}_{1}$ and $T = 1/{b}_{1}$, then all minors of the matrix $M$ are positive (all eigenvalues of $M$ are also positive) and so ${Y}^{T}MY>0.$ Now, it is clear that

Therefore, the healthy equilibrium point ${S}_{1}$ is globally asymptotically stable if

In biological terms, it means that the tumor cells will be killed by an external injection of the adoptive immunotherapy if

The vector field plot along with four trajectories with different initial points (with the same set of parameter values taken in Table 1) and the equilibrium point ${S}_{1}$ has been shown in Figure 1.

From the vector field plot, it is seen that any trajectory starting from any starting point in the basin of attraction converges to the tumor-free equilibrium point ${S}_{1}$ (indicated by the black dot) indicating that it is a globally stable point for the system. Biologically, this indicates the fact that the body is recovering from the tumor regardless of the initial condition which contains tumor growth.

7.   Optimal control

This section is devoted to the study of the model after we administer adoptive immunotherapy treatment at a specific time. From a biomedical perspective, we have used the concept of optimal control in the model. In the earlier discussed case, the amount of injected adoptive immunotherapy remains the same even when the tumor size gets reduced. This in turn might infuse detrimental effects to the patients' immune system and other diseases may attack the patient. So, practically, the amount of injected immunotherapy should be decreased when the size of the tumor gets smaller. For this purpose, we should look into the problem with a control strategy that can lessen the health hazard for the patient. Therefore, we propose and analyze the optimal control problem applicable to model (2.1) to determine the optimal dose of adoptive immunotherapy to control the tumor. We determine control inputs ${v}_{1}$ and ${v}_{2}$ of cellular immunotherapy which are included in the first and second equations of the model (2.1); to be supplied from an external source at different times.

So, let us assume that the time-dependent form of our considered model has the following form:

with the initial conditions

The objective function which is to be minimized is defined as:

where ${\alpha }_{1}$, ${\alpha }_{2}$, ${\beta }_{1}$, and ${\beta }_{2}$ are non-negative constants. It should be mentioned that ${ \beta }_{1}$ and${ \beta }_{2}$ represent the weight factors of the objective function and are used for balancing the size of the terms. The optimal combination of control variables ${v}_{1}$ and ${v}_{2}$ will be adequate to minimize the tumor density and negative side effects over a fixed time. The first two terms of the integrand function are the total number of tumor cells and immune cells. The third and fourth terms of the integrand show the effect of adoptive immunotherapy on the body. Here, we have used the problem of optimal control for the model to reduce the burden of tumor cells, to reduce the time period for recovery of the patient, to reduce side effects due to immunotherapy, and also to increase the effectiveness of immunotherapy to strengthen the immune system.

Here, we establish an optimal control ${v}_{1}^{*}$ and ${v}_{2}^{*}$ such that

where $∆ = \left\{{v}_{1}, {v}_{2}:\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}, \mathrm{ }0\le {v}_{1}, {v}_{2}\le 1, t\in \left[0, {t}_{f}\right]\right\}$ is the admissible control set.

7.1.   Existence of optimal control

In this sub-section, the existence of optimal control of the system (7.1) is discussed. The property of supersolutions $\overline{E}$, $\overline{T}$, and $\overline{H}$ of the model (7.1) is that trajectories are given by

are bounded ^[9]. We rewrite (7.4) as follows:

Since it is a linear system with bounded coefficients and the time limit is limited, we conclude that supersolutions $\overline{E}$, $\overline{T}$, and $\overline{H}$ of the above system are uniformly bounded. Using the theorem proposed by Lukes ^[31], we found that the admissible control class and the corresponding state equations are nonempty with initial conditions given in (7.2). Also, by the definition of the set $∆$, the control set $∆$ is convex and closed. Since the state solutions are bounded, hence, the continuity of R.H.S of the state system (7.1) holds and is bounded above by a sum of the bounded control and state.

Now, we have to show that the convexity of the integrand of ${\mathit{\Omega}} \left({v}_{1}, {v}_{2}\right)$ on $∆$ and bounded below by ${\tau }_{1}\left({v}_{1}^{2}+{v}_{2}^{2}\right)-{\tau }_{2}$ with ${\tau }_{1}, {\tau }_{2}>0$.

Let $u$ and $w$ are distinct elements of ${\mathit{\Omega}}$ and $0\le p\le 1.$ We have to show

where

and $u$ and $w$ are two control vectors and $p\in \left(\mathrm{0, 1}\right).$

By substituting (7.7) into (7.6), we get

which implies that

and

This shows that${\tau }_{1}({v}_{1}^{2}+{v}_{2}^{2})-{\tau }_{2}$ is a lower bound of ${\mathit{\Omega}} ({v}_{1}, {v}_{2})$.

Therefore, there exists an optimal control ${v}_{1 }^{*}$ and ${v}_{2 }^{*}$for which ${\mathit{\Omega}} ({v}_{1}, {v}_{2})$ is minimized. From the above analysis, we establish the following theorem.

Theorem 7.1. For given objective functional

where

is subject to the system (7.1) with the initial conditions $E\left(0\right) = {E}_{0}$, ${T}{(0)} = {{T}}_{{0}}$, and ${H}{(0)} = {{H}}_{{0}}$, there exists an optimal control ${v}_{1 }^{*}$and ${v}_{2}^{*}$ such that ${\mathit{\Omega}} \left({v}_{1}^{*}, {v}_{2}^{*}\right) = \mathrm{min}\left\{{\mathit{\Omega}} \left({v}_{1}, {v}_{2}\right):{v}_{1}, {v}_{2}\in ∆\right\}$.

7.2.   Characterization of the optimal control

Now, we implement the procedure of applying the Pontryagin minimum principle and Hamiltonian function. We introduce the three co-state variable ${p}_{s}, s = 1, 2, 3$ and so the Hamiltonian function is given by

By substituting (7.1) and (7.2) into (7.9), we find

The Hamiltonian equations are

where ${p}_{s}\left(t\right), s = 1, 2, 3$ are the adjoint functions to be determined suitably.

Adjoint equations and forms of transversality conditions are standard results from the Pontryagin minimum principle ^[32,33]. In the case of our considered system, an adjoint system can be obtained in the form of:

where ${p}_{s}\left({t}_{f}\right) = 0, s = 1, 2, 3$ are the transversality conditions.

The optimal control functions are determined by the circumstances $\frac{\partial M}{\partial {v}_{i}} = 0, i = \mathrm{1, 2}$. Hence, we get

Using the bounds for the control variables ${v}_{1 }^{*}$and ${v}_{2}^{*}$ from (7.13), we get

In the compact notation, let us consider

and

From (2.1), (7.12), (7.14), and (7.15), we get the subsequent optimal system as

subject to the conditions $E\left(0\right) = {E}_{0}$, $T\left(0\right) = {T}_{0}$, and $H\left(0\right) = {H}_{0}$ and ${p}_{s}\left({t}_{f}\right) = 0, s = \mathrm{1, 2}, 3$.

Theorem 7.2. Considering optimal control variables ${v}_{1 }^{*}$and ${v}_{2}^{*}$ and corresponding state variables ${E}^{*}\left(t\right)$, ${T}^{*}\left(t\right)$, and ${H}^{*}\left(t\right)$, there exist ongoing specific adjoint variables ${p}_{s}\left(t\right), s = 1, 2, 3,$ satisfying the following system:

subject to the transversality conditions ${p}_{s}\left({t}_{f}\right) = 0, s = 1, 2, 3$.

In addition, the following properties hold:

Next, we proceed to numerically solve the proposed model and the optimal control problem.

8.   Numerical simulation

This section is devoted to numerical solutions of the system (2.1), first without drug administration, then with the introduction of external immunotherapy and ultimately optimal control model of (2.1) as defined in (7.1). We consider the parameter values in Table 1. The numerical solutions of the model are found using MATLAB, while the optimal system was solved using a fourth-order Runge-Kutta method.

The optimal system (7.1) is associated with the adjoint system (7.14) and (7.15) and separated boundary conditions at times $t = 0$ and $t = {t}_{f}$. The forward method is used to solve the optimal system (7.1) and the backward method is used to solve the respective adjoint system (7.16) for ${t}_{f} = 50$. The variables associated with optimal systems and in the objective functions have different scales. Hence, they are balanced by choosing weight constants ${\alpha }_{1} = 2$, ${\alpha }_{2} = 0.5$, ${\beta }_{1} = 0.5$, and ${\beta }_{2} = 5$ in the objective function given in (7.3).

Case 1: In Figure 2, we have portrayed the dynamics of effector cells (a), tumor cells (b), and the normal cells (c) in (2.1) arising out of their mutual interaction, when no external therapy (${v}_{1} = {v}_{2} = 0$) is used. It gives us an insight into what happens to the system without treatment and from this we get an idea about designing various optimal control policies.

From Figure 2, we note that without treatment, the density of effector cells, tumor cells, and normal cells goes up and down and achieves their stable state after some time interval. As the tumor cell load gradually increases, effector cells also increase and the normal cells decrease and after some oscillation attain their steady-state and reach the unhealthy equilibrium point${ S}_{2}^{*}$. It is also clear from this that removal of tumor cells is not possible without a treatment strategy such as adoptive immunotherapy.

Case 2: In Figure 3, we have plotted the dynamics of the system (2.1) with adoptive immunotherapy treatment but without optimal control.

From Figure 3, it is seen that when treatment is started, though for a brief period, normal cells initially decrease and tumor cells grow, the pattern changes after that time period and the tumor cell load decrease, but the effector cells and normal cells grow and ultimately settle down at the steady-state which is the tumor-free equilibrium point ${S}_{1}$. This shows that with the incorporation of therapeutic strategies such as adoptive immunotherapy, tumor cells can be eradicated from the body of the patient.

Case 3: In Figure 4, we have portrayed the dynamics of the system (2.1) with adoptive immunotherapy treatment when optimal control is used.

From Figure 4, it is seen that the optimal treatment strategies reduce the burden of tumor cells and increase the number of effector cells without causing damage to normal cells after a certain time of introduction of treatment. From Figure 4(b), we can conclude that the incorporation of optimal control to eradicate the tumor cells is more effective as it makes the patients' body tumor-free in less time without putting the patients' health at any risk. From this perspective, we can deduce from the optimal control diagram (Figures 4 and 5) that all efforts at the onset of the disease to reduce the proliferation of tumor cells should be kept under optimal control. Figures 5(a) and 5(b) further show that control inputs ${v}_{\text{1}}$ and ${v}_{\text{2}}$ of drugs can be reduced rather than keeping those constant with the decrease in the number of tumor cells.

9.   Conclusions

In this paper, we have considered a basic model related to cancer which was proposed by De Pillis et al. showing interactions between effector-tumors-normal cells in the human body. Mathematical analysis of the model demonstrated that if the system is left to itself, cancerous tumor cells cannot be eradicated. Such a result tempted the authors to incorporate external immunotherapy in the form of injection showing that the tumorous cells can be wiped out after a certain interval of timekeeping external immunotherapy input constant throughout the time. Of course, this work suffers the drawback that an unreasonable amount of external immunotherapy may put the patients' health at risk. Generally, it is expected that after the introduction of immunotherapy, when tumor cells decrease in number after a certain time, the amount of immunotherapy should also be lessened accordingly instead of keeping it fixed. Another important issue that should be kept in mind is that "can we wipe out the cancerous cells in a lesser interval of time?" To answer these issues we implemented an optimal treatment strategy with the introduction of adoptive immunotherapy. Satisfactory answers were found with regard to the above-mentioned important issues. Appropriate numerical simulations were incorporated to justify our conclusions.

Conflict of interest

The authors declare no conflict of interest.