Optimal Control Theory and Calculus of Variations in Mathematical Models of Chemotherapy of Malignant Tumors

Abstract

This paper is devoted to the analysis of mathematical models of chemotherapy for malignant tumors growing according to the Gompertz law or the generalized logistic law. The influence of the therapeutic agent on the tumor dynamics is determined by a therapy function depending on the time-varying concentration of the drug in the patient’s body. The case of a non-monotonic therapy function with two maxima is studied. It reflects the use of two different therapeutic agents. The state variables of the dynamics are the tumor volume and the amount of the therapeutic agent able to suppress malignant cells (concentration of the drug in the body). The treatment protocol (the rate of administration of the therapeutic agent) is the control in the dynamics. The optimal control problem for this models is considered. It is the problem of the construction of treatment protocols that provide the minimal tumor volume at the end of the treatment. The solution of this problem was obtained by the authors in previous works via the optimal control theory. The form of the considered therapy functions provides a specific structure for the optimal controls. The managerial insights of this structure are discussed. In this paper, the structure of the viability set is described for the model according to the generalized logistic law. It is the set of the initial states of the model for which one can find a treatment protocol that guarantees that the tumor volume remains within the prescribed limits throughout the treatment. The description of the viability set’s structure is based on the optimal control theory and the theory of Hamilton–Jacobi equations. An inverse problem of therapy is also considered, namely the problem of reconstruction of the treatment protocol and identification of the unknown parameter of the intensity of the tumor growth. Reconstruction is carried out by processing information about the observations of the tumor volume dynamics and the measurements of the drug concentration in the body. A solution to this problem is obtained through the use of a method based on the calculus of variations. The results of the numerical simulations are presented herein.

1. Introduction

This paper is devoted to the study of mathematical models of chemotherapy of malignant tumors with different laws of the growth of the tumor volume (the generalized logistic law and Gompertz’s law). The research focuses on solid tumors, such as breast, lung, liver, colon, pancreas and prostate cancers. As known, a malignant tumor formed in an organ is capable of uncontrolled or poorly controlled rapid growth of its cells and can penetrate into adjacent tissues, damage them and penetrate further into other organs, resulting in metastases. Spatially homogeneous mathematical models of the growth of solid avascular tumors are considered (on this theme, see, for example, [1,2]).

Usually, in such mathematical models, the total volume of the tumor is considered. The tumor’s dynamics are described using ordinary or partial differential equations. In this paper, we consider tumors with dynamics described by ordinary differential equations that obey the generalized logistic law or the Gompertz law [2,3,4].

Chemotherapy is one of the most effective ways to fight cancer. The effect of a therapeutic agent on a tumor is defined by adding a therapy function to equations of the tumor’s dynamics. This function depends on the time-varying concentration of the drug in the body. The case of a non-monotonic therapy function with two maxima is studied. It reflects the use of two different therapeutic agents. This case has not been previously considered. The infusion and the dissipation of the drug in the body are described in the considered model by ordinary differential equations.

Let us also note another well-known mathematical model of cancer therapy: the Lotka–Volterra predator–prey model. In this model, the therapeutic agent plays the role of the predator, and the tumor plays the role of the victim. For further information on this model, see, for example, [5,6,7].

The tumor growth dynamics are considered over a certain fixed period of time. At the final time instant (the control point), the total tumor volume is estimated. The state variables of the dynamics are the tumor volume and the amount of the therapeutic agent able the to suppress malignant cells (concentration of the drug in the body). The treatment protocol (the rate of administration of the therapeutic agent) is the control in the dynamics.

The optimal control problem is considered. The problem is to construct a treatment protocol (a control) for any initial state of the model (the tumor volume and the concentration of the drug in the body at the beginning of the treatment) such that the tumor volume is minimal at the control point. An optimal treatment strategy and optimal protocols were suggested in [8,9]. Their construction is based on the results of the optimal control theory [10] and the theory of generalized solutions of partial differential equations of the first order [11,12]. The optimal treatment strategy is a feedback depending on the current state of the model: the tumor volume and the concentration of the drug in the body. The optimal treatment protocol is formed on the basis of information about the current state of the model according to the optimal strategy. For any fixed initial state of the model, the optimal treatment protocol ensures the minimal tumor volume at the end of the treatment. It is has been proven that the optimal protocol has a piecewise constant structure with no more than one switch. Construction of optimal therapy strategies and optimal treatment protocols has been previously discussed, for example, in [3,5,6,7,13,14,15,16,17,18]. The innovation of the research presented in this paper is that a new type of therapy function is considered, namely a non-monotonic therapy function with two maxima. It reflects the use of two different therapeutic agents. Such a therapy function provides a specific structure of the optimal control. The managerial insights of such an optimal control are discussed in this paper.

The viability set [19] is described for the model according to the generalized logistic law. It is the set of initial states of the model for which one can find a treatment protocol that guarantees that the tumor volume remains within the prescribed limits and is compatible with life throughout the treatment. In this paper, we propose and substantiate the description of the structure of the viability set for the model according to the generalized logistic law. The justifications of the constructions are based on the theory of generalized solutions of partial differential equations of the first order [11,12]. The viability set for the dynamics according to the Gompertz law was described in [20].

Research on models of chemotherapy of malignant tumors is also of interest to solve inverse problems, for example, identification of the parameter of the model characterizing the intensity of the tumor’s growth and reconstruction of the treatment protocol (the control) for a given patient using the observations of the tumor volume dynamics and the current measurements of the concentration of the drug in the body.

The results of the numerical simulations of the solution of such inverse problems are presented in this paper. They were obtained using a numerical method [21] based on the theory of calculus of variations [22].

It should be noted that a process as complex as the growth and suppression of malignant cells is presented in the models in a simplified, idealized form. Still, the obtained results may be of use for the analysis of experimental data or for the study of new effective therapy methods.

2. Materials and Methods

2.1. Dynamics of the Model

The state variables of the dynamics are the tumor volume ($m\left(t\right)$) and the amount of the chemotherapy drug ($h\left(t\right)$) that is able to suppress malignant cells. The therapy function ($f\left(h\right)$) describes the effect of the drug on the malignant cells. The control (the treatment protocol) ($u\left(t\right)$) is the amount of the drug injected within the time unit.

The dynamics of the malignant tumor and the concentration of the drug are described by the known model [2,3,4]:

$$\begin{array}{c}{\displaystyle \frac{dm\left(t\right)}{dt}=g\left(m\left(t\right)\right)-\gamma m\left(t\right)f\left(h\left(t\right)\right),}\\ {\displaystyle \frac{dh\left(t\right)}{dt}=-\alpha h\left(t\right)+u\left(t\right),}\\ m\left(t_0\right)=m_0,h\left(t_0\right)=h_0,t\in [t_0,T].\end{array}$$

(1)

where $0\le t_0<T<\infty$ are the starting and control points, $\gamma >0$ is the effectiveness coefficient of the therapy and $\alpha >0$ is the dissipation coefficient. The feasible controls are measurable (piecewise constant) bounded functions:

$$u(·):[t_0,T]\to \mathbf{U}\triangleq [0,Q],$$

(2)

where Q is the maximum allowed value.

In the model (1), the function $g\left(m\right)$ describes the law of the tumor’s growth. In this paper, two different models of this law are considered.

