Global Stability and Thermal Optimal Control Strategies for Hyperthermia Treatment of Malignant Tumors

Abstract

Malignant tumor (cancer) is the leading cause of death globally and the annual cost of managing cancer is trillions of dollars. Although, there are established therapies including radiotherapy, chemotherapy and phototherapy for malignant tumors, the hypoxic environment of tumors and poor perfusion act as barriers to these therapies. Hyperthermia takes advantage of oxygen deficiency and irregular perfusion in the tumor environment to destroy malignant cells. Despite successes recorded with hyperthermia, there are concerns with the post-treatment condition of patients as well as the required thermal dose to prevent harm. The investigation of the dynamics of tumor-induced immune suppression with hyperthermia treatment using mathematical analysis and optimal control theory is potentially valuable in the development of hyperthermia treatment. The role of novel tumor-derived cytokines in counterattacking immune cells is considered in this study as a mechanism accounting for the aggressiveness of malignant tumors. Since biological processes are not instantaneous, a discrete time delay is used to model biological processes involved in tumor inhibitory mechanisms by secretion, the elaboration of suppressive cells, and effector cell differentiation to produce suppressive cells. Analytical results obtained using Lyapunov’s function indicate the conditions required for global stability of the tumor-present steady-state. A thermal optimal control strategy is pursued based on optimal control theory, and the best strategy to avoid adverse outcomes is obtained. We validate the analytical results numerically and demonstrate the impact of both inadequate and excessive heat on the dynamics of interactive cell functioning.

1. Introduction

Malignant tumors evolve from the growth of mutated cells which require more energy to survive than normal cells. Moreover, the body’s blood vessels are unable to match the tumor cell’s oxygen and nutrient needs. This phenomenon results in stimulation of the growth of additional blood vessels which are chaotic in structure compared to normal vessels. The insufficiency of oxygen makes the tumor environment hypoxic; however, poor perfusion makes it difficult to kill malignant tumor cells with ionizing radiation or chemotherapy [1]. Hyperthermic treatment takes advantage of oxygen deficiency and irregular perfusion to damage the plasma, the cell skeleton, and the cell nucleus [2] of malignant tumors.

Results from in vitro and in vivo studies have shown that hyperthermia in the range of 40–$44{}^{\circ }$C exerts various effects, including the direct killing of tumor cells, alteration of the tumor microenvironment, induction of heat shock proteins, activation of the immune response, induction of the apoptotic cascade, changing cell cycle regulatory signaling path-ways and alteration of blood flow, oxygen and nutrient distribution in the tumor site, and improvement in therapeutic outcome when applied with other treatments [2,3,4]. These make it successful in treating a wide range of tumors, including tumors of the head and neck, breast, brain, bladder, cervix, rectum, lung, esophagus, vulva and vagina, as well as skin [3].

Hyperthermia has long been thought to suppress the immune system by inducing tolerance [5]. However, more recent studies have demonstrated an increase in immunological attacks against tumors after heat induction, which is suggested to occur through activation of heat-shock proteins and subsequent modulation of innate and adaptive immune responses against tumor cells [5,6,7]. Sufficiently high heating has been found to result in protein denaturation, while milder heating only results in protein inactivation [3,8,9,10]. Despite the results obtained for hyperthermia treatment of malignant tumors, the long-term effect of hyperthermia treatment on patients with malignant tumors, as well as the thermal dose required to prevent adverse consequences, such as the breakdown of DNA and protein denaturation, are major issues of concern in the development of new applications of hyperthermia, such as heat-controlled gene therapy, heat-enhanced immunotherapy and vaccination [11].

The modeling of the tumor immunotherapy process has received significant attention from researchers in recent decades. The dynamics of various approaches to the re-energizing of the immune system have been explored in order to represent the mechanisms of action of tumor immunotherapies, as in [12,13,14,15,16,17,18,19]. Delay models of immunotherapy with chemotherapy have received attention, as in [20,21,22]. Rihan et al. [23] proposed a delay model incorporating optimal control variables to describe the dynamics of tumour-immune interactions in the presence of chemotherapy treatments. The specific goal of their research was to minimize the presence of tumor cells and to maximize the presence of effector cells. Based on the assumption that chemotherapy drugs kill cells at a different killing rate, it was suggested that treatment options are needed that will kill tumor cells and, at the same time, benefit immune cells.

Arcerio et al. [24] proposed a model describing tumor–effector cell interaction, as well as the pro-tumor role of the suppressive T-cells TG-$\beta$. They included siRNA treatment in their model to reduce the pro-tumor effect of TG-$\beta$ in the tumor microenvironment. Their findings indicated the occurrence of sustained oscillation which resulted in dormancy in the case of passive tumors, uncontrollable tumor growth for aggressive tumors, and controlled oscillatory tumor behavior when siRNA treatment was incorporated. It was suggested in [24] that effector cell contact with tumor cells in the tumor micro-environment reduces tumor concentration with interleukin-12 acting as a mediating variable. Recent findings have shown that such contact leads to immune cell modulation and subsequent production of cytokines by tumors to counterattack effector cells [25,26], and that interleukin-12 is secreted by effector T-cells [27].

Mathematical models have been used to gain insights into the hyperthermia treatment of malignancy. Rybinski et al. [28] adopted a heat-shock protein (HSP) model to estimate heat intensity and thermotolerance in the hyperthermia treatment of malignancy. Their paper highlighted the need for more extensive research based on experimental data using varying temperatures. Suleman et al. [29] reviewed the mathematical tools that have been used to gain insight into the hyperthermia treatment of malignacy and highlighted the need for more optimization research to improve treatment outcomes.

To obtain a more biologically accurate model of tumor–immune system interaction which can represent immune suppression and the uncontrollable growth of tumors, this investigation sought to modify the model in [24] by incorporating the ability of tumors to counterattack effector cells in the tumor micro-environment using tumor-derived cytokines, as shown recently in [25,26], and by incorporating hyperthermia treatment in the modified model to bring about tumor cell death and improvement of immune cell functioning. The remainder of this paper is structured as follows: In the next section, the model proposed in [24] is modified by incorporation of the ability of tumors to counterattack effector cells in the tumor micro-environment using tumor-derived cytokines, using a discrete time delay to represent the non-instantaneous modulation of immune cells. In addition, the future effect of hyperthermia treatment on the malignant tumor patient is investigated using the Lyapunov function. Section 3 presents a formulation of thermal optimal control, the nature and uniqueness of the thermal control, as well as its convexity and necessary conditions. A numerical simulation and its results are presented in Section 4. Section 5 concludes the paper with a summary of the findings.

2. The Model

To incorporate a tumors’ ability to counterattack effector cells using tumor-derived cytokines, we modified the model in [24] with the novel tumors’ production of cytokines to counterattack effector cells and therefore make the interaction dual-aggressive. Since the processes leading to immune cell modulation are not instantaneous, a discrete time delay is used to model the time it takes tumors to secrete inhibitory agents (${\tau }_1$) and develop suppressive cells (${\tau }_2$), as well as the time involved in the differentiation of effector cells to suppressive T-cells (${\tau }_3$). The purpose was to "factor in" the complexities involved in immune cell modulation into the dynamics of the interacting cells. The modified non-autonomous and non-dimensionalized version of the model in [24], as derived in Appendix A, is given below:

$$\begin{array}{cc}\hfill \frac{dE}{dt}& =\alpha E\left(t\right)+{\gamma }_{1}E\left(t\right)T\left(t\right)-{\upsilon }_{3}E\left(t\right)S\left(t\right)+{\omega }_{1}H\left(t\right)E\left(t\right)-\theta E\left(t\right),\hfill \\ \hfill \frac{dT}{dt}& =\delta T\left(t\right)+{\gamma }_{2}E\left(t\right)T(t-{\tau }_1)+{\upsilon }_{4}T\left(t\right)S\left(t\right)-{\omega }_{2}H\left(t\right)T\left(t\right)-\vartheta T\left(t\right),\hfill \\ \hfill \frac{dS}{dt}& ={\upsilon }_{1}E(t-{\tau }_2)+{\upsilon }_{2}T(t-{\tau }_3)-{\omega }_{3}H\left(t\right)S\left(t\right)-{\rm Y}S\left(t\right),\hfill \\ \hfill \frac{dH}{dt}& =\sigma H\left(t\right)-\Phi H\left(t\right).\hfill \end{array}$$