• The Gompertz law:

  $$g_1\left(m\right)=rm-\theta mln\left(m\right),$$

  (3)

  where the parameter $r>0$ characterizes the speed of the tumor’s growth, and $\theta >0$ characterizes the speed of tumor suppression.
• The generalized logistic law:

  $$g_2\left(m\right)=rm\left(1-{\left[\frac{m}{\theta }\right]}^{\beta }\right),$$

  (4)

  where the parameter $r>0$ again characterizes the speed of the tumor’s growth, $\theta >0$ characterizes the threshold tumor volume and $0<\beta \le 1$ changes the steepness of the model’s curve.

These models are described in [2,3,4].

In the both models, the following restrictions are considered:

$$0\le m_0\le M,0\le h_0\le L,$$

where M is the threshold tumor volume compatible with life, and L is the maximum allowed amount of the drug in the body.

2.2. The Therapy Function

The therapy function ($f\left(h\right)$) describes the effect of the drug on the malignant cells and is an important part of the considered model (1). In this paper, we consider a non-monotonic, continuously differentiable, positive therapy function such that its derivative ($f^{\prime }\left(h\right)$) has three different roots:

$$0<{\widehat{h}}_1<{\widehat{h}}_2<{\widehat{h}}_3\le L,f^{\prime }\left({\widehat{h}}_{1,2,3}\right)=0.$$

We assume that $f\left(h\right)$ satisfies the following conditions:

A1.

$f^{\prime }\left(h\right)>0\forall h<{\widehat{h}}_1,f^{\prime }\left(h\right)<0\forall h>{\widehat{h}}_3$.

A2.

$0<\alpha {\widehat{h}}_{1,2,3}<Q$.

A3.

$f\left({\widehat{h}}_1\right)=f\left({\widehat{h}}_3\right)$.

A4.

$f^{\prime }\left(h\right)<0\forall h\in ({\widehat{h}}_1,{\widehat{h}}_2),f^{\prime }\left(h\right)>0\forall h\in ({\widehat{h}}_2,{\widehat{h}}_3)$.

It follows from conditions $A1$–$A4$ that the roots ${\widehat{h}}_1$ and ${\widehat{h}}_3$ are the maximum points and that root ${\widehat{h}}_2$ is the minimum point of the therapy function ($f\left(h\right)$).

2.3. Optimal Control Problem

The optimal control problem consists of constructing feasible (piecewise constant) feedback that minimizes the terminal therapy quality index (the value function) ($\sigma \left(m\right(T\left)\right)$). It is different for the two considered models.

• For the Gompertz law (3):

  $${\sigma }_1\left(m_1\left(T\right)\right)={m_1}^2(T;t_0,h_0,m_0,u(·))\to \underset{u(·)}{min}.$$

  (5)
• For the generalized logistic law (4):

  $${\sigma }_2\left(m_2\left(T\right)\right)={m_2}^{\beta }(T;t_0,h_0,m_0,u(·))\to \underset{u(·)}{min}.$$

  (6)

  where $m_{1,2}\left(t\right)=m_{1,2}(t;t_0,h_0,m_0,u(·))$ are the solutions of the system (1), that are generated by a feasible control ($u\left(t\right)$) for $g\left(m\right)=g_1\left(m\right)$ (3) and $g\left(m\right)=g_2\left(m\right)$ (4), respectively, for the boundary conditions $(t_0,h_0,m_0)$.

The system (1) can be integrated analytically for both $g_1$ and $g_2$. The solutions are:

• For the Gompertz law (3):

  $${\displaystyle m_1\left(t\right)=m_0^{-\theta (t-t_0)}exp\left[\underset{t_0}{\overset{t}{\int }}{e}^{-\theta (t-\tau )}\left(r-\gamma f\left(h\right(\tau \left)\right)\right)d\tau \right].}$$
• For the generalized logistic law (4):

  $${\displaystyle m_2\left(t\right)=\frac{{\displaystyle \theta m_{0}exp\left[\underset{t_0}{\overset{t}{\int }}[r-\gamma f\left(h\left(\tau \right)\right)]d\tau \right]}}{{\displaystyle {\left({\theta }^{\beta }+\beta m_0^{\beta }r\underset{t_0}{\overset{t}{\int }}exp\left[\beta \underset{\tau }{\overset{T}{\int }}[r-\gamma f\left(h\left(y\right)\right)]dy\right]d\tau \right)}^{\frac{1}{\beta }}}}.}$$

  (7)

In each case, $h\left(t\right)=h(t;t_0,h_0,u(·))$ is the solution of the second independent equation of (1).

2.4. Value Function and Optimal Synthesis

We introduce the value function in problems (1), (5) and (1), (6). It matches each boundary condition ($(t_0,h_0,m_0)\in [0,T]\times [0,L]\times [0,M]$) with the optimal values $\mathrm{V}\mathrm{a}{\mathrm{l}}_{1,2}$ respectively for (5) and (6):

$$\mathrm{V}\mathrm{a}{\mathrm{l}}_{1,2}(t_0,h_0,m_0)=\underset{u(·)}{min}{\sigma }_{1,2}\left(m_{1,2}\left(T\right)\right)$$

According to [8,9],:

• For the Gompertz law (3) [8]:

  $${\displaystyle \mathrm{V}\mathrm{a}{\mathrm{l}}_1(t_0,h_0,m_0)=m_0^2{e}^{-\theta (T-t_0)}exp\left[2\frac{r}{\theta }(1-e^{-\theta (T-t_0)})\right]exp\left[-2\gamma V_1(t_0,h_0)\right];}$$

  (8)
• For the generalized logistic law (4) [8]:

  $${\displaystyle \mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)=\frac{{\theta }^{\beta }{m}_0^{\beta }{e}^{-\beta \gamma V_2(t_0,h_0)}}{{\theta }^{\beta }{e}^{-\beta \left(r(T-t_0)\right)}+\beta m_0^{\beta }r\underset{t_0}{\overset{T}{\int }}exp\left[-\beta (r(\tau -t_0)+\gamma V_2(\tau ,h^0\left(\tau \right)))\right]d\tau }.}$$

  (9)

Here,

$$V_{1,2}(t_0,h_0)=\underset{u(·)}{max}{J}_{1,2}(t_0,h_0,u(·))$$

(10)

is the optimal value in the reduced optimal control problem:

$${\displaystyle \frac{dh\left(t\right)}{dt}=-\alpha h\left(t\right)+u\left(t\right),h\left(t_0\right)=h_0.}$$

For the Gompertz law (3), the cost value functional is

$${\displaystyle J_1(t_0,h_0,u(·))=\underset{t_0}{\overset{T}{\int }}{e}^{-\theta (T-\tau )}f\left(h(t;t_0,h_0,u(·))\right)dt\to \underset{u(·)}{max}.}$$

(11)

For the generalized logistic law (4),

$${\displaystyle J_2(t_0,h_0,u(·))=\underset{t_0}{\overset{T}{\int }}f\left(h(t;t_0,h_0,u(·))\right)dt\to \underset{u(·)}{max}.}$$

(12)

The value functions (10) are constructed piecewise from several functions (${\phi }_i(·)$):

$$\begin{array}{c}{V}_{1,2}(t,h)=\left\{\begin{array}{c}{\phi }_1(t,h),(t,h)\in G_1,\hfill \\ {\phi }_2(t,h),(t,h)\in G_2,\hfill \\ {\phi }_3(t,h),(t,h)\in {\mathrm{\Pi }}_1,\hfill \\ {\phi }_4(t,h),(t,h)\in {\mathrm{\Pi }}_2,\hfill \\ {\phi }_5(t,h),(t,h)\in {\mathrm{\Pi }}_3.\hfill \end{array}\right.\end{array}$$

(13)

The geometric regions are

$$\begin{array}{c}{G}_1=\{(t,h):t\in [0,T],h={\widehat{h}}_1\},G_2=\{(t,h):t\in [0,T],h={\widehat{h}}_3\},\\ {\mathrm{\Pi }}_1=[0,T]\times [0,{\widehat{h}}_1),{\mathrm{\Pi }}_2=[0,T]\times ({\widehat{h}}_3,L],\\ {\mathrm{\Pi }}_3=[0,T]\times ({\widehat{h}}_1,x\left(t\right)],{\mathrm{\Pi }}_4=[0,T]\times (x\left(t\right),{\widehat{h}}_3).\end{array}$$

(14)

For visualization, see Figure 1. The graph of the $x(·)$ function is the Rankine–Hugoniot line

$$\mathrm{\Gamma }=\{(t,x\left(t\right)):t\in [0,T],x\left(T\right)={\widehat{h}}_2\},$$

as defined by the points where

$${\phi }_5(t,x\left(t\right))={\phi }_6(t,x\left(t\right)).$$

It is shown in [9] that the $\mathrm{\Gamma }$ line satisfies the Rankine–Hugoniot conditions with a boundary condition of $x\left(T\right)={\widehat{h}}_2$ and has the following form:

$$\frac{dx}{dt}=-\alpha x\left(t\right)+Q\frac{s_h^6\left(t\right)}{s_h^6\left(t\right)-s_h^5\left(t\right)},$$

(15)

where

$${\displaystyle s_h^5\left(t\right)=\frac{\partial {\phi }_5(t,x\left(t\right))}{\partial h},s_h^6\left(t\right)=\frac{\partial {\phi }_6(t,x\left(t\right))}{\partial h}.}$$

The functions ${\phi }_i(t,h),i=\overline{1,6}$ are constructed by the Cauchy characteristic method for the auxiliary linear Hamilton–Jacobi equations with boundary conditions of the special form [9].

It was proven in [8,20] that the optimal strategy in problems (1), (5) and (1), (6) has the following form:

$$u^0(t,h)=\left\{\begin{array}{cc}\alpha {\widehat{h}}_1,\hfill & (t,h)\in G_1,\hfill \\ \alpha {\widehat{h}}_3,\hfill & (t,h)\in G_2,\hfill \\ Q,\hfill & (t,h)\in {\mathrm{\Pi }}_1,\hfill \\ 0,\hfill & (t,h)\in {\mathrm{\Pi }}_2,\hfill \\ 0,\hfill & (t,h)\in {\mathrm{\Pi }}_3,\hfill \\ Q,\hfill & (t,h)\in {\mathrm{\Pi }}_4\setminus \mathrm{\Gamma }.\hfill \end{array}\right.$$

(16)

Note that all optimal controls ($u^0\left(t\right)=u^0(t,h^0\left(t\right))$) satisfy the Pontryagin maximum principle.

Remark 1. 

The process of construction of the optimal treatment strategy (16) is explained and justified in [8,20]. The optimal strategy is a feedback, so the optimal protocols are defined by the current state of the system $\left(m\right(t),h(t\left)\right)$. Their values are chosen following (16).

When the state of the system is in each of the regions (${\mathrm{\Pi }}_1,{\mathrm{\Pi }}_2,{\mathrm{\Pi }}_3,{\mathrm{\Pi }}_4$), the control is either minimal (0) or maximal (Q). When the state of the model transfers from one region to another, the control is switched. The specifics are that the border between regions ${\mathrm{\Pi }}_3$ and ${\mathrm{\Pi }}_4$ is defined by the Rankine-Hugoniot line (Γ). The conditions that determine this line are non-trivial (15).

The other feature is that the optimal control does not always have a maximal or minimal value. If the state of the system belongs to lines $G_1$ and $G_2$, the optimal control corresponds to special regimes $u^0(t,h)=\alpha {\widehat{h}}_1$ and $u^0(t,h)=\alpha {\widehat{h}}_3$. These regimes are stable and provide the optimal effect of the drug, which is achieved at the maximum points of the therapy function.

It follows from Pontryagin’s maximum principle, which provides the basis for the optimal strategy construction, that this control provides the optimal terminal therapy quality index (5) and (6). In other words, it minimizes the tumor volume ($m\left(T\right)$) at the control point.

The Rankine-Hugoniot line (Γ) plays an important role in the construction of the optimal protocols. The selection of the therapeutic agent depends on what regions the initial state of the model belongs to: ${\mathrm{\Pi }}_2,{\mathrm{\Pi }}_4,G_2$ or ${\mathrm{\Pi }}_1,{\mathrm{\Pi }}_3,G_1$.

2.5. The Viability Set

We now define the viability set for the generalized logistic law (1), (6).

Let $m_2^0\left(t\right)=m(t;t_0,h_0,m_0,u^0(·)))$ be the solution of system (1), (6) with boundary conditions $(t_0,h_0,m_0)$ that is generated by the optimal control ($u^0\left(t\right)$).

Definition 1. 

The viability set for problem (1), (6) according to the generalized logistic law is a set of the points ($(t_0,h_0,m_0)\in [0,T]\times [0,L]\times [0,M_2]$) such that

$$m_2^0(t;t_0,h_0,m_0,u^0(·))\le M_2,t\in [t_0,T],$$

where $M_2$ is the threshold tumor volume compatible with life specified for the model according to the generalized logistic law.

The parameter $M_2$ is the stable equilibrium point of the dynamics (1) when the drug effectiveness (that is, the therapy function ($f\left(h\right)$)) is maximal [3]. Thus, consider the following equation:

$${\displaystyle \frac{dm}{dt}=rm\left(1-{\left[\frac{m}{\theta }\right]}^{\beta }\right)-\gamma mF,}$$

(17)

where

$$F=f\left({\widehat{h}}_1\right)=f\left({\widehat{h}}_3\right)=\underset{h\in [0,L]}{max}f\left(h\right).$$

Let us find the roots of the right side of (17). The trivial zero root provides an unstable equilibrium. Stable equilibrium takes place when $m=\tilde{m}$,

$$\tilde{m}=\theta {\left[1-\frac{\gamma F}{r}\right]}^{\frac{1}{\beta }}.$$

(18)

Therefore, we assume

$$M_2=\tilde{m}.$$

We introduce set $W_2$:

$$W_2=\left\{(t_0,h_0,m_0)\in [0,T]\times [0,L]\times [0,M_2]:\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)\in [0,{M_2}^{\beta }]\right\}.$$

(19)

We will show that this set is the maximal viability set for problem (1), (6).

Definition 2. 

Set W is the maximal viability set for problem (1), (6) if, for any point ($w=(t_0,h_0,m_0)\notin W$) and any measurable function ($u(·):[t_0,T]\to [0,Q]$), there exists $t^{*}\in (t_0,T)$ such that