(1)

with initial conditions:

$$\begin{array}{c}\hfill E_0={\varphi }_1,T_0={\varphi }_2,S_0={\varphi }_3,H_0={\varphi }_{4}fort\in [-{\tau }_i,0].\end{array}$$

(2)

where E represents the effector cells, T represents the tumor cells, S represents the suppressive T-cells and H is the hyperthermia induction. The model parameters are described in

Table 1. System (1) Parameter Descriptions.

Parameter	Descriptions
$\alpha$	Activation rate of effector cells.
${\gamma }_1$	Surveillance rate of effector cells against tumor cells.
$\theta$	Apoptosis rate of effector cells.
$\delta$	Proliferation rate of tumor cells.
${\gamma }_2$	Inhibitory rate of tumor cells.
$\vartheta$	Death rate of tumor cells.
${\upsilon }_1$	Differentiation rate of effector cells to suppressive T-cells.
${\upsilon }_2$	Elaboration rate of suppression of cells by tumors.
${\upsilon }_3$	Suppression rate of effector cells by suppressive T-cells.
${\upsilon }_4$	Rate at which suppressive T-cells aid tumor escape.
${\rm Y}$	Death rate of suppressive T-cells.
${\omega }_1$	Rate at which heat boosts immune cell performance.
${\omega }_2$	Tumor shrinking rate due to heat induction.
${\omega }_3$	Rate at which heat reduces suppressive T-cells.
$\sigma$	Hyperthermia induction application rate.
$\Phi$	Control rate of hyperthermia-induction

as given below:

3. Qualitative Analysis

3.1. Existence of Non-Negative Solutions

The existence of the non-negativity of a solution of the system (1) is given in the lemma below.

Lemma 1.

Suppose $E\left(t\right)$, $T\left(t\right)$, $S\left(t\right)$ and $H\left(t\right)$ are positive and continuous under the initial conditions in Equation (2). Then, there exist non-negative solutions for $E\left(t\right),T\left(t\right),S\left(t\right),andH\left(t\right)$.

Proof of Lemma 1.

Solving for $E\left(t\right),T\left(t\right),S\left(t\right)andH\left(t\right)$ in system (1), we have

$$\begin{array}{cc}\hfill E\left(t\right)& ={\varphi }_1{e}^{{\int }_0^t(\alpha +{\gamma }_{1}T\left(t\right)-{\upsilon }_{3}S\left(t\right)+{\omega }_{1}H\left(t\right)-\theta )dt},\hfill \\ \hfill T\left(t\right)& ={\varphi }_2{e}^{{\int }_0^t(\vartheta +{\omega }_{2}H\left(t\right)-\delta -{\upsilon }_{4}S\left(t\right))dt}+{\int }_0^t\left({\gamma }_{2}T(t-{\tau }_1)E\left(t\right)\right)dt,\hfill \\ \hfill S\left(t\right)& ={\varphi }_3{e}^{({\rm Y}+{\omega }_{3}H)t}+{\int }_0^t({\upsilon }_{1}E(t-{\tau }_2)+{\upsilon }_{2}T(t-{\tau }_3))dt,\hfill \\ \hfill H\left(t\right)& =[{\varphi }_4\left(\rho \right)+M]e^{(\sigma -\Phi )t}.\hfill \end{array}$$

(3)

If M is a non-negative constant and assumptions in the lemma (1) hold, then $E\left(t\right),T\left(t\right),$$S\left(t\right)andH\left(t\right)$ are non-negative. □

3.2. Existence of Positive Steady States for Tumors’ Present

Theorem 1.

Assume $E^{*},T^{*},S^{*}{H}^{*}$ of system (1) are positive, then there exists a tumor-present steady-state at $E^{*}>0,T^{*}>0,S^{*}>0,H^{*}>0$.

Proof of Theorem 1.

If the assumption in Theorem (1) is true, then $E^{*},T^{*},S^{*},H^{*}$ of System (1) satisfy the following equations. Setting $\frac{dE}{dt}=0$, $\frac{dT}{dt}=0$, $\frac{dS}{dt}=0$ and $\frac{dH}{dt}=0$ in System (1) yield

$$\begin{array}{cc}\hfill \alpha E^{*}+{\gamma }_1{E}^{*}{T}^{*}-{\upsilon }_3{E}^{*}{S}^{*}+{\omega }_1{H}^{*}{E}^{*}-\theta E^{*}& =0\hfill \\ \hfill \delta T^{*}+{\gamma }_2{E}^{*}{T}^{*}+{\upsilon }_4{T}^{*}{S}^{*}-{\omega }_2{H}^{*}{T}^{*}-\vartheta T^{*}& =0\hfill \\ \hfill {\upsilon }_1{E}^{*}+{\upsilon }_2{T}^{*}-{\omega }_3{H}^{*}{S}^{*}-{\rm Y}{z}^{*}& =0\hfill \\ \hfill \sigma H^{*}-\Phi H^{*}& =0.\hfill \end{array}$$

(4)

Solving for $E^{*},T^{*}{S}^{*},H^{*}$ in Equation (4) gives

$$\begin{array}{cc}\hfill If& E^{*}\ne 0,then{T}^{*}=\frac{\theta +{\upsilon }_3{S}^{*}-\alpha -{\omega }_1{H}^{*}}{{\gamma }_1}\hfill \\ \hfill If& T^{*}\ne 0,then{E}^{*}=\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_4{S}^{*}}{{\gamma }_2}\hfill \\ \hfill If& S^{*}\ne 0,then{S}^{*}=n=\frac{{\gamma }_1{\upsilon }_1(\vartheta +{\omega }_2{H}^{*}-\delta )+{\gamma }_2{\upsilon }_2(\theta -\alpha -{\omega }_1{H}^{*})}{{\gamma }_1({\gamma }_2({\omega }_3{H}^{*}+{\rm Y})+{\upsilon }_1{\upsilon }_4)-{\gamma }_2{\upsilon }_2{\upsilon }_3}\hfill \\ \hfill If& H^{*}\ne 0,then\sigma =\Phi .\hfill \end{array}$$

(5)

Therefore, there exists a tumor-present steady-state

$$\begin{array}{c}\hfill (E^{*},T^{*},S^{*},H^{*})=\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{{\gamma }_2},\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{{\gamma }_1},n,H^{*}\right).\end{array}$$

(6)

□

3.3. Global Stability Analysis of Tumor-Present Steady-State

Biomedical interests are concerned with the likelihood of post-treatment conditions of a malignant tumor patient treated with hyperthermia. In this section, we seek to determine the conditions for the global stability of the hyperthermia treatment of tumors and the alleviation of suppressed immune cells using the Lyapunov function. The characteristic Equation of System (1) is thus;

$$\begin{array}{c}\hfill \left|\begin{array}{cccc}\lambda -(\alpha +{\gamma }_1{T}^{*}-\theta -{\upsilon }_3{S}^{*}+{\omega }_1{H}^{*})& {\gamma }_1{E}^{*}& -{\upsilon }_3{E}^{*}{e}^{-\lambda {\tau }_1}& {\omega }_1{E}^{*}\\ {\gamma }_2{T}^{*}& \lambda -(\delta -\vartheta +{\gamma }_2{E}^{*}{e}^{-\lambda {\tau }_2}+{\upsilon }_4{S}^{*}-{\omega }_2{H}^{*})& {\upsilon }_4{T}^{*}& -{\omega }_2{T}^{*}\\ {\upsilon }_1{E}^{*}{e}^{-\lambda {\tau }_3}& {\upsilon }_2{T}^{*}{e}^{-\lambda {\tau }_4}& \lambda +{\rm Y}+{\omega }_3{H}^{*}& {\omega }_3{S}^{*}\\ 0& 0& 0& \lambda -(\sigma -\Phi )\end{array}\right|=0\end{array}$$