$$m(t^{*};t_0,h_0,m_0,u(·))>M_2.$$

3. Results

3.1. Construction of Set $W_2$

In this section, we describe the structure of set $W_2$.

First, we transform the expression (9) for $\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)$ into a simpler form. It follows from (12) that in expression (9),

$${\displaystyle V_2(t_0,h_0)=\underset{u(·)}{max}{J^2}_{t_0,h_0}\left(u(·)\right)=\underset{t_0}{\overset{T}{\int }}f\left(h^0\left(t\right)\right)dt,}$$

where $h^0\left(t\right)=h(t;t_0,h_0,u^0(·))$ is the solution of the second independent equation of the system (1) that is generated by the optimal feedback control ($u^0(t,h)$ (16)). Also note that $\mathrm{V}\mathrm{a}{\mathrm{l}}_2(T,h_0)=0$. Then,

$$\begin{array}{c}{\displaystyle \mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)=\frac{{\theta }^{\beta }{m}_0^{\beta }exp\left[-\beta \gamma V_2(t_0,h_0)\right]}{{\theta }^{\beta }{e}^{-\beta \left(r(T-t_0)\right)}+\beta m_0^{\beta }r\underset{t_0}{\overset{T}{\int }}exp\left[-\beta (r(\tau -t_0)+\gamma V_2(\tau ,h^0\left(\tau \right)))\right]d\tau }}\\ {\displaystyle =\frac{{\theta }^{\beta }{m}_0^{\beta }exp\left[\underset{t_0}{\overset{T}{\int }}\beta [r-\gamma f\left(h^0\left(\tau \right)\right)]d\tau \right]}{{\theta }^{\beta }+\beta m_0^{\beta }r\underset{t_0}{\overset{T}{\int }}exp\left[{\displaystyle \beta \underset{\tau }{\overset{T}{\int }}[r-\gamma f\left(h^0\left(y\right)\right)]dy}\right]d\tau },}\end{array}$$

(20)

where $h^0(·)$ is the solution of system (1) that is generated by the optimal feedback control ($u^0(t,h)$) (16).

Note that $(t_0,h_0,m_0)\in W_2$ if

$$m_2^0(T;t_0,h_0,m_0,u^0(·))={\left[\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)\right]}^{-\frac{1}{\beta }}\le M_2.$$

Thus, it follows from (20) that $m_0\in W_2$ if

$$\begin{array}{c}{\displaystyle m_0\le \frac{\theta M_2}{{\left[M_2^{\beta }\beta r\underset{t_0}{\overset{T}{\int }}exp\left[\beta \underset{\tau }{\overset{T}{\int }}\left[r-\gamma f\left(h^0\left(y\right)\right)\right]dy\right]d\tau -{\theta }^{\beta }exp\left[\beta \underset{t_0}{\overset{T}{\int }}\left[r-\gamma f\left(h^0\left(\tau \right)\right)\right]d\tau \right]\right]}^{\frac{1}{\beta }}}}\\ \triangleq \tilde{m}\left(h^0(·)\right).\end{array}$$

(21)

We now successively consider the regions (14) of set $[0,T]\times [0,L]$ to describe the construction of set $W_2$.

(1)

Consider the following points:

$$(t_0,h_0,m_0)\subset G_1\times [0,M_2],$$

(22)

where set $G_1$ is defined in (14).

It follows from the form of the optimal control (16) that in system (1),

$$u^0(t,h)=\alpha {\widehat{h}}_1,h^0\left(t\right)=h^0(t;t_0,h_0,u^0(·))=h^0(t;t_0,{\widehat{h}}_1,\alpha {\widehat{h}}_1)\equiv {\widehat{h}}_1.$$

Then, it follows from (21) that

$$\begin{array}{c}\{(t_0,h_0,m_0):(t_0,h_0)\in G_1,m_0\in [0,\tilde{m}\left({\widehat{h}}_1\right)]\}\subset W_2.\end{array}$$

(23)

(2)

Consider the following points:

$$(t_0,h_0,m_0)\subset G_2\times [0,M_2],$$

(24)

where set $G_2$ is defined in (14).

It follows from the form of the optimal control (16) that in system (1),

$$u^0(t,h)=\alpha {\widehat{h}}_3,h^0\left(t\right)=h^0(t;t_0,h_0,u^0(·))=h^0(t;t_0,{\widehat{h}}_3,\alpha {\widehat{h}}_3)\equiv {\widehat{h}}_3.$$

Then, it follows from (21) that

$$\begin{array}{c}\{(t_0,h_0,m_0):(t_0,h_0)\in G_2,m_0\in [0,\tilde{m}\left({\widehat{h}}_3\right)]\}\subset W_2.\end{array}$$

(25)

(3)

Consider the following points:

$$(t_0,h_0,m_0)\subset {\mathrm{\Pi }}_1\times [0,M_2],$$

(26)

where set ${\mathrm{\Pi }}_1$ is defined in (14).

There exist two cases:

(3.a) The graph of $h^0\left(t\right)$ does not reach ${\widehat{h}}_1$ on $[t_0,T]$.

Then, it follows from the form of the optimal control (16) that the solution of (1) is

$$h^0(t;t_0,h_0,u^0(·))=h^0(t;t_0,h_0,Q)=\left(h_0-\frac{Q}{\alpha }\right)e^{-\alpha (t-t_0)}+\frac{Q}{\alpha }.$$

(27)

It follows from (27) that case 3.a occurs when

$$h_0<h_3^{*}\left(t_0\right)\triangleq \left({\widehat{h}}_1-\frac{Q}{\alpha }\right)e^{\alpha (T-t_0)}+\frac{Q}{\alpha }.$$

Therefore, it follows from (21) that

$$\begin{array}{c}\{(t_0,h_0,m_0):t_0\in [t_0,T],h_0\in [0,h_3^{*}\left(t_0\right)),\\ m_0\in \left[0,\tilde{m}\left(\left(h_0-\frac{Q}{\alpha }\right)e^{-\alpha (t-t_0)}+\frac{Q}{\alpha }\right)\right]\}\subset W_2.\end{array}$$

(28)

(3.b) The graph of $h^0\left(t\right)$ reaches ${\widehat{h}}_1$:

$$\exists t_3^{*}\in (t_0,T]:h^0(t_3^{*};t_0,h_0,u^0(·))=h^0(t_3^{*};t_0,h_0,Q)={\widehat{h}}_1.$$

It follows from (27) that this case occurs when

$$h_0\ge h_3^{*}\left(t_0\right)$$

and

$$t_3^{*}=\frac{1}{\alpha }ln\frac{h_0-\frac{Q}{\alpha }}{{\widehat{h}}_1-\frac{Q}{\alpha }}+t_0.$$

It follows from the optimality principle that

$$\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)=\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_3^{*},{\widehat{h}}_1,m_2^{*}\left(t_3^{*}\right)),$$

where

$$m_2^{*}\left(t\right)\triangleq m_2^0(t;t_0,h_0,m_0,u^0(·))=m_2^0(t;t_0,h_0,m_0,Q).$$

Therefore, in case 3.b, $(t_0,h_0,m_0)\in W_2$ when

$$m_2^0(T;t_0,h_0,m_0,u^0(·))=m_2^0(T;t_3^{*},{\widehat{h}}_1,m_2^{*}\left(t_3^{*}\right),\alpha {\widehat{h}}_1)\le M_2.$$