(7)

Substituting $(E^{*},T^{*},S^{*},H^{*})=(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{\gamma },\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{\beta },n,H^{*})$ into Equation (7), we have

$$\begin{array}{cc}\hfill D(\lambda ,\tau )=& [\lambda -(\sigma -\Phi )][{\lambda }^3+(({\rm Y}+{\omega }_3{H}^{*})-(\delta -\vartheta +(\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n)\hfill \\ & e^{-\lambda {\tau }_1}+{\upsilon }_{4}n-{\omega }_2{H}^{*}\left)\right){\lambda }^2+((\vartheta -\delta -(\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n)e^{-\lambda {\tau }_1}-{\upsilon }_{4}n+{\omega }_2{H}^{*})\hfill \\ & ({\rm Y}+{\omega }_3{H}^{*})+{\upsilon }_2\left(\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{\beta }\right){\upsilon }_4\left(\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{\beta }\right)e^{-\lambda {\tau }_3})\lambda \hfill \\ & -\beta \left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{\gamma }\right)(\gamma \left(\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{\beta }\right)(\lambda +{\rm Y}+{\omega }_3{H}^{*})\hfill \\ & -{\upsilon }_1\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{\gamma }\right)e^{-\lambda {\tau }_2}{\upsilon }_4\left(\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{\beta }\right))\hfill \\ & -{\upsilon }_3\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{\gamma }(\gamma \left(\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{\beta }\right){\upsilon }_2\left(\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{\beta }\right)e^{-\lambda {\tau }_3}\hfill \\ & -{\upsilon }_1\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{\gamma }\right)e^{-\lambda {\tau }_2}(\lambda -(\delta -\vartheta +(\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n)e^{-\lambda {\tau }_1}\hfill \\ & +{\upsilon }_{4}n-{\omega }_2{H}^{*}\left)\right)\left)\right]=0\hfill \end{array}$$

(8)

Simplifying Equation (8), we have

$$\begin{array}{cc}\hfill D(\lambda ,\tau )=& {\lambda }^4+(A_{11}-A_{12}{e}^{-\lambda {\tau }_1}){\lambda }^3+(A_{21}-A_{22}{e}^{-\lambda {\tau }_1}+A_{23}{e}^{-\lambda {\tau }_2}+A_{24}{e}^{-\lambda {\tau }_3}){\lambda }^2\hfill \\ & -[A_{31}-A_{32}{e}^{-\lambda {\tau }_2}+(A_{33}{e}^{-\lambda {\tau }_3}+A_{34}{e}^{-\lambda ({\tau }_1+{\tau }_2)})]\lambda -(\sigma -\Phi )[{\lambda }^3\hfill \\ & +(A_{11}-A_{12}{e}^{-\lambda {\tau }_1}){\lambda }^2+(A_{21}-A_{22}{e}^{-\lambda {\tau }_1}+A_{23}{e}^{-\lambda {\tau }_2}+A_{24}{e}^{-\lambda {\tau }_3})\hfill \\ & \lambda -A_{31}+A_{32}{e}^{-\lambda {\tau }_2}-(A_{33}{e}^{-\lambda {\tau }_3}+A_{34}{e}^{-\lambda ({\tau }_1+{\tau }_2)})]=0,\hfill \end{array}$$

(9)

where,

$$\begin{array}{cc}\hfill A_{11}& =({\rm Y}+{\omega }_3{H}^{*})+(\vartheta -\delta -{\upsilon }_{4}n+{\omega }_2{H}^{*})\hfill \\ \hfill A_{12}& =\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n\hfill \\ \hfill A_{21}& =(\vartheta -\delta -{\upsilon }_{4}n+{\omega }_2{H}^{*})({\rm Y}+{\omega }_3{H}^{*})-(\vartheta -\delta -{\upsilon }_{4}n+{\omega }_2{H}^{*})(\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*})\hfill \\ \hfill A_{22}& =(\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n)({\rm Y}+{\omega }_3{H}^{*})\hfill \\ \hfill A_{23}& ={\upsilon }_1{\upsilon }_3{\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{{\gamma }_2}\right)}^2\hfill \\ \hfill A_{24}& ={\upsilon }_2{\upsilon }_4{\left(\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{{\gamma }_1}\right)}^2\hfill \\ \hfill A_{31}& =(\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n)(\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*})({\rm Y}+{\omega }_3{H}^{*})\hfill \\ \hfill A_{32}& ={\upsilon }_1{\upsilon }_3{\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{{\gamma }_2}\right)}^2(\vartheta -\delta -{\upsilon }_{4}n+{\omega }_2{H}^{*})\hfill \\ & +{\upsilon }_1{\upsilon }_4{\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{{\gamma }_2}\right)}^2(\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*})\hfill \\ \hfill A_{33}& ={\upsilon }_2{\upsilon }_3(\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n){\left(\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{{\gamma }_1}\right)}^2\hfill \\ \hfill A_{34}& ={\upsilon }_1{\upsilon }_3{\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{{\gamma }_2}\right)}^2(\vartheta -\delta -{\upsilon }_{4}n+{\omega }_2{H}^{*})\hfill \end{array}$$

(10)

and the scalar equation for the characteristic Equation (9) supposing $\lambda \equiv x$ is thus:

$$\begin{array}{cc}\hfill x^4=& -\left(A_{11}-A_{12}[t-{\tau }_1]\right)x^3-\left(A_{21}-A_{22}[t-{\tau }_1]+A_{23}[t-{\tau }_2]+A_{24}[t-{\tau }_3]\right)x^2\hfill \\ & +\left(A_{31}-A_{32}[t-{\tau }_2]+A_{33}[t-{\tau }_3]+A_{34}[t-[{\tau }_1+{\tau }_2]]\right)x\hfill \\ & +(\sigma -\Phi )[x^3+\left(A_{11}-A_{12}[t-{\tau }_1]\right)x^2+(A_{21}-A_{22}[t-{\tau }_1]+A_{23}[t-{\tau }_2]\hfill \\ & +A_{24}[t-{\tau }_3])x-[A_{31}-A_{32}[t-{\tau }_2]+A_{33}[t-{\tau }_3]+A_{34}[t-[{\tau }_1+{\tau }_2]]].\hfill \end{array}$$

(11)

The global stability condition for tumor-present steady-state with hyperthermia of incorporation $(E^{*},T^{*},S^{*},H^{*})$ of System (1) is ascertained using Lyapunov’s function. Firstly, definitions for the basic conditions for positive definite global for Lyapunov’s function are summarized below:

Definition 1.

Suppose $V(t,x_t)$ is a continuous scalar function $V:{\mathbb{R}}^{\mathbb{n}}-\mathbb{R}$

1.

If, $V(t,x_t)\to \infty$ as $\left|x\right|\to \infty$ then $V(t,x_t)$ is radially unbounded.

2.

If

$$\begin{array}{cc}\hfill V(t,x^{*})& =0,\hfill \\ \hfill V(t,x_t)& ⪈0,forx\ne x^{*}\hfill \end{array}$$

(12)

then $V(t,x_t)$ is positive definite globally.

Definition 2.

If the function $V(t,x_t)$ is radially unbounded and positive definite globally such that it has a globally negative time derivative, then

$$\begin{array}{c}\hfill V^{{}^{\prime }}(t,x_t)<0,\forall x_t\in {\mathbb{R}}^{*}\end{array}$$

(13)

and the invariant set is defined as follow;

$$\begin{array}{c}\hfill \mathfrak{J}=\{x_t\in {\mathbb{R}}^{*}:V^{{}^{\prime }}(t,x_t)=0\}.\end{array}$$

(14)

It suffices to say from the above that, if $\mathfrak{J}$ contains only $x^{*}$, then the steady state is globally stable. Therefore, the global stability of System (1) is defined by Theorem 2 below.

Theorem 2.

The tumor-present steady-state with hyperthermia incorporation for System (1) is globally asymptotically stable for ${\tau }_i$ where $i=1,2,3$, if, and only if,

$$\begin{array}{c}\hfill \sigma \le \Phi ,(\vartheta +{\omega }_2{H}^{*})\ge \delta +{\upsilon }_{4}n,({\rm Y}+{\omega }_3{H}^{*})\ge 0,\alpha +{\omega }_1{H}^{*}\ge \theta +{\upsilon }_{3}n.\end{array}$$

Proof of Theorem 2.

Suppose the Lyapunov function is defined as

$$\begin{array}{c}\hfill V(t,x_t)=\frac{x^2\left(t\right)}{2}\end{array}$$

(15)

Differentiating (15) along the solution of (11) gives

$$\begin{array}{cc}\hfill V^{{}^{\prime }}{x}^4=& -\left(A_{11}-A_{12}[t-{\tau }_1]\right)x^4-\left(A_{21}-A_{22}[t-{\tau }_1]+A_{23}[t-{\tau }_2]+A_{24}[t-{\tau }_3]\right)x^3\hfill \\ & +\left(A_{31}-A_{32}[t-{\tau }_2]+A_{33}[t-{\tau }_3]+A_{34}[t-[{\tau }_1+{\tau }_2]]\right)x^2\hfill \\ & +(\sigma -\Phi )[x^4+\left(A_{11}-A_{12}[t-{\tau }_1]\right)x^3+(A_{21}-A_{22}[t-{\tau }_1]+A_{23}[t-{\tau }_2]\hfill \\ & +A_{24}[t-{\tau }_3])x^2+[A_{31}-A_{32}[t-{\tau }_2]+A_{33}[t-{\tau }_3]+A_{34}[t-[{\tau }_1+{\tau }_2]]]x].\hfill \end{array}$$

(16)

Since $x^4(t-{\tau }_{})$ appears in (16), we need to generate a term like $-x^4(t-{\tau }_i)$ in (16) where $i=1,2,3$. Making use of Lemma 5.3 in [30], we assume

$$\begin{array}{c}\hfill V(t,x_t)=\frac{x^2\left(t\right)}{2}+\overline{\alpha }{\int }_{t-{\tau }_i}^t{x}^4\left(\varrho \right)d\varrho .\end{array}$$

where $\overline{\alpha }$ is a non-negative constant. Then Equation (16) becomes,

$$\begin{array}{cc}\hfill V^{{}^{\prime }}\left(x_t\right)=& -\left(A_{11}-A_{12}[t-{\tau }_1]-\overline{\alpha }\right)x^4-\left(A_{21}-A_{22}[t-{\tau }_1]+A_{23}[t-{\tau }_2]+A_{24}[t-{\tau }_3]\right)x^3\hfill \\ & +\left(A_{31}-A_{32}[t-{\tau }_2]+A_{33}[t-{\tau }_3]+A_{34}[t-[{\tau }_1+{\tau }_2]]\right)x^2\hfill \\ & +(\sigma -\Phi )[x^4+\left(A_{11}-A_{12}[t-{\tau }_1]\right)x^3+(A_{21}-A_{22}[t-{\tau }_1]+A_{23}[t-{\tau }_2]\hfill \\ & +A_{24}[t-{\tau }_3])x^2+[A_{31}-A_{32}[t-{\tau }_2]+A_{33}[t-{\tau }_3]+A_{34}[t-[{\tau }_1+{\tau }_2]]]x]\hfill \\ & -\overline{\alpha }{x}^4(t-{\tau }_i),\hfill \end{array}$$

(17)

where $i=1,2,3.$ It is obvious that (17) has all negative coefficients if $\sigma <\Phi$, and

$$\begin{array}{c}\hfill N\equiv A_{31}-A_{32}(t-{\tau }_2)+A_{33}(t-{\tau }_3)+A_{34}(t-[{\tau }_1+{\tau }_2])<0\end{array}$$

(18)

Computing (18) using definitions in (10) yields

$$\begin{array}{cc}\hfill N\equiv & (\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n)(\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*})({\rm Y}+{\omega }_3{H}^{*})\hfill \\ & -[{\upsilon }_1{\upsilon }_3{\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{{\gamma }_2}\right)}^2(\vartheta -\delta -{\upsilon }_{4}n+{\omega }_2{H}^{*})\hfill \\ & +{\upsilon }_1{\upsilon }_4{\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{{\gamma }_2}\right)}^2(\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*})](t-{\tau }_2)\hfill \\ & +{\upsilon }_2{\upsilon }_3(\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n){\left(\frac{\theta +{\upsilon }_{3}n-\alpha -{\omega }_1{H}^{*}}{{\gamma }_1}\right)}^2(t-{\tau }_3)\hfill \\ & +{\upsilon }_1{\upsilon }_3{\left(\frac{\vartheta +{\omega }_2{H}^{*}-\delta -{\upsilon }_{4}n}{{\gamma }_2}\right)}^2(\vartheta -\delta -{\upsilon }_{4}n+{\omega }_2{H}^{*})(t-[{\tau }_1+{\tau }_2]).\hfill \end{array}$$

(19)

Clearly, $N<0$ if

$$\begin{array}{c}\hfill \vartheta +{\omega }_2{H}^{*}>\delta +{\upsilon }_{4}n,{\rm Y}+{\omega }_3{H}^{*}>0,and\alpha +{\omega }_1{H}^{*}>\theta +{\upsilon }_{3}n.\end{array}$$

(20)

Therefore, $V^{{}^{\prime }}\left(x_t\right)<0$ if $\sigma <\Phi$, and

$$\begin{array}{c}\hfill \vartheta +{\omega }_2{H}^{*}>\delta +{\upsilon }_{4}n,{\rm Y}+{\omega }_3{H}^{*}>0,\alpha +{\omega }_1{H}^{*}>\theta +{\upsilon }_{3}n.\end{array}$$

(21)

Likewise, $V^{{}^{\prime }}\left(x_t\right)=0$, if, and only if, $\sigma =\Phi$, $\overline{\alpha }=0$ and

$$\begin{array}{cc}\hfill \vartheta +{\omega }_2{H}^{*}& =\delta +{\upsilon }_{4}n,\hfill \\ \hfill {\rm Y}+{\omega }_3{H}^{*}& =0,\hfill \\ \hfill \alpha +{\omega }_1{H}^{*}& =\theta +{\upsilon }_{3}n.\hfill \end{array}$$

(22)

The condition for having $V^{{}^{\prime }}\left(x_t\right)=0$ translates to having $E^{*}=T^{*}=S^{*}=H^{*}=0$. Hence, the largest invariant set with respect to the System (1) is

$$\begin{array}{c}\hfill \mathfrak{J}=\{(E^{*},T^{*},S^{*},H^{*})\in {\mathbb{R}}^{*}:V^{{}^{\prime }}(t,x_t)=0\}.\end{array}$$

(23)

Since (23) contains only the tumor-present steady-state $(E^{*},T^{*},S^{*},H^{*})$, the tumor-present steady-state with hyperthermia treatment is globally stable for ${\tau }_i$ where $i=1,2,3.$ □

In the next section, thermal optimal control is formulated for hyperthermia incorporation in the dynamics of tumor-immune cells, as described in (1). The theory of optimal control is employed to obtain thermal optimal control $\left(H\right)$ for efficient prevention of adverse effects of heat on normal cells.

3.4. Formulation of Thermal Optimal Control

In this section, thermal optimal control strategies for the hyperthermia treatment of malignant tumors are formulated to curb excessive high heating leading to adverse effects on normal cells. The need for heat optimization is advocated in the hyperthermia treatment of the malignant tumor, as excessively high heating will result in protein denaturation [3,8,9,10,11]. In the context of suggestions for heat optimization in the hyperthermia treatment of a malignant tumor, determination of the heat optimal control strategy for the prevention of protein denaturation which renders effector cells incapable of recognizing and binding to antigens was performed.