However, it follows from case 1 that the latter expression is true when

$$(t_3^{*},{\widehat{h}}_1,m_2^{*}\left(t_3^{*}\right))\in G_1\times [0,\tilde{m}\left({\widehat{h}}_1\right)].$$

Therefore,

$$\begin{array}{c}\{(t_0,h_0,m_0):t_0\in [0,T],h_0\in [h_3^{*}\left(t_0\right),{\widehat{h}}_1),m_2^{*}\left(t_3^{*}\right)\in [0,\tilde{m}({\widehat{h}}_1))]\}\subset W_2.\end{array}$$

(29)

Remark 2. 

Hereafter, the value of $m_2^0(t;t_0,h_0,m_0)$ can be calculated according to Formula (7).

(4)

Consider the following points:

$$(t_0,h_0,m_0)\subset {\mathrm{\Pi }}_2\times [0,M_2],$$

(30)

where set ${\mathrm{\Pi }}_2$ is defined in (14).

Similarly to case 3, there exist two cases:

(4.a) The graph of $h^0\left(t\right)$ does not reach ${\widehat{h}}_3$ on $[t_0,T]$.

Then, it follows from the form of the optimal control (16) that the solution of (1) is

$$h^0(t;t_0,h_0,u^0(·))=h^0(t_4^{*};t_0,h_0,0)=h_0{e}^{-\alpha (t-t_0)}.$$

(31)

It follows from (31) that this case occurs when

$$h_0>h_4^{*}\left(t_0\right)\triangleq {\widehat{h}}_3{e}^{\alpha (T-t_0)}.$$

Therefore, it follows from (21) that

$$\begin{array}{c}\{(t_0,h_0,m_0):t_0\in [t_0,T],h_0\in (h_4^{*}\left(t_0\right),L],m_0\in \left[0,\tilde{m}\left(h_0{e}^{-\alpha (t-t_0)}\right)\right]\}\subset W_2.\end{array}$$

(32)

(4.b) The graph of $h^0\left(t\right)$ reaches ${\widehat{h}}_3$ on $[t_0,T]$:

$$\exists t_4^{*}\in (t_0,T]:h^0(t_4^{*};t_0,h_0,u^0(·))=h^0(t_4^{*};t_0,h_0,0)={\widehat{h}}_3.$$

It follows from (31) that this case occurs when

$$h_0\le h_4^{*}\left(t_0\right)$$

and

$${\displaystyle t_4^{*}=\frac{1}{\alpha }ln\frac{h_0}{{\widehat{h}}_1}+t_0.}$$

It follows from the optimality principle that

$$\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)=\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_4^{*},{\widehat{h}}_3,m_2^{*}\left(t_4^{*}\right)),$$

where

$$m_2^{*}\left(t\right)=m_2^0(t;t_0,h_0,m_0,u^0(·))=m_2^0(t;t_0,h_0,m_0,0).$$

Therefore, in case 4.b, $(t_0,h_0,m_0)\in W_2$ when

$$m_2^0(T;t_0,h_0,m_0,u^0(·))=m_2^0(T;t_4^{*},{\widehat{h}}_3,m_2^{*}\left(t_4^{*}\right),\alpha {\widehat{h}}_3)\le M_2.$$

However, it follows from case 1 that the latter expression is true when

$$(t_4^{*},{\widehat{h}}_3,m_2^{*}\left(t_4^{*}\right))\in G_2\times [0,\tilde{m}\left({\widehat{h}}_3\right)].$$

Therefore,

$$\begin{array}{c}\{(t_0,h_0,m_0):t_0\in [0,T],h_0\in [h_4^{*}\left(t_0\right),{\widehat{h}}_3),m_2^0\left(t_4^{*}\right)\in [0,\tilde{m}\left({\widehat{h}}_3\right))]\}\subset W_2.\end{array}$$

(33)

(5)

Consider the following points:

$$(t_0,h_0,m_0)\in {\mathrm{\Pi }}_3\times [0,M_2],$$

(34)

where set ${\mathrm{\Pi }}_3$ is defined in (14).

There exist two cases:

(5.a) The graph of $h^0\left(t\right)$ does not reach ${\widehat{h}}_1$ on $[t_0,T]$.

Then, it follows from the form of the optimal control (16) that the solution of (1) is

$$h^0(t;t_0,h_0,u^0(·))=h^0(t;t_0,h_0,0)=h_0{e}^{-\alpha (t-t_0)}.$$

(35)

It follows from (35) that this case occurs when

$$h_0>h_5^{*}\left(t_0\right)\triangleq {\widehat{h}}_1{e}^{\alpha (T-t_0)}.$$

Therefore, it follows from (21) that

$$\begin{array}{c}\{(t_0,h_0,m_0):t_0\in [t_0,T],h_0\in (h_5^{*}\left(t_0\right),x\left(t_0\right)],m_0\in \left[0,\tilde{m}\left(h_0{e}^{-\alpha (t-t_0)}\right)\right]\}\subset W_2,\end{array}$$

(36)

where $x\left(t_0\right)$ is the starting point of $\mathrm{\Gamma }$ curve (see Figure 1).

(5.b) The graph of $h^0\left(t\right)$ reaches ${\widehat{h}}_1$ on $[t_0,T]$:

$$\exists t_5^{*}\in (t_0,T]:h^0(t_5^{*};t_0,h_0,u^0(·))=h^0(t;t_0,h_0,0)={\widehat{h}}_1.$$

It follows from (35) that this case occurs when

$$h_0\le h_5^{*}\left(t_0\right)$$

and

$${\displaystyle t_5^{*}=\frac{1}{\alpha }ln\frac{h_0}{{\widehat{h}}_1}+t_0.}$$

It follows from the optimality principle that

$$\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)=\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_5^{*},{\widehat{h}}_1,m_2^{*}\left(t_5^{*}\right)),$$

where

$$m_2^{*}\left(t\right)=m_2^0(t;t_0,h_0,m_0,u^0(·))=m_2^0(t;t_0,h_0,m_0,0).$$

Therefore, in case 2.a, $(t_0,h_0,m_0)\in W_2$ when

$$m_2^0(T;t_0,h_0,m_0,u^0(·))=m_2^0(T;t_5^{*},{\widehat{h}}_1,m_2^{*}\left(t_5^{*}\right),\alpha {\widehat{h}}_1)\le M_2.$$

However, it follows from case 1 that the latter expression is true when

$$(t_5^{*},{\widehat{h}}_1,m_2^{*}\left(t_5^{*}\right))\in G_1\times [0,\tilde{m}\left({\widehat{h}}_1\right)].$$

Therefore,

$$\begin{array}{c}\{(t_0,h_0,m_0):t_0\in [0,T],h_0\in [h_5^{*}\left(t_0\right),{\widehat{h}}_1),m_2^{*}\left(t_5^{*}\right)\in [0,\tilde{m}\left({\widehat{h}}_1\right))]\}\subset W_2.\end{array}$$

(37)

(6)

Consider the following points:

$$(t_0,h_0,m_0)\subset {\mathrm{\Pi }}_4\times [0,M_2],$$

(38)

where set ${\mathrm{\Pi }}_4$ is defined in (14).