The next section provides a formulation of the objective function according to the multiple time-delays optimal control problem in [31]. The objective is to minimize only the control (hyperthermia dosage) for effective elimination of a tumor without causing an adverse effect on normal cells. Hence, the optimal control problem is formulated based on the related work in [32] as

$$\begin{array}{c}\hfill mimimize\mathbb{J}\left(H\right)=\frac{1}{2}{\int }_0^{t_f}{H}^{2}dt.\end{array}$$

bounded within $\sigma \le H\left(t\right)\le \Phi$ and subject to Equation (1) with initial condition values,

$$\begin{array}{cc}\hfill E\left(t_0\right)=0.001& ,T\left(t_0\right)=1,S\left(t_0\right)=0.0001.\hfill \end{array}$$

(24)

where $\sigma$ is the heat application rate and $\Phi$ is the heat control rate such that the heating induction is bounded between $\sigma$ and $\Phi$. To generate sufficient conditions for the existence of a finite objective functional value at the optimal control and states, the existence of solutions for thermal optimal controls are presented in Section 3.5.

3.5. Existence of Thermal Optimal Control Solution

The sufficient condition for the existence of a solution to thermal optimal control captured in Equation (24) is summarized by the Theorem 3.

Theorem 3.

If H is continuously differentiable in the constrain functions $(\frac{dE}{dt},\frac{dT}{dt},\frac{dS}{dt},\frac{dH}{dt})$ and convex on [$\sigma ,\Phi$] in objective functional $\mathbb{J}\left(\mathbb{H}\right)$ of (24), then there exists a unique thermal optimal solution $\overline{H}$ minimizing $\mathbb{J}\left(H\right)$ with $\mathbb{J}\left(\overline{H}\right)$ bounded.

Proof of Theorem 3.

To establish the existence and uniqueness of the solution of thermal optimal control as well as its convexity, it is required to state and verify if the following properties hold:

1.

The convexity condition for $\mathbb{J}\left(H\right)$ holds, that is $\frac{d^2\mathbb{F}}{dH}>0$.

2.

The thermal optimal solution of (24) is non-empty, that is $\mathbb{J}\left(\overline{H}\right)\le \mathbb{J}\left(H\right)$.

3.

The thermal optimal solution of (24) is unique.

4.

The thermal optimal solution of (24) is bounded.

Property 1: verifies if the thermal control H genuinely minimizes the objective functional $\mathbb{J}\left(H\right)$.

Defining the Hamiltonian $\mathbb{F}$ as

$$\begin{array}{cc}\hfill \mathbb{F}=& \frac{1}{2}{H}^2+{\lambda }_E(\alpha E\left(t\right)+{\gamma }_{1}E\left(t\right)T\left(t\right)-{\upsilon }_{3}E\left(t\right)S\left(t\right)+{\omega }_{1}H\left(t\right)E\left(t\right)-\theta E\left(t\right))\hfill \\ & +{\lambda }_T(\delta T\left(t\right)+{\gamma }_{2}E\left(t\right)T(t-{\tau }_1)+{\upsilon }_{4}T\left(t\right)S\left(t\right)-{\omega }_{2}H\left(t\right)T\left(t\right)-\vartheta T\left(t\right))\hfill \\ & +{\lambda }_S({\upsilon }_{1}E(t-{\tau }_2)+{\upsilon }_{2}T(t-{\tau }_3)-{\omega }_{3}H\left(t\right)S\left(t\right)-{\rm Y}S\left(t\right))+{\lambda }_H(\sigma -\Phi )H\left(t\right).\hfill \end{array}$$

(25)

Taking the second derivative of $\mathbb{F}$ with respect to H gives

$$\begin{array}{c}\hfill \frac{d^2\mathbb{F}}{d{H}^2}=1\end{array}$$

(26)

The convexity condition states that, if the second derivative of the twice differentiable Hamiltonian $\mathbb{F}$ with respect to the control H is non-negative then the optimal control problem (OPC) is minimization. Thus, the convexity condition holds since $\frac{d^2\mathbb{F}}{d{H}^2}>0$.Property 2: is verified using the definition in [32], then the tangent line property for the convexity condition is

$$\begin{array}{c}\hfill \frac{1}{2}{\overline{H}}^2-\frac{1}{2}{H}^2\le \frac{1}{2}(H{\overline{H}}^2-H^2){\mathbb{F}}_H\left(H^2\right).\end{array}$$

(27)

for every $\sigma \le t_0$, $t_1\le \Phi$. This yields

$$\begin{array}{cc}\hfill \mathbb{J}\left(\overline{H}\right)-\mathbb{J}\left(H\right)& =\frac{1}{2}{\int }_{t_0}^{t_1}({\overline{H}}^2-H^2)dt,\hfill \\ & \le \frac{1}{2}{\int }_{t_0}^{t_1}({\overline{H}}^2-H^2){\mathbb{F}}_H\left(H^2\right)dt.\hfill \end{array}$$

(28)

Using the optimality condition,