There exist two cases:

(6.a) The graph of $h^0\left(t\right)$ does not reach ${\widehat{h}}_3$ on $[t_0,T]$.

Then, it follows from the form of the optimal control (16) that the solution of (1) is

$$h^0(t;t_0,h_0,u^0(·))=h^0(t;t_0,h_0,Q)=\left(h_0-\frac{Q}{\alpha }\right)e^{-\alpha (t-t_0)}+\frac{Q}{\alpha }.$$

(39)

It follows from (39) that this case occurs when

$$h_0<h_6^{*}\left(t_0\right)\triangleq \left({\widehat{h}}_3-\frac{Q}{\alpha }\right)e^{\alpha (T-t_0)}+\frac{Q}{\alpha }.$$

Therefore, it follows from (21) that

$$\begin{array}{c}\{(t_0,h_0,m_0):t_0\in [t_0,T],h_0[0,h_6^{*}\left(t_0\right)),\\ m_0\in \left[0,\tilde{m}\left(\left(h_0-\frac{Q}{\alpha }\right)e^{-\alpha (t-t_0)}+\frac{Q}{\alpha }\right)\right]\}\subset W_2.\end{array}$$

(40)

(6.b) The graph of $h^0\left(t\right)$ reaches ${\widehat{h}}_3$:

$$\exists t_6^{*}\in (t_0,T]:h^0(t_6^{*};t_0,h_0,u^0(·))=h^0(t;t_0,h_0,Q)={\widehat{h}}_3.$$

It follows from (39) that this case occurs when

$$h_0\ge h_6^{*}\left(t_0\right)$$

and

$$t_6^{*}=\frac{1}{\alpha }ln\frac{h_0-\frac{Q}{\alpha }}{{\widehat{h}}_3-\frac{Q}{\alpha }}+t_0.$$

It follows from the optimality principle that

$$\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)=\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_6^{*},{\widehat{h}}_3,m_2^{*}\left(t_6^{*}\right)),$$

where

$$m_2^{*}\left(t\right)=m_2^0(t;t_0,h_0,m_0,u^0(·))=m_2^0(t;t_0,h_0,m_0,Q).$$

Therefore, in case 6.b, $(t_0,h_0,m_0)\in W_2$ when

$$m_2^0(T;t_0,h_0,m_0,u^0(·))=m_2^0(T;t_6^{*},{\widehat{h}}_3,m_2^{*}\left(t_6^{*}\right),\alpha {\widehat{h}}_3)\le M_2.$$

However, it follows from case 1 that the latter expression is true when

$$(t_6^{*},{\widehat{h}}_3,m_2^{*}\left(t_6^{*}\right))\in G_1\times [0,\tilde{m}\left({\widehat{h}}_3\right)].$$

Therefore,

$$\begin{array}{c}\{(t_0,h_0,m_0):t_0\in [0,T],h_0\in [h_6^{*}\left(t_0\right),{\widehat{h}}_3),m_2^{*}\left(t_6^{*}\right)\in [0,\tilde{m}\left({\widehat{h}}_3\right))]\}\subset W_2.\end{array}$$

(41)

Note that it is obvious that

$$W_2{|}_{t=T}=\{(T,h_0,m_0):(T,h_0,m_0)\in W_2\}=T\times [0,L]\times [0,M_2].$$

Therefore, expressions (23), (25), (28), (29), (32), (33), (36), (37), (40) and (41) describe the structure of viability set $W_2$ (19) for problem (1), (6) according to the generalized logistic law.

Remark 3. 

It was proven in [20] that the viability set for problem (1), (5) according the Gompertz law has the following form:

$$W_1=\left\{(t_0,h_0,m_0)\in [0,T]\times [0,M_1]\times [0,L]:{\mathrm{V}\mathrm{a}\mathrm{l}}_1(t_0,h_0,m_0)\le {M_1}^2\right\},$$

$$M_1=e^{{\displaystyle \frac{r-\gamma F}{\theta }}}.$$

3.2. The Maximal Viability Set

Set $W_2$, as described in the previous section, is the maximal viability set for problem (1), (6), and the following theorem is true.

Theorem 1. 

The following statements are true:

1.

$\forall (t_0,h_0,m_0)\in W_2{m}_0\le M_2$.

2.

Set $W_2$ (19) is weakly invariant with respect to the differential inclusion: $\dot{w}\in Y\left(w\right)$, where

$$w=(t,h,m)\to Y\left(w\right)={\left(1,g_2\left(m\right)-\gamma mf\left(h\right),-\alpha h+[0,Q]\right)}^{\top }.$$

(42)

3.

Set $W_2$ is the maximal viability set for problem (1), (6).

Proof. 

(1) Consider a point ($(t_0,h_0,m_0)\in W_2$). Let

$$m_2^0\left(T\right)=m_2^0(T;t_0,h_0,m_0,u^0(·))\le M_2.$$

Then, one can express $m_0$ from (20) as:

$$m_2^0\left(T\right)=\frac{\theta m_{0}exp\left[{\displaystyle \underset{t_0}{\overset{T}{\int }}{F}^0\left(\tau \right)d\tau }\right]}{{\left[{\displaystyle {\theta }^{\beta }+\beta m_0^{\beta }r\underset{t_0}{\overset{T}{\int }}exp\left[{\displaystyle \beta \underset{\tau }{\overset{T}{\int }}{F}^0\left(y\right)dy}\right]d\tau }\right]}^{\frac{1}{\beta }}}\le M_2,$$

where $F^0\left(\tau \right)\triangleq r-\gamma f\left(h^0\left(\tau \right)\right)$, and $h^0\left(t\right)$ is the solution of (1) generated by the optimal feedback ($u^0\left(t\right)=u^0(t,h)$) (16).

Also note the expression for $M_2$ (18):

$$M_2=\theta {\left[1-\frac{\gamma F}{r}\right]}^{\frac{1}{\beta }},$$

(43)

where $F=f\left({\widehat{h}}_1\right)=f\left({\widehat{h}}_3\right)=\underset{h\in [0,L]}{max}f\left(h\right)$. Thus,

$${\displaystyle {m_2^0\left(T\right)}^{\beta }\left({\theta }^{\beta }+\beta m_0^{\beta }r\underset{t_0}{\overset{T}{\int }}exp\left[\underset{\tau }{\overset{T}{\int }}\beta F^0\left(y\right)dy\right]d\tau \right)={\theta }^{\beta }{m}_0^{\beta }exp\left[\underset{t_0}{\overset{T}{\int }}\beta F^0\left(\tau \right)d\tau \right],}$$

and

$${\displaystyle m_0^{\beta }=\frac{{m_2^0\left(T\right)}^{\beta }{\theta }^{\beta }}{{\theta }^{\beta }exp\left[\underset{t_0}{\overset{T}{\int }}\beta F^0\left(\tau \right)d\tau \right]-{m_2^0\left(T\right)}^{\beta }\left(\beta r\underset{t_0}{\overset{T}{\int }}exp\left[\underset{\tau }{\overset{T}{\int }}\beta F^0\left(y\right)dy\right]d\tau \right)}.}$$

(44)

The following estimates are true:

$$\begin{array}{c}{F}^0\left(y\right)=r-\gamma f\left(h^0\left(y\right)\right)\ge r-\gamma F>0,\beta F^0\left(y\right)\ge \beta (r-\gamma F),\\ exp\left[{\displaystyle \beta \underset{\tau }{\overset{T}{\int }}{F}^0\left(y\right)dyd\tau }\right]\ge exp\left[\beta (r-\gamma F)(T-\tau )\right].\end{array}$$

(45)

Apply (43) and (45) to (44):

$$\begin{array}{c}{\displaystyle m_0^{\beta }\le \frac{{M_2}^{\beta }{\theta }^{\beta }}{{\theta }^{\beta }exp\left[\underset{t_0}{\overset{T}{\int }}\beta F^0\left(\tau \right)d\tau \right]-{\theta }^{\beta }\left(1-\frac{\gamma F}{r}\right)\beta r\underset{t_0}{\overset{T}{\int }}exp\left[\underset{\tau }{\overset{T}{\int }}\beta F^0\left(y\right)dy\right]d\tau }}\\ {\displaystyle \le \frac{{M_2}^{\beta }}{e^{\beta (r-\gamma F)(T-t_0)}-(r-\gamma F)\beta \underset{t_0}{\overset{T}{\int }}{e}^{\beta (r-\gamma F)(T-\tau )}d\tau }={M_2}^{\beta }.}\end{array}$$

(46)

Therefore, for any $(t_0,h_0,m_0)\in W_2$, $m_0\le M_2$.

(2) It follows from the construction of $W_2$ that all the trajectories ($m_2^0(t;t_0,h_0,m_0,u(·))$) that start from the inner points ($(t_0,h_0,m_0)$) of $W_2$ stay inside this set for any $t\in [t_0,T]$.

We assume that such a trajectory corresponds to a point $(t^{*},h^{*},m^{*})$ on the boundary of $W_2$. Then, it follows from the optimality principle that

$$\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t^{*},h^{*},m^{*})=\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t,h^0\left(t\right),m_2^0\left(t\right))\le M_2,t\in [t^{*},T].$$

In other words, the trajectory does not leave $W_2$. Thus, set $W_2$ is weakly invariant with respect to the differential inclusion (42).

(3) Consider the point $(\overline{t},\overline{h},\overline{m})\notin W_2$. Then, it follows from (19) that the value function $\mathrm{V}\mathrm{a}{\mathrm{l}}_2(\overline{t},\overline{h},\overline{m})$ for the optimal control $u^0(·)$ is greater than $M_2$:

$$m_2^0(T;\overline{t},\overline{h},\overline{m},u^0(·))>M_2.$$

However, due to the continuity of the trajectory ($m_2^0(t;\overline{t},\overline{h},\overline{m},u^0(·))$), there exists $t^{*}\in [\overline{t},T)$ such that

$$m_2^0(t^{*};\overline{t},\overline{h},\overline{m},u^0(·))>M_2.$$

For any point ($(t_0,h_0,m_0)\notin W_2$) and any measurable function ($u(·):[t_0,T]\to [0,Q]$), according to the definition of the optimal result ($\mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)$),

$${\sigma }_2\left(m_2(T;t_0,h_0,m_0,u\left(t\right))\right)\ge \mathrm{V}\mathrm{a}{\mathrm{l}}_2(t_0,h_0,m_0)>M_2.$$

Again, it follows from the continuity of the trajectory ($m_2(t;t_0,h_0,m_0,u(·))$) that

$$m_2(t^{*};t_0,h_0,m_0,u(·))>M_2,t^{*}<T.$$

In other words, the tumor volume exceeds the threshold compatible with life at moment $t^{*}$, which is earlier than the control point (T).

Therefore, $W_2$ is the maximal viability set for problem (1), (6). □

3.3. The Inverse Problem

For models of chemotherapy for malignant tumors, it is often of interest to solve inverse problems, for example, to identify the growth parameter (r) in the model of this tumor and to reconstruct the treatment protocol ($u(·)$) for the patient using the current measurements of the tumor volume ($m\left(t\right)$) and the amount of drug ($h\left(t\right)$) in the body of the patient.

Consider a model according to the Gompertz law:

$$\begin{array}{c}{\displaystyle \frac{dm}{dt}=g_1\left(m\right)-\gamma mf\left(h\right),}\\ {\displaystyle \frac{dh}{dt}=-\alpha h+u\left(t\right),}\\ t\in [0,T],m\left(0\right)=m_0,h\left(0\right)=h_0.\end{array}$$

(47)

The known parameters are:

$$\begin{array}{c}{g}_1\left(m\right)=rm-\theta mln\left(m\right),\\ f\left(h\right)=0.15(sin[2.2h-0.5\pi ]+1),\\ \theta ={10}^{-2},\gamma =1,\alpha =0.03,\\ T=8,m_0=5*{10}^7,h_0=0.\end{array}$$

The feasible controls are measurable (piecewise constant) functions satisfying the following restrictions:

$$u\left(t\right)\in \mathbf{U}=[0,Q].$$

(48)

The parameters to reconstruct are:

$$r,u(·).$$

(49)

It is assumed that the process ($m^{*}(·),h^{*}(·)$) generated by some unknown feasible control ($u^{*}(·)$) is observed at the following time interval:

$$t\in [0,T],0<T<\infty .$$

Discrete inaccurate measurements ($m_i^{\delta },h_i^{\delta }$) of the process arrive in real time with a uniform time step ($▵t$). The measurement error is $\delta =({\delta }_m,{\delta }_h)$:

$$\begin{array}{c}\begin{array}{c}|m_i^{\delta }-m^{*}\left(t_i\right)|\le {\delta }_m,\hfill \\ |h_i^{\delta }-h^{*}\left(t_i\right)|\le {\delta }_h\hfill \end{array},t_i=i▵t,i=0,\dots ,N,N=\lceil T/▵t\rceil .\end{array}$$

(50)

The inverse problem is to construct approximations ($r^{\delta }$) of parameter r and piecewise constant approximating functions ($u^{\delta }(·)$) such that:

• $u^{\delta }\left(t\right)\in \mathbf{U},t\in [0,T]$;
• $u^{\delta }\left(t\right)\stackrel{\delta \to 0,▵t\to 0}{\to }{u}^{*}\left(t\right)$ almost everywhere on $[0,T]$;
• ${\displaystyle r^{\delta }\stackrel{\delta \to 0,▵t\to 0}{\to }r}$;
• The trajectories ($m^{\delta }(·),h^{\delta }(·)$) of the system (47) that are generated by the reconstructed parameters ($r^{\delta }$ and $u^{\delta }(·)$) converge uniformly to the observed process:

  $$\begin{array}{c}\underset{t\in [0,T]}{max}|m^{\delta }\left(t\right)-m^{*}\left(t\right)|\stackrel{\delta \to 0,▵t\to 0}{\to }0,\\ \underset{t\in [0,T]}{max}|h^{\delta }\left(t\right)-h^{*}\left(t\right)|\stackrel{\delta \to 0,▵t\to 0}{\to }0.\end{array}$$

  (51)

An approach for constructing such approximations was suggested and justified in [21] using the constructions from auxiliary problems of calculus of variations. These constructions are the solutions of Hamiltonian systems. The feature of the approach is the use of non-convex auxiliary variational problems, namely, the integrands of the functionals in the auxiliary problems are d. c. functions [23]. Such an approach proves the stability of the methods based on this approach with respect to the measurement and approximation errors. This approach reduces the inverse problems for dynamical systems to integration of linear ODE systems.