$$\begin{array}{c}\hfill {\mathbb{F}}_H{H}^2=-{\lambda }_E{\omega }_{1}E+{\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S-{\lambda }_H(\sigma -\Phi ).\end{array}$$

(29)

where ${\lambda }_H$ is the thermal control range, and $({\lambda }_E,{\lambda }_T,{\lambda }_S)$ are the thermal effects on effector cells, malignant tumors and suppressive T-cells, respectively. Substituting Equation (29) into Equation (28) gives,

$$\begin{array}{cc}\hfill \mathbb{J}\left(\overline{H}\right)-\mathbb{J}\left(H\right)& \le -\frac{1}{2}{\int }_{t_0}^{t_1}({\overline{H}}^2-u^2)({\lambda }_E{\omega }_1(\overline{E}-E)-{\lambda }_T{\omega }_2(\overline{T}-T)-{\lambda }_S{\omega }_3(\overline{S}-S)+{\lambda }_H(\sigma -\Phi ))dt.\hfill \end{array}$$

(30)

Taking into account that H is also convex in constraint functions $(\frac{dE}{dt},\frac{dT}{dt},\frac{dS}{dt},\frac{dH}{dt})$ suggest that

$$\begin{array}{cc}\hfill \mathbb{J}\left(\overline{H}\right)-\mathbb{J}\left(H\right)& \le 0.\hfill \end{array}$$

(31)

Thus, there exists a thermal optimal solution $\overline{H}$ for (24) since $\mathbb{J}\left(\overline{H}\right)\le \mathbb{J}\left(H\right)$.Property 3 (Uniqueness): this is verified by contradiction. Suppose that H is not a thermal optimal, such that there exists a control $\overline{H_1}$ on the interval $[\overline{t},t]$ with $\mathbb{J}\left(\overline{H_1}\right)\le \mathbb{J}\left(\overline{u}\right)$. Constructing a new thermal control $H_1$ on the interval $[t_0,t]$ yields

$$H_1\left(t\right)=\left\{\begin{array}{cc}\overline{H}\left(t\right)\hfill & for{t}_0\le t\le \overline{t}\hfill \\ \\ \overline{H_1}\left(t\right)\hfill & for\overline{t}\le t\le t_1\hfill \end{array}\right.$$

(32)

Hence,

$$\begin{array}{cc}\hfill \mathbb{J}{H}_1-\mathbb{J}\left(\overline{H}\right)=& \left({\int }_{t_0}^{\overline{t}}\frac{1}{2}\left(H_1^2\right)dt+\mathbf{J}{H}_1\right)-\left({\int }_{t_0}^{\overline{t}}\frac{1}{2}\left({\overline{H}}^2\right)dt+\mathbb{J}\left(\overline{H}\right)\right)\hfill \\ \hfill =& \mathbb{J}\left(\overline{H_1}\right)-\mathbb{J}\left(\overline{H}\right).\hfill \end{array}$$

(33)

Since $H_1$ exists on the extended interval, $\mathbb{J}{H}_1-\mathbb{J}\left(\overline{H}\right)>0$. This contradicts the condition for existence of the thermal optimal solution obtained. Thus, there exists no such control $H_1$ and $\overline{H}$ is the unique thermal control for (24).

Property 4 (Boundedness): The boundedness of $\overline{H}$ follows directly from the definition of optimal control problem (24). □

3.6. Characterization of Thermal Optimal Control

The characterizations for obtaining the optimality condition, the adjoint differential equation functions, and the transversality condition for thermal optimal control are summarized in Theorem 4 below.

Theorem 4.

(Necessary Condition). If $\overline{H}$ is a thermal optimal control for a bounded control H in (24) with associated states $(\overline{E},\overline{T},\overline{S})$ and $({\lambda }_E,{\lambda }_T,{\lambda }_S,{\lambda }_H)$ piece-wise differentiable functions, then the necessary conditions for $\overline{H}$, $(\overline{E},\overline{T},\overline{S})$ and $({\lambda }_E,{\lambda }_T,{\lambda }_S,{\lambda }_H)$ are:

1 

Optimality condition for H at $\sigma \le H\left(t\right)\le \Phi$ is

$$\overline{H}=\left\{\begin{array}{cc}\sigma ,\hfill & when\frac{\partial \mathbb{F}}{\partial \overline{H}}>0,\hfill \\ \sigma \le \overline{H}\le \Phi ,\hfill & when\frac{\partial \mathbb{F}}{\partial \overline{H}}=0,\hfill \\ \Phi ,\hfill & when\frac{\partial \mathbf{F}}{\partial \overline{H}}<0.\hfill \end{array}\right.$$

(34)

2 

Adjoint equations are:

$$\begin{array}{cc}\hfill {\lambda }_E^{{}^{\prime }}\left(t\right)& =-\frac{\partial \mathbb{F}}{\partial E},\hfill \\ \hfill {\lambda }_T^{{}^{\prime }}\left(t\right)& =-\frac{\partial \mathbb{F}}{\partial T},\hfill \\ \hfill {\lambda }_S^{{}^{\prime }}\left(t\right)& =-\frac{\partial \mathbb{F}}{\partial S},\hfill \\ \hfill {\lambda }_H^{{}^{\prime }}\left(t\right)& =-\frac{\partial \mathbb{F}}{\partial H}forall{t}_i\in t.\hfill \end{array}$$

(35)

Proof of Theorem 4.

For a bounded H within $\sigma \le H\left(t\right)\le \phi$, the optimality condition becomes

$$\begin{array}{c}\hfill \frac{\partial \mathbb{F}}{\partial H}=\overline{H}+{\lambda }_E{\omega }_{1}E-{\lambda }_T{\omega }_{2}T-{\lambda }_S{\omega }_{3}S+{\lambda }_H(\sigma -\Phi ).\end{array}$$

(36)

The sign of (36) is determined by the difference between $\overline{H}+{\lambda }_H(\sigma -\Phi )+{\lambda }_E{\omega }_{1}E$ and ${\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S$. Hence, the following are the possible cases:

1.

If $\overline{H}+{\lambda }_H(\sigma -\Phi )+{\lambda }_E{\omega }_{1}E>{\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S$, then

$$\begin{array}{c}\hfill \sigma =\frac{\partial \mathbb{F}}{\partial H}=\overline{H}+{\lambda }_H(\sigma -\Phi )+{\lambda }_E{\omega }_{1}E-{\lambda }_T{\omega }_{2}T-{\lambda }_S{\omega }_{3}S>0.\end{array}$$

(37)

2.

If $\overline{H}+{\lambda }_H(\sigma -\Phi )+{\lambda }_E{\omega }_{1}E={\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S$, then

$$\begin{array}{c}\hfill \sigma \le \overline{H}\le \Phi =\frac{\partial \mathbb{F}}{\partial H}=\overline{H}+{\lambda }_H(\sigma -\Phi )+{\lambda }_E{\omega }_{1}E-{\lambda }_T{\omega }_{2}T-{\lambda }_S{\omega }_{3}S=0.\end{array}$$

(38)

3.

If $\overline{H}+{\lambda }_H(\sigma -\Phi )+{\lambda }_E{\omega }_{1}E<{\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S$, then

$$\begin{array}{c}\hfill \Phi =\frac{\partial \mathbb{F}}{\partial H}=\overline{H}+{\lambda }_H(\sigma -\Phi )+{\lambda }_E{\omega }_{1}E-{\lambda }_T{\omega }_{2}T-{\lambda }_S{\omega }_{3}S<0.\end{array}$$

(39)

Case 1 describes the state where the thermal effect on effector cell rejuvenation and optimal heat availability is greater than the thermal effect on tumor and suppressive cells. This case depicts thermal excess for effector cell restoration.

Case 2 portrays the state where the thermal effect on tumor cells and suppressive cells is equal to the thermal effect on effector cell rejuvenation and the optimal heat available. This case depicts moderate thermal levels for both the elimination of tumor cells and the restoration of immune cells.

Case 3 interprets the state where the thermal effect on effector cell rejuvenation and the optimal thermal effect available is less than the thermal effect on tumor and suppressive cells. This case depicts thermal inadequacy for effector cells.

However, Case 2 depicts the required control strategies for effective elimination of a tumor and restoration of suppressed effector cells without causing an adverse effect. Thus, Case 2 matches the objective function, and will be adopted as the thermal optimal control in subsequent analysis. Therefore, from Equation (38)

$$\begin{array}{c}\hfill \overline{\mathbf{H}}=-{\lambda }_E{\omega }_{1}E+{\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S-{\lambda }_H(\sigma -\Phi ).\end{array}$$

(40)

The adjoint equations are:

$$\begin{array}{cc}\hfill {\lambda }_E^{{}^{\prime }}\left(t\right)& =-\frac{\partial \mathbf{F}}{\partial E}=-\left({\lambda }_E(\alpha +{\gamma }_{1}T-{\upsilon }_{3}S+{\omega }_{1}H-\theta )+{\lambda }_T\gamma T+{\lambda }_S{\upsilon }_1{e}^{-{\lambda }_E{\tau }_2}\right),\hfill \\ \hfill {\lambda }_T^{{}^{\prime }}\left(t\right)& =-\frac{\partial \mathbf{F}}{\partial T}=-\left({\lambda }_T(\delta +{\gamma }_{2}E{e}^{-{\lambda }_T{\tau }_1}+{\upsilon }_{4}S-{\omega }_{2}H-\vartheta )+{\lambda }_E\beta E+{\lambda }_S{\upsilon }_2{e}^{-{\lambda }_T{\tau }_3}\right),\hfill \\ \hfill {\lambda }_S^{{}^{\prime }}\left(t\right)& =-\frac{\partial \mathbf{F}}{\partial S}=-\left(-{\lambda }_E{\upsilon }_{3}E-{\lambda }_S({\omega }_{3}H+{\rm Y})+{\lambda }_T{\upsilon }_{4}T\right),\hfill \\ \hfill {\lambda }_H^{{}^{\prime }}\left(t\right)& =-\frac{\partial \mathbb{F}}{\partial H}=-\left(\overline{H}+{\lambda }_H(\sigma -\Phi )+{\lambda }_E{\omega }_{1}E-{\lambda }_T{\omega }_{2}T-{\lambda }_S{\omega }_{3}S\right)\hfill \end{array}$$

(41)

Solving (41) using integrating factors, we have

$$\begin{array}{cc}\hfill {\lambda }_E\left(t\right)& =c_1{exp}^{-(\alpha +{\gamma }_{1}T-{\upsilon }_{3}S+{\omega }_{1}H-\theta )t}-\left({\lambda }_T\gamma T+{\lambda }_S{\upsilon }_1{e}^{-{\lambda }_E{\tau }_2}\right)t\hfill \\ \hfill {\lambda }_T\left(t\right)& =c_2{exp}^{-(\delta T\left(t\right)+{\gamma }_{2}E\left(t\right)T(t-{\tau }_1)+{\upsilon }_{4}S-{\omega }_{2}H-\vartheta )t}-\left({\lambda }_E\beta E+{\lambda }_S{\upsilon }_2{e}^{-{\lambda }_T{\tau }_3}\right)t\hfill \\ \hfill {\lambda }_S\left(t\right)& =c_3{exp}^{-({\omega }_{3}H\left(t\right)S\left(t\right)+{\rm Y}S\left(t\right))t}+\left({\lambda }_E{\upsilon }_{3}E-{\lambda }_T{\upsilon }_{4}T\right)t\hfill \\ \hfill {\lambda }_S\left(t\right)& =c_3{exp}^{-({\omega }_{3}H\left(t\right)S\left(t\right)+{\rm Y}S\left(t\right))t}+\left({\lambda }_E{\upsilon }_{3}E-{\lambda }_T{\upsilon }_{4}T\right)t\hfill \\ \hfill {\lambda }_H\left(t\right)& =c_4{exp}^{-(\sigma -\phi )t}+\left({\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S-\overline{H}-{\lambda }_E{\omega }_{1}E\right)t\hfill \end{array}$$

(42)

For ${\lambda }_E\left(0\right)=c_1,{\lambda }_T\left(0\right)=c_{2,},{\lambda }_S\left(0\right)=c_3$ and ${\lambda }_H\left(0\right)=C_4$. For transversality conditions, we substitute the initial conditions $E\left(t_0\right)=0.001,T\left(t_0\right)=1,S\left(t_0\right)=0.0001$ into (42) and implement $\overline{H}:\sigma =\Phi .$ to obtain

$$\begin{array}{cc}\hfill {\lambda }_E\left(0.001\right)& =c_1{exp}^{-(\alpha E\left(t\right)+\beta T-{\upsilon }_{3}S+{\omega }_{1}H-\theta )0.001}-\left({\lambda }_T\gamma T+{\lambda }_S{\upsilon }_1{e}^{-{\lambda }_E{\tau }_2}\right)\left(0.001\right)=0,\hfill \\ \hfill {\lambda }_T\left(1\right)& =c_2{exp}^{-(\delta T\left(t\right)+\gamma E\left(t\right)T(t-{\tau }_1)+{\upsilon }_{4}S-{\omega }_{2}H-\vartheta )}-\left({\lambda }_E\beta E+{\lambda }_S{\upsilon }_2{e}^{-{\lambda }_T{\tau }_3}\right)=0,\hfill \\ \hfill {\lambda }_S\left(0.0001\right)& =c_3{exp}^{-({\omega }_{3}H\left(t\right)S\left(t\right)+{\rm Y}S\left(t\right))0.0001}+\left({\lambda }_E{\upsilon }_{3}E-{\lambda }_T{\upsilon }_{4}t\right)=0\hfill \\ \hfill {\lambda }_H\left(0\right)& =c_4+\left({\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S-\overline{H}-{\lambda }_E{\omega }_{1}E\right)t\hfill \end{array}$$

(43)

Solving for constants $c_1,c_2$ and $c_3$ in (43) gives

$$\begin{array}{cc}\hfill c_1& =\left({\lambda }_T\gamma T+{\lambda }_S{\upsilon }_1{e}^{-{\lambda }_E{\tau }_2}\right)\left(0.001\right){exp}^{(\alpha E\left(t\right)+\beta T-{\upsilon }_{3}S+{\omega }_{1}H-\theta )0.001},\hfill \\ \hfill c_2& =\left({\lambda }_E\beta E+{\lambda }_S{\upsilon }_2{e}^{-{\lambda }_T{\tau }_3}\right){exp}^{(\delta T\left(t\right)+\gamma ,E\left(t\right)T(t-{\tau }_1)+{\upsilon }_{4}S-{\omega }_{2}H-\vartheta )}\hfill \\ \hfill c_3& =-0.0001({\lambda }_E{\upsilon }_{3}E-{\lambda }_T{\upsilon }_{4}E){exp}^{({\omega }_{3}H\left(t\right)S\left(t\right)+{\rm Y}S\left(t\right))0.0001}.\hfill \\ \hfill c_4& ={\lambda }_H\left(0\right)-\left({\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S-\overline{H}-{\lambda }_E{\omega }_{1}E\right)t\hfill \end{array}$$

(44)

Substituting $c_1,c_2,c_3$ and $c_4$ into (42) yields

$$\begin{array}{cc}\hfill {\lambda }_E\left(t\right)& =(\left({\lambda }_T\gamma T+{\lambda }_S{\upsilon }_1{e}^{-{\lambda }_E{\tau }_2}\right)\left(0.001\right)[{exp}^{(\alpha E\left(t\right)+\beta T-{\upsilon }_{3}S+{\omega }_{1}H-\theta )[0.001-t]}-1]\hfill \\ \hfill {\lambda }_T\left(t\right)& =(\left({\lambda }_E\beta E+{\lambda }_S{\upsilon }_2{e}^{-{\lambda }_T{\tau }_3}\right)[{exp}^{-(\delta T\left(t\right)+{\gamma }_{2}E\left(t\right)T(t-{\tau }_1)+{\upsilon }_{4}S-{\omega }_2\mathbb{H}-\vartheta )[1-t]}-1]\hfill \\ \hfill {\lambda }_S\left(t\right)& =[1-0.0001{exp}^{({\omega }_3\mathbf{H}\left(t\right)S\left(t\right)+{\rm Y}S\left(t\right))[0.0001-t]}]\left({\lambda }_E{\upsilon }_{3}E-{\lambda }_T{\upsilon }_{4}T\right)t\hfill \\ \hfill {\lambda }_H\left(t\right)& ={\lambda }_H\left(0\right){exp}^{-(\sigma -\phi )t}+\left({\lambda }_T{\omega }_{2}T+{\lambda }_S{\omega }_{3}S-\overline{H}-{\lambda }_E{\omega }_{1}E\right)t[1-{exp}^{-(\sigma -\phi )t}]\hfill \end{array}$$

(45)

Hence, the thermal optimal control is

$$\begin{array}{c}\hfill \overline{H}:\sigma \le \overline{H}\le \Phi .\end{array}$$

(46)

To obtain optimal states $(\overline{E},\overline{T},\overline{S},\overline{H})$, inserting $\overline{H}$ in the constraint functions $(\frac{dE}{dt},\frac{dT}{dt},\frac{dS}{dt},\frac{dH}{dt})$ and solving for $(E,T,S,H)$ using integrating factors with initial conditions $E\left(t_0\right)=0.001,$ $T\left(t_0\right)=1,S\left(t_0\right)=0.0001$ and $H\left(t_0\right)$ is at ${37}^{\circ }$ body temperature. Then

$$\begin{array}{cc}\hfill E\left(t\right)& =0.001e^{{\int }_0^t(\alpha +\beta T\left(t\right)-{\upsilon }_{3}S\left(t\right)+{\omega }_1\overline{H}-\theta )dt},\hfill \\ \hfill T\left(t\right)& =e^{{\int }_0^t(\vartheta +{\omega }_2\overline{H}-\delta -{\upsilon }_{4}S\left(t\right))dt}+{\int }_0^t\left(\gamma T(t-{\tau }_1)E\left(t\right)\right)dt,\hfill \\ \hfill S\left(t\right)& =0.0001e^{({\rm Y}+{\omega }_3\overline{H})t}+{\int }_0^t({\upsilon }_{1}E(t-{\tau }_2)+{\upsilon }_{2}T(t-{\tau }_3))dt,\hfill \\ \hfill H\left(t\right)& =[{37}^{\circ }+M]e^{(\sigma \le \Phi )t}.\hfill \end{array}$$

where M is a constant. Thus, the thermal optimal control is

$$\begin{array}{c}\hfill \overline{H}:\sigma \le \overline{H}\le \Phi ,\end{array}$$

(47)

and the optimal states are,

$$\begin{array}{cc}\hfill E\left(t\right)& =0.001e^{{\int }_0^t(\alpha +\beta T\left(t\right)-{\upsilon }_{3}S\left(t\right)+{\omega }_1\overline{H}-\theta )dt},\hfill \\ \hfill T\left(t\right)& =e^{{\int }_0^t(\vartheta +{\omega }_2\overline{H}-\delta -{\upsilon }_{4}S\left(t\right))dt}+{\int }_0^t\left(\gamma T(t-{\tau }_1)E\left(t\right)\right)dt,\hfill \\ \hfill S\left(t\right)& =0.0001e^{({\rm Y}+{\omega }_3\overline{H})t}+{\int }_0^t({\upsilon }_{1}E(t-{\tau }_2)+{\upsilon }_{2}T(t-{\tau }_3))dt,\hfill \\ \hfill H\left(t\right)& =[{37}^{\circ }+M]e^{(\sigma \le \Phi )t}.\hfill \end{array}$$

□

4. Numerical Results

In this section, the numerical simulation of the global stability solution and the thermal optimal control strategies for tumor elimination and immune cell restoration for System (1) around the tumor-present steady-state with hyperthermia treatment are performed with a MATLAB $DDE23$ routine using the explicit Runge–Kutta scheme [33] with the parameter values in

Table 2. Model Parameter Values.

Parameter	Value	Unit	References
$\alpha$	0.1811	day${}^{-1}$	[34]
${\gamma }_1$	0.9	ng mL${}^{-1}$	[35,36]
$\theta$	1	day${}^{-1}$	[36,37]
$\delta$	0.1–1	day${}^{-1}$	[34,36]
${\gamma }_2$	3.5	ng mL${}^{-1}$	[38]
$\vartheta$	0.9	day${}^{-1}$	[35,38]
${\upsilon }_1$	1.8	day${}^{-1}$	[39,40,41]
${\upsilon }_2$	1.1	ng mL${}^{-1}$	[35,42]
${\upsilon }_3$	0.1–2.9	day${}^{-1}$	[36]
${\upsilon }_4$	0.27	day${}^{-1}$	[43]
${\rm Y}$	0.7	day${}^{-1}$	[39,44]
${\omega }_1$	0.362	day${}^{-1}$	[45]
${\omega }_2$	0.7	ng mL${}^{-1}$	[46]
${\omega }_3$	0.7	day${}^{-1}$	[46]

. This aims to validate the analytical results as well as to provide a graphical view of the cell dynamics during the implementation of different thermal control optimal strategies.

Figure 1 describes the time evolution of interacting cells for System (1) when the global stability conditions obtained analytically in Section 2 are simulated for ${\tau }_i=0$ and ${\tau }_i>0$ where $i=1,2,3$. (a) When ${\tau }_i=0$, all the interacting cell time evolution curves slide as time progresses to attain asymptotic stability. (b) When ${\tau }_i>0$, all the interacting cell time evolution curves rise higher than in the case (a) but slide aftermath as time progresses to attain asymptotic stability.

Figure 2 depicts thermal optimal control strategy in $\mathbf{Case}\mathbf{1}$ when the heat induction is in excess, that is when $\sigma =\overline{H}$. The heat concentration under this strategy increases as time progresses. However, the heat increase incites non-regulating growth of effector cells and both tumor cells and suppressive T-cells are seen to be overlapped and dormant in the early stage of treatment. The overlapped cells later split to grow independently at a later stage of treatment. The non-regulation of effector cells and adversities leading to uncontrollable tumor growth characterized by this case. This strategy created adverse effects and immune cell non-regulation occasioned by DNA damage due to excessive heating.

Figure 3 illustrates the optimal control strategy in $\mathbf{Case}\mathbf{3}$, when the heat induction is inadequate, that is when $\phi =\overline{H}$. The heat concentration under this strategy decreases as time progress (see Figure 3a). This situation produces tumor cells and suppressive T-cells relapse as the treatment commences. Effector cells are seen to be active and progress to attain asymptotic stability. However, the suppressive T-cell time evolution curve is seen to be higher than that of the effector cells. This strategy will eliminate tumor cells, rejuvenate effector cells, but account for a higher concentration of suppressive T-cells beyond the effector cells available.

Figure 4 illustrates the optimal control strategy in $\mathbf{Case}\mathbf{2}$ when the heat induction is moderate, that is when $\sigma \le \overline{H}\le \phi$. The heat concentration under this strategy remains steady and unchanged as time progress (see Figure 4a). This situation creates complete tumor elimination as the tumor regresses to zero. Effector cells rejuvenate and subsequently attain an asymptotic stable state. The suppressive T-cells are seen ascending without inhibition or alignment with effector cells (see Figure 4b). Thus, moderate thermal optimal control guarantees total tumor elimination, rejuvenation of suppressed effector cells, and control of anti-effector cell effects on suppressive T-cells without causing adverse effects.

5. Conclusions

Here, we investigated the aftermath of hyperthermia treatment of a malignant tumor patient, and the adverse effects caused by excessive heating leading to protein denaturation and DNA damage. This was achieved by investigation of the global stability status of the tumor-present steady-state with hyperthermia treatment, and by formulation of a thermal optimal control problem to minimize heat-induction so as to reduce the adverse effect of heat on normal cells.

The global stability analysis was performed using Lyapunov’s function. This revealed the conditions for global stability of the tumor-present steady-state with hyperthermia incorporation. The global stability conditions suggest that the tumor-present steady-state with hyperthermia treatment is globally stable if hyperthermia treatment leads to the following:

1.

The rate of tumor death is greater than or equal to the rate of its progression due to self-proliferation and influence of suppressive T-cells.

2.

The rate of suppressive T-cell apoptosis remains zero.

3.

The rate at which effector cells proliferate is greater than or equal to their rate of apoptosis and regulation by T-regs.

The thermal optimal control problem to minimize heat induction was initially studied to determine the existence and uniqueness of thermal optimal control for hyperthermia treatment. The convexity of the problem is well-established. The characterization of the thermal optimal control problem enhances the optimality condition, and the thermal optimal control impacts on the interactive cells in the course of treatment (adjoint ${\lambda }_{E,T,S}$). The optimality conditions indicate optimal control strategies which may lead to the following:

1.

Thermal optimal control leading to excessive heating that generates adversity and immune cell non-regulation occasioned by DNA damage.

2.

Thermal optimal control leading to inadequate heating that would eliminate tumors, rejuvenate effector cells, but account for a higher concentration of suppressive T-cells beyond the effector cells available.

3.

Thermal optimal control leading to moderate heating that would lead to the elimination of tumors and the restoration of effector cells without adversity.

The control leading to both the elimination of tumor cells and the restoration of effector cells without adversity was used to obtain the adjoint equations and optimal solutions for Equation (24).

The numerical simulation results validated the analytical results and suggested the cell dynamics when implementing thermal optimal control strategies. When there is inadequate heating, it was found that the heat concentration keeps dropping; elimination of tumor cells and rejuvenation of effector cells is also manifest in this case. However, a higher concentration of suppressive T-cells beyond the available effector cells is seen (see Figure 3). In the case of excessive heating during the hyperthermia treatment of malignant tumors shown in Figure 2, the heat concentrations keep rising as the treatment continues. This strategy leads to tumor elimination but causes immune cell non-regulation occasioned by excessive heat-induced DNA damage. For moderate heating, the heat concentration remains steady as treatment progresses. Tumor elimination and the rejuvenation of immune cells without causing adversity are also characterized by this case (see Figure 4).

This investigation addresses concerns regarding the aftermath of hyperthermia treatment on malignant tumor patients, and the adverse outcomes arising from excessive heating leading to protein denaturation and DNA damage. It is anticipated that the results obtained will inform the application of hyperthermia in the treatment of malignant tumors. They also confirm the efficacy of hyperthermia treatment and provide insight into possible adverse effects in the case of inadequate or excessive heat application during hyperthermia treatment. Hyperthermia treatment is a complex process, so the mathematical study of the thermal chemistry leading to effector cell performance improvement, especially regarding the role of heat-shock proteins (HSPs) in the restoration of suppressed cells, should benefit and enrich thermotherapy.