In this paper, we present the results of the numerical simulation for a method based on this approach, as explained and justified in [21].

This method was used in two numerical simulations. In the simulations, the parameters to reconstruct were:

$$r=0.18,u^{*}\left(t\right)=\left\{\begin{array}{cc}1,\hfill & t\in [0,4)\bigcup [6,8],\hfill \\ 0.1,\hfill & t\in [4,6).\hfill \end{array}\right.$$

(52)

The graph of the therapy function $f\left(h\right)$ is shown in Figure 2.

To obtain the inaccurate measurements ($m_i^{\delta },h_i^{\delta }$) of process $m^{*}(·),h^{*}(·)$ generated by $u^{*}(·)$, this process was constructed numerically and randomly perturbed with uniform distribution for two sets of the parameters of observation (the measurement error ($\delta$) and the time step ($▵t$)). In the first simulation, $({\delta }_m,{\delta }_h)=(5\ast {10}^5,3\ast {10}^{-2}),▵t=T/20=0.4$. In the second simulation, $({\delta }_m,{\delta }_h)=({10}^5,0.6\ast {10}^{-2}),▵t=T/100=0.08$. The approximations of the presumably unknown parameters (52) were constructed based on these sets of data using the method proposed in [21].

The result of the identification in the first simulation are presented as follows. The approximation ($r^{\delta }$) of parameter r (52), which characterizes the speed of the tumor’s growth, is

$$r^{\delta }\approx 0.152,|r-r^{\delta }|\approx 0.028.$$

The results of the reconstruction of the control $u^{*}(·)$ in the first simulation are shown in Figure 3.

In this figure (and in Figure 4), the black graph is the control ($u^{*}\left(t\right)$) (52) to be reconstructed. This control was applied to the model’s system (47) to simulate the inaccurate measurements of the basic (observed) trajectory. The red graph is the piecewise constant control ($u^{\delta }\left(t\right)$), which approximates the reconstructed control ($u^{*}\left(t\right)$). It was constructed using the variational method suggested in [21].

The process ($m^{\delta }\left(t\right),h^{\delta }\left(t\right)$) generated by $r^{\delta }$ and $u^{\delta }(·)$ in the first simulation is shown in Figure 5. In this figure (and in Figure 6), the blue graph is the “observed” process ($\{m^{*}\left(t\right),h^{*}\left(t\right)\}$). It is the trajectory of the model’s system (47) generated by the “unknown” control ($u^{*}(·)$) (52) (to be reconstructed). It is supposed that the inaccurate measurements of this trajectory are known. To simulate them, this trajectory was perturbed at discrete points with random error (with uniform distribution). These data are indicated by black dots. The red graph is the process $\{m^{\delta }\left(t\right),h^{\delta }\left(t\right)\}$. It is the trajectory of the model’s system (47) generated by the control ($u^{\delta }(·)$) (52) (which is the constructed approximation of the “unknown” control ($u^{*}(·)$)). The aim of these figures is to show that the reconstructed trajectory of the system remains close to the observed trajectory (see condition 4 (51) in the inverse problem statement).

The results of identification in the second simulation are reported as follows. The approximation ($r^{\delta }$) of parameter r (52), which characterizes the speed of the tumor’s growth, is

$$r^{\delta }\approx 0.174,|r-r^{\delta }|\approx 0.006.$$

The results of reconstruction of the control ($u^{*}(·)$) in the second simulation are shown in Figure 4. In this figure, the black graph is the control ($u^{*}\left(t\right)$) (52) to be reconstructed. The red graph is the piecewise constant control ($u^{\delta }\left(t\right)$), which approximates the reconstructed control ($u^{*}\left(t\right)$).

The process ($m^{\delta }(·),h^{\delta }(·)$) generated by $r^{\delta }$ and $u^{\delta }(·)$ in the second simulation is shown in Figure 6. In this figure, the blue graph is the “observed” process ($\{m^{*}\left(t\right),h^{*}\left(t\right)\}$). The black dots are the inaccurate measurements of this process. The red graph is the process ($\{m^{\delta }\left(t\right),h^{\delta }\left(t\right)\}$) generated by the control ($u^{\delta }(·)$) (52). The aim of these figures is to show that the reconstructed trajectory of the system remains close to the observed trajectory (see condition 4 (51) in the inverse problem statement).

4. Discussion

In this paper, models of chemotherapy for malignant tumors were studied. In particular, the problem of constructing the optimal control with the aim of minimizing the tumor volume at a given time point was considered. Possible directions for future research may include minimization of the total cost of the drug injected into the patient’s body. For such cases, the cost functional (which is a treatment quality criterion) should have the following form:

$$\underset{t_0}{\overset{T}{\int }}{u}^2\left(t\right)dt+\sigma \left(m\left(T\right)\right).$$

Another problem to consider in the future is minimization of the total volume of the drug in the body for the entire period of treatment, i.e., minimization of the quality criterion

$$\underset{t_0}{\overset{T}{\int }}{h}^2\left(t\right)dt+\sigma \left(m\left(T\right)\right).$$

All these problems can be solved by following the scheme proposed in this paper. The viability sets in these problems can also be constructed.

The models of chemotherapy considered in this paper are the simplest ones. They were previously studied in [3] for the case of a therapy function that is either linear or has a single maximum point. In this paper, we studied a model with a therapy function that has two maximum points. It describes a case in which two different therapeutic agents are used. One maximum point corresponds to maximum efficiency of the first agent and the second to maximum efficiency of the second agent.

The future development of the research presented in this paper will include the study of models tumor growth, which include the dynamics of the healthy cells and the body’s immune system and their interaction with the cancer cells and the drug (see, for example, [6,7,16,17]).

The results of the numerical construction of the viability sets will be presented in the future works.

It should be noted that the studied models are simplified idealizations of real processes. Therefore, the obtained results are informative and recommendatory in nature and require further experimental verification.

5. Conclusions

This paper presents the results of a study of a chemotherapy model of a malignant tumor growing according to the Gompertz law and the generalized logistic law. The optimal therapy strategy and optimal treatment protocols, which were constructed in previous works, were discussed with the aim of minimizing the tumor volume at the control point. We also discussed managerial insights with respect to their structure.

In the model of chemotherapy, a non-monotonic therapy function with two maxima is considered. This function defines the effect of the therapeutic agents on the tumor. Such a therapy function leads to a specific form of the optimal controls. They are piecewise constant functions with no more than one switch. The specifics are that the Rankine-Hugoniot line plays an important role in the construction of the optimal protocols. The position of the initial state of the model with respect to this line affects the construction of the optimal treatment protocols.

The structure of the maximal viability set for the considered model with the generalized logistic law is described. It is the set of the initial states of the model (the tumor volume and the concentration of the drug in the body) such that there exists a treatment protocol (a control) that guarantees that the tumor volume will not exceed the threshold amount compatible with life.

An inverse problem of reconstruction of the protocol of the ongoing treatment and identification of the parameter of the intensity of the tumor’s growth based on the observations of the dynamics of the treatment process is considered. Such problems have never been considered before. A solution to this problem is suggested. The results of the numerical simulations for this problem are presented herein.