1. Introduction

Anti-angiogenic therapies are cancer treatments, proposed at first by Folkman [13] in the early seventies, consolidated along the nineties by several discoveries on the main principles regulating tumor angiogenesis [20] and widely debated in several theoretical and experimental studies throughout the last decade. Anti-angiogenic treatments aim at inhibiting the development of the vascular network necessary to support tumor growth during the vascular phase, so providing a way to control the heterogeneous and growth-unconstrained tumor population throughout the control of the homogeneous and growth-constrained population of endothelial cells [7,12]. Tumors have the capability to develop resistance to the conventional chemotherapeutic drugs mainly because of the rapidity of tumor cells in evolving towards new resistant phenotypes. Due to the indirect action of the anti-angiogenic drugs, the outcome of the therapy should not be impaired by the capability of tumor cells to generate resistant phenotypic variants [16,17], and the effectiveness of anti-angiogenic treatments on the control and possible remission of experimental tumors has been demonstrated [17]. Moreover, anti-angiogenic therapies have limited side effects respect to the conventional chemotherapies and radiotherapies. Conventional chemotherapies may also have anti-angiogenic effects on the vascular network [1,18,19].

The scope of this note is to investigate and design a closed-loop model-based anti-angiogenic therapy. The adopted model refers to [14], where a quantitative model describing the growth of experimental tumors under the control of the vascular network is presented. In [14], Hahnfeldt and coworkers introduced the concept of the carrying capacity of the vasculature, i.e. the tumor volume potentially sustainable by the vasculature, in order to account for the vascular control on the tumor growth. As the carrying capacity is strictly dependent on the vasculature extension, its dynamics can be assumed to represent the dynamics of the vascular network. Moreover, the effect on the carrying capacity of stimulatory and inhibitory angiogenic signals produced by the tumor itself and of administered anti-angiogenic drugs are explicitly accounted by the model formulation. The paper represents one of the first attempts to model, with a minimal number of parameters, the interplay between the dynamics of the tumor volume and of the carrying capacity, with or without administration of anti-angiogenic drugs. In [14] the predictions of the model have also been successfully compared with experimental data on anti-angiogenically treated and untreated Lewis lung tumors in mice.

Besides experimental frameworks, the model proposed in [14] has been widely exploited in theoretical studies in order to predict the effectiveness of new anti-angiogenic therapies and some model extensions have also been proposed in the related scientific literature. In [8] some model modifications are proposed and conditions for the eradication of the tumor under a periodic anti-angiogenic treatment are provided. In [23,24] problems addressing the optimal scheduling of a given amount of angiogenic inhibitors are presented. Anti-angiogenic therapies are nowadays investigated for their use in combination with chemotherapy, within the framework of a multidrug cancer therapy, and the Hahnfeldt model is still a landmark to synthesize and validate both combined or pure anti-angiogenic drug therapies. In [11] the optimal scheduling problem of a combined radiotherapy and anti-angiogenic treatment is formulated by exploiting a suitable modified version of the original model. Further extensions are also proposed in [9], where the authors aim at describing the interplay between the populations of tumor cells and endothelial cells subject to a combined therapy of chemotherapeutic and anti-angiogenic drugs. In [10] an optimal scheduling problem for the combined chemotherapeutic/anti-angiogenic treatment is presented. The reader may refer also to the very recent [21], where a modified version of [14] accounting for multiple control delays related to pharmacodynamic/pharmacokinetic drug absorption is exploited, in order to investigate local controllability and optimal control on a finite treatment horizon utilizing anti-angiogenic therapy combined to chemotherapy. On the other hand, pure anti-angiogenic approaches are still investigated in the recent literature, like in [22,26], where optimal and robust control strategies are applied to the linearized version of [14].

By suitably exploiting the Hahnfeldt model [14], tumor volume reduction is here pursued according to a closed-loop, model-based approach, with the control strategy making use of the feedback linearization theory [15]. To this end only available tumor volume measurements are exploited, with the carrying capacity estimated by means of a state observer for nonlinear systems [4,6], and the regulator synthesized as a feedback from the observed state. As known, the purpose of anti-angiogenic treatment is to keep the tumor volume below a safe level and for this reason the closed-loop control law is designed to track a desired volume level smaller than the safety value, possibly starting from a high level situation.

Preliminary results have been presented in [5], where a continuous stream of measurements has been considered. Instead, in a realistic therapeutic setting, measurements are acquired at discrete sampling times. In order to comply with this issue without resorting to an approximate discretization scheme, we exploit recent developments in the theory of time-delay systems [2] where sampled measurements are modeled as continuous-time outputs affected by time-varying delays [25], according to which the continuous-discrete control problem can be reformulated as a continuous problem with delayed measurements. We pursue this approach to design the state observer for the tumor-growth control algorithm.

It is proven that the control scheme allows to set independently the control and the observer parameters, thanks to the structural properties of the tumor growth model that guarantee the separability of estimation and feedback control algorithms.

Simulations are carried out in order to mimic a real experimental framework on mice. These results seem extremely promising with respect to both the reduction of the tumor growth level and to the bound on the average drug administration, required to make feasible the proposed drug delivery therapy. More in details, simulations provide very good performances according to the measurements sampling intervals suggested by the experimental literature, and show a noticeable level of robustness against the observer initial estimate, as well as against the uncertainties affecting the model parameters.

The paper is organized as follows. The next section is devoted to detail the tumor growth mathematical model chosen to design the closed-loop control law, which is defined in Section 3, where the main results are provided. Simulations are reported in Section 4. Conclusions follow.

---

2. The tumor growth model

The model under investigation is a nonlinear model accounting for angiogenic stimulation and inhibition [14], and is given by the following Ordinary Differential Equations (ODE) system:

$\begin{array}{ccl} \dot x_1& = &-\lambda x_1 \ln\left(\frac{x_1}{x_2}\right), \\ \dot x_2& = &bx_1 - ( \mu + d x_1^{2/3})x_2 -c x_2 x_3, \end{array}$

(1)

where $x_1$, [mm$^3$], and $x_2$, [mm$^3$], denote the tumor volume and the carrying capacity of the vasculature, respectively. The first equation describes the phenomenology of tumor growth slowdown, as the tumor grows and resorts its available support with $\lambda$, [day$^{-1}$], denoting the tumor growth rate constant. In the second equation the first term ($bx_1$) represents the stimulatory capacity of the tumor upon the inducible vasculature, the second term accounts for spontaneous loss ($\mu x_2$) and for tumor-dependent endogenous inhibition ($dx_1^{2/3} x_2$) of previously generated vasculature, and the third term refers to the vasculature inhibitory action performed by an exogenous drug administration ($cx_2x_3$) with $x_3$, [mg/kg], denoting the serum level of the administered angiogenic inhibitor. As far as the parameters, $b$, [day$^{-1}$], is the vascular birth rate constant, $d$, [day$^{-1}$mm$^{-2}$], represents the endothelial cell death, $\mu$, [day$^{-1}$], is the rate constant of the spontaneous vascular inactivation, and $c$, [day$^{-1}$(mg/kg)$^{-1}$], is the sensitivity to the drug. According to the model literature [14], without loss of generality, parameter $\mu$ will be set equal to zero in the following, because it was found to be negligible (i.e., constitutive endothelial cell loss does not play a major role in this system).

Being the anti-angiogenic drug not directly administered in vein, a further compartment is considered to account for drug diffusion:

$x_3(t) = \int_0^t e^{-\eta(t-t')} u(t')dt',$

(2)

with $u$, [day$^{-1}$(mg/kg)], being the actual control law and $\eta$, [day$^{-1}$], being the degradation rate constant. As a matter of fact, the whole system (1)-(2) (with $\mu = 0$) may be written in a compact ODE form:

$\begin{array}{ccl} \dot x_1& = &-\lambda x_1 \ln\left(\frac{x_1}{x_2}\right), \\ \dot x_2& = &bx_1 - d x_1^{2/3}x_2 -c x_2 x_3, \\ \dot x_3& = &-\eta x_3 + u, \end{array}$

(3)

where $u$ stands for the anti-angiogenic drug administration rate, supposed to be delivered continuously.

We assume that the only available measurement, exploited to design the control law, is the size of the tumor, that is, the first component of the state vector $x$. These measurements are acquired at discrete sampling times. With this assumption, system (3) can be represented as

$\begin{align}\dot x = & f(x)+g(x)u, \end{align}$

(4)

$\begin{align} y_k = & h\big(x(kT)\big), \;\;\;\; k = 1, 2, \ldots, \end{align}$

(5)

where $T$ is the sampling time and $f:\mathbb{R}^3\rightarrow \mathbb{R}^3$, $g:\mathbb{R}^3\rightarrow \mathbb{R}^{3\times 1}$ and $h:\mathbb{R}^3\rightarrow \mathbb{R}$ are defined by

$\begin{align} f(x) = \begin{bmatrix}-\lambda x_1 \ln\left(\frac{x_1}{x_2}\right)\\ bx_1 - dx_1^{2/3}x_2 -c x_2 x_3\\ -\eta x_3 \end{bmatrix} &, \;\;\;\; g(x) = \begin{bmatrix}0\\ 0\\ 1\end{bmatrix}, \end{align}$

(6)

$\begin{align}h(x) = x_1 = C_b x &, \;\;\;\; C_b = [1\ 0\ 0]. \end{align}$

(7)

---

3. The control algorithm

Without a control input (i.e. with $u(t) = 0$), the open loop system described by (3) exhibits an asymptotically stable non-trivial equilibrium point, given by

$x_{1, ss} = x_{2, ss} = (b/d)^{\frac{3}{2}} = \rho.$

(8)

The goal of the proposed control scheme is to design a feedback control law $u$ that allows to reduce the tumor volume (state variable $x_1$) down to a desired level $r$, with $0 < r < \rho$. By requiring $\dot{x}_1 = \dot{x}_2 = \dot{x}_3 = 0$ in (3) it follows that the desired setting point $x_1 = r$ implies

$x_2 = r, x_3 = \frac{b-dr^{2/3}}{c}, u = \frac{\eta}{c}(b-dr^{2/3}).$

(9)

The aim of steering system (3) to the equilibrium (9) can be achieved by properly exploiting the output feedback linearization theory to synthesize a state-feedback control law and by using a state observer to provide an estimate of the state. Indeed, we do not have a complete knowledge of the state, but we assume to measure only the tumor size, i.e. the first component of $x_1$. Because of the discrete feature of the output, measurements are preliminary restated as continuous delayed outputs, and then an exponential observer for such systems is applied [2]. A proper post-processing of the available measurement is as well considered, in order to obtain a separation principle that allows to independently design the control and observer parameters. In summary, we define the following delayed output:

$\begin{align} {\bar{h}}(x) = & \ln(x_1), \label{eq:hbar} \end{align}$

(10)

$\begin{align}\overline{y}^\delta(t) = & {\bar{h}}\big(x(t-\delta(t))\big) = \ln \Big(x_1\big(t-\delta(t)\big)\Big) = \ln(y_k), \;\;\;\; t\in[kT, (k+1)T)\label{eq:delayed measurements} \end{align}$

(11)

where the last equality is obtained from $\delta(t) = t-kT$. The basic idea, originally formulated in [25], is that the sampled output can be considered as a delayed measurement with a non-continuous delay function that grows together with $t$ in $[kT, (k+1)T)$ and is reset to $0$ at each sampling point. In this way $t-\delta(t)$ is a constant time point $kT$ in the interval $[kT, (k+1)T)$ between pairs of measurements. It is worth noticing that system (4) endowed with the post-processing function ${\bar{h}}$ has full observation degree since

$\begin{align} \label{FRD} L_g{\bar{h}}(x) = L_gL_f{\bar{h}}(x) = 0, \end{align}$

(12)

$\begin{align}L_gL_f^2{\bar{h}}(x) =-\lambda c\neq 0, \end{align}$

(13)

where $L_f^j{\bar{h}}(x)$ stands for the Lie derivative of order $j$ of the scalar function ${\bar{h}}$ along the vector field $f:\mathbb{R}^n\rightarrow\mathbb{R}^n$, [15]

$L_f {\bar{h}} = \frac{d{\bar{h}}}{dx}f, \;\;\;\; L_f^j {\bar{h}} = \frac{dfL_f^{j-1}{\bar{h}}}{dx}f, \;\;\;\;j = 2, 3, \ldots.$

(14)

In order to define the observer-based control, we use a nonlinear change of coordinates based on the drift-observability map [6], which in our case is

$z = \Theta(x) = \begin{bmatrix} {\bar{h}}(x)-\ln(r) \\ L_f {\bar{h}}(x) \\ L^2_f {\bar{h}}(x) \end{bmatrix} = \begin{bmatrix} \ln(x_1)-\ln(r) \\ -\lambda \ln{\frac{x_1}{x_2}} \\ \lambda^2\ln{\frac{x_1}{x_2}} + b\lambda\frac{x_1}{x_2}-d\lambda x_1^{2/3}-c\lambda x_3 \end{bmatrix},$

(15)

and we consider the hypotheses needed by the observer based on this map.

$H_1$) The system is Globally Uniformly Lipschitz Drift-Observable (GULDO), i.e. $\Theta$ is a diffeomorphism in the domain of interest $D$ and the maps $\Theta$, $\Theta^{-1}$ are uniformly Lipschitz in $D$.

$H_2$) The function $L_f^3{\bar{h}}\big(\Theta^{-1}(z)\big)$ is uniformly Lipschitz in $D$.

$H_3$) The function $L_gL_f^2{\bar{h}}(x)$ is a constant.

$H_4$) The system has full observation degree.

In [2] it has been proven that, if $H_1$-$H_4$ are satisfied, then there exists a bound $\Delta$ for the time-varying delay and a positive $\beta_c>0$ such that there exist matrices $K_o$, $K_c$ and a scalar $\beta>\beta_c$ ensuring that the observer-based control law defined by

$\label{eq:output_control} u(t) = \frac{ L_f^3{\bar{h}}(\hat x)-K_c\hat z}{\lambda c}, \;\;\;\;\hat z = \Theta(\hat x),$

(16)

where $\hat x(t)$ is the observer estimate computed as

$\dot{\hat x} = f(\hat x) + g(\hat x)u(t) + J_{\Theta}(\hat x)^{-1}e^{-\beta \delta(t)}K_o\Big( \overline{y}^\delta(t) - {\bar{h}}\big(\hat x(t-\delta(t))\big)\Big), \label{eq:observer2}$

(17)

provides a closed-loop system with the origin of the state space $(z, \hat z)$ globally asymptotically stable for any initial condition of the system and of the observer estimate, and for any delay $\delta(t) < \Delta$. From (15) it is easy to check that $z =[0,0,0]^T$ is the desired setting point defined in (9).

Hypotheses $H_3$ and $H_4$ have been previously discussed and we have shown they are satisfied. Actually, the post-processing output (11) has been introduced to obtain $H_3$.

As for $H_1$-$H_2$, the natural domain of interest is

$\label{IniDomain} [0, \, \rho]\times[0, \, \rho]\times [0, \overline{x}_3],$

(18)

where $\overline{x}_3$ can be set equal to ${\bar{u}}/\eta$ with ${\bar{u}}$ an upper bound for the input $u$. Indeed, model (3) associated to the real system under consideration has a physical meaning only for $x_1, x_2>0$. Past literature has shown that system (3) has positive solutions for positive initial conditions and positive input function (see, e.g., [26]). Moreover, the model equations have been derived in order to exhibit, without a control action, a monotone increase of the pair $(x_1, x_2)$ to a plateau $(\rho, \rho)$ when the initial carrying capacity is larger than the tumor size, i.e. $0 < x_1(0) < x_2(0)$, [14]. By adding a positive input we introduce a decreasing term in the carrying capacity equation, thus preventing the higher uncontrolled steady-state to be reached. In summary, (18) is an invariant domain for the system solutions.

Unfortunately, (18) cannot be chosen as the domain of interest $D$. Indeed, to comply with $H_1$, consider the Jacobian of $\Theta$:

$J_{\Theta} = \frac{d\Theta}{dx} = \begin{bmatrix} \frac{1}{x_1}&0&0 \\ -\frac{\lambda}{x_1}&\frac{\lambda}{x_2}&0 \\ J_{31}&J_{32}&-c\lambda \end{bmatrix},$

(19)

where

$J_{31}(x) = \lambda \left( \frac{\lambda}{x_1} + \frac{b}{x_2} - \frac{2}{3}\frac{d}{x_1^{1/3}}\right), \;\;\;\; J_{32}(x) = - \frac{\lambda\left(b x_1 + \lambda x_2 \right)}{x_2^2}.$

(20)

The properties needed by $H_1$ require to restrict the domain of interest to the set

$D = [\varepsilon , \, \rho]\times[\varepsilon , \, \rho]\times [0, \overline{x}_3],$

(21)

where $\varepsilon$ is any arbitrarily small value satisfying $0 < \varepsilon < r < \rho$. Such a domain of interest $D$ contains the desired equilibrium point (9) and complies with $H_2$, since:

$L_f^3{\bar{h}}(x) = \lambda^2a(x)\log\left(\frac{x_1}{x_2} \right) - \lambda \left(\frac{bx_1+x_2}{x_2^2}\right)\xi(x)+c\eta\lambda x_3,$

(22)

with

$a(x) = \lambda+b\frac{x_1}{x_2}-\frac{2}{3}dx_1^{2/3},\;\;\;\; \xi(x) = bx_1 - d x_1^{\frac{2}{3}}x_2 -c x_2x_3,$

(23)

and $\Theta^{-1}$ is uniformly Lipschitz in $D$.

In order to make $D$ invariant for system (3), we have to ensure that all the system trajectories are contained in $D$ for any input. To this aim the following issues will be considered.

$h_1$) The time-horizon is limited. This point allows to weaken the invariant property of $D$ to ensure that $x(t)\in D$ only for $t\in [0, t_f]$, where $t_f$ is the time period in which the anti-angiogenic therapy is applied: it is not, clearly, a loss of generality from a clinical viewpoint.

$h_2$) The initial condition $x_1(0)>r$ is known whilst, regards to $x_2(0)$ we assume that $x_2(0)>x_1(0)$ holds. Also this point is not a loss of generality, since $x_1$ is the measured variable, and condition $x_2(0)>x_1(0)$ has been shown to exist according to a physically meaningful initialization.

These additional hypotheses allow to state the following proposition.

Proposition 1. Consider system (4) with initial conditions complying with $h_2$. Then, according to a physiologically meaningful choice of the model parameters, for any chosen time instant $t_f>0$ and for any positive $\varepsilon$ such that $(2/3)^3 < \varepsilon < r$, there can be found an upper bound ${\bar{u}}$ for the input $u(t)$, such that the domain

$D = [\varepsilon , \, \rho]\times[\varepsilon , \, \rho]\times \left[0, {\bar{u}}/\eta\right]$

(24)

is invariant in the interval $[0, t_f]$ for any initial condition in $D$, with $x_2(0)>x_1(0)$, and any input $u(t)\leq {\bar{u}}$.

Proof. We have already stated that the domain $[0, \, \rho]\times[0, \, \rho]\times \left[0, {\bar{u}}/\eta\right]$ is invariant for any positive initial condition and any positive input upper bounded by ${\bar{u}}$. The main step of the proof is to show that for any chosen time instant $t_f>0$ and for any positive $\varepsilon$ such that $(2/3)^3 < \varepsilon < r$, there can be found a constant input $u(t) = {\bar{u}}$, such that the evolution $X_{(X_0, {\bar{u}})}$ associated to the initial condition $X_0 = (x_1(0), x_2(0), 0)\in D$ with $x_2(0)>x_1(0)$ is lower bounded by $X_{(\widetilde X_0, {\bar{u}})} = (\tilde x_1(t), \tilde x_2(t), \tilde x_3(t))$, the evolution associated to the initial condition $\widetilde X_0 = (x_1(0), x_1(0), 0)\in D$. The continuity of the solution ensures that there exists a non-trivial interval $(0, t_0]$, with $0 < t_0\leq t_f$ such that $x_1(t)> \tilde x_1(t)$ and $x_2(t) > \tilde x_2(t)$ in $(0, t_0$]. Now, assume that $X_{(X_0, {\bar{u}})}$ is not lower bounded by $X_{(\widetilde X_0, {\bar{u}})}$ in $[0, t_f]$. Then, it would mean that there exists a time instant $\bar t \in (t_0, t_f]$ such that

ⅰ) $x_1(\bar t) = \tilde x_1(\bar t)$ and $x_2(\bar t)>\tilde x_2(\bar t)$ with $\dot x_1(\bar t)\leq \dot {\tilde x}_1(\bar t)$, or

ⅱ) $x_2(\bar t) = \tilde x_2(\bar t)$ and $x_1(\bar t)>\tilde x_1(\bar t)$ with $\dot x_2(\bar t)\leq \dot {\tilde x}_2(\bar t)$.

We will show that ${\bar{u}}$ can be set in order to make both these statements false. Indeed, assume ⅰ) is true. That means:

$-\lambda x_1(\bar t) \ln\left(\frac{x_1(\bar t)}{x_2(\bar t)}\right)\leq -\lambda \tilde x_1(\bar t) \ln\left(\frac{\tilde x_1(\bar t)}{\tilde x_2(\bar t)}\right)\;\;\;\; \Longrightarrow\;\;\;\; x_2(\bar t) \leq \tilde x_2(\bar t)$

(25)

which is a contradiction. Similarly, assume ⅱ) is true. That means:

$bx_1(\bar t) - d x_1^{2/3}(\bar t)x_2(\bar t) -c x_2(\bar t) x_3(\bar t) \leq b\tilde x_1(\bar t) - d \tilde x_1^{2/3}(\bar t)\tilde x_2(\bar t) -c \tilde x_2(\bar t) \tilde x_3(\bar t)$

(26)

and, after straightforward simplifications:

$\label{ineq} \frac{x_1(\bar t) - \tilde x_1(\bar t)}{ x_1^{2/3}(\bar t) - \tilde x_1^{2/3}(\bar t)} \leq {\frac{d\tilde x_2(\bar t)}{b}}.$

(27)

Because of the sub-linear growth of $x^{2/3}$, it comes that the left-hand-side of (27) is greater than 1 whenever $\tilde x_1>(2/3)^3$ which is rather a small value if compared to the usual amounts related to $x_1$ and $x_2$, about 3 orders of magnitude higher [14]. Therefore, according to the constraint on $x_1(0)>\varepsilon > (2/3)^3$, it is reasonable to find a value for $\bar u$ that ensures $\tilde x_1>(2/3)^3$. On the other hand, a proper choice of ${\bar{u}}$ can set the right-hand-side of (27) small enough to make inequality (27) false. Indeed, according to a physiologically meaningful choice of the model parameters suggested by Hanhfeldt [14], we have that the ratio $b/d\simeq 0.0015$ (see also Table 1). Again, it is reasonable to find (e.g. by simulation) a value for $\bar u$ that ensures an upper bound for $\tilde x_2$ such that $d\tilde x_2(\bar t)/b < 1$: indeed, the usual target for $r$ is smaller than $1/0.0015\simeq 667$. Therefore, inequality (27) is false since, according to the choice of ${\bar{u}}$ we have

$\frac{x_1(\bar t) - \tilde x_1(\bar t)}{ x_1^{2/3}(\bar t) - \tilde x_1^{2/3}(\bar t)}>1 \;\;\;\;{\rm and} \;\;\;\;\frac{d\tilde x_2(\bar t)}{b} < 1.$

(28)

In summary, we can tune ${\bar{u}}$ so that, for fixed inputs $u(t) = {\bar{u}}$, the solution $X_{(X_0, {\bar{u}})}$ is lower bounded by $X_{(\widetilde X_0, {\bar{u}})}$. By allowing any other smaller input $u(t)\leq {\bar{u}}$ we have that the solution $X_{(X_0, u(t))}$ would still be lower bounded by $X_{(\widetilde X_0, {\bar{u}})}$. Indeed, we would have smaller values of $x_3(t)$ and, by consequence, higher values of $x_1(t)$ and $x_2(t)$.

This proves the invariance of $D$ for any chosen $\varepsilon >(2/3)^3$ provided that $u(t)\leq {\bar{u}}$.

Due to the presence of input saturation the regulation theorem provided by [2] cannot be straightforwardly applied. Instead, the following result is provided.

Theorem 3.1. Consider the closed-loop system (4), (11), (16)-(17). Consider any positive $\varepsilon$ such that $(2/3)^3 < \varepsilon < r$. Then, there exists a positive value ${\bar{u}}$ such that, under the hypothesis that the feedback control law (16) complies with the bound $u(t)\leq {\bar{u}}$, then, there exists a maximum value of the sampling interval $T$ such that

(ⅰ) it is possible to design the observer parameters $(K_o, \beta)$ in (17) in order to ensure exponential convergence to zero of the observer error when $t_f\rightarrow \infty$;

(ⅱ) it is possible to design the control gain $K_c$ in order to ensure asymptotic convergence to zero of $z(t)$ for $t_f\rightarrow\infty$.

Proof. According to Proposition 1 there exists an upper bound ${\bar{u}}$ such that the domain $D$ defined in (24) is invariant for the system solution, provided that $u(t)\leq {\bar{u}}$. Moreover, as previously stated $H_1$-$H_4$ are satisfied in $D$. According to the drift-observability map, it is therefore possible to write system (4) endowed with output $\overline{y}^\delta(t)$ in the $z$ coordinates,

$\begin{align}\label{eq:zmodel} \dot{z}(t) = & A_b z(t) + B_b\left( (L_f^3{\bar{h}}\left(\Theta^{-1}(z(t))\right)-c\lambda u(t)\right), \end{align}$

(29)

$\begin{align} \overline{y}^\delta(t) = & C_b z(t-\delta(t)), \end{align}$

(30)

where $A_b$, $B_b$ are the Brunowski matrices

$\label{ABBrun} A_b = \begin{bmatrix} 0&1&0\\ 0&0&1\\ 0&0&0\end{bmatrix}, \;\;\;\; B_b = \begin{bmatrix}0\\ 0\\ 1\end{bmatrix}.$

(31)

Let us define the error $e_z = z-\hat z$. By substituting the control input in the error dynamics we obtain

$\begin{array}{ccl} \dot e_z(t) & = & A_b e_z(t) + B_b\Big(L_f^3{\bar{h}}\big(\Theta^{-1}(z(t))\big)- L_f^3{\bar{h}}\big(\Theta^{-1}(\hat z(t))\big)\Big) \\ & &- e^{-\beta \delta(t)}K_oC_be_z\big(t-\delta(t)\big). \end{array}$

(32)

Denoting $\gamma_L$ the Lipschitz coefficient of $L_f^3{\bar{h}}\big(\Theta^{-1}(\cdot)\big)$, which is bounded by $H_2$,

$\left\| L_f^3{\bar{h}}\Big(\Theta^{-1}\big(z(t)\big)\Big)- L_f^3{\bar{h}}\Big(\Theta^{-1}\big(\hat z(t)\big)\Big)\right\| \leq \gamma_L \left\| z(t)-\hat z(t)\Vert = \gamma_L \Vert e_z(t)\right\|.$

(33)

Eq. (32) does not contain $K_c$, and the convergence to $0$ of the observer error $e_z(t)$ is independent from the choice of the control gain. Indeed, by choosing $K_o$ to assign eigenvalues with arbitrary negative real parts to $A_b-K_oC_b$ the existence of $T\geq\delta(t)$ that ensures the exponential stability of $e_z(t)$ has been proved in [2]. The bound on $T$ depends on $\beta$, $\gamma_L$ and the norm of $K_o$. As a consequence, for any $K_o$ that makes stable $A_b-K_oC_b$ there is a maximum delay value $T$, that depends also on $\beta$, for which (32) is exponentially stable.

In closed loop, assuming $u(t)\leq{\bar{u}}$ the equation for $z(t)$ is obtained by replacing (16) in (29)

$\dot z = A_b z+B_b K_c z -B_bK_ce_z + B_b\Big(L_f^3{\bar{h}}\big(\Theta^{-1}(z)\big)- L_f^3{\bar{h}}\big(\Theta^{-1}(z-e_z)\big)\Big).$

(34)

If $K_o$ and $\beta$ have been designed in order to ensure the exponential convergence to $0$ of $e_z$, the solution of (34) converges exponentially to the solution of

$\dot \zeta = \left(A_b +B_b K_c\right)\zeta, \;\;\;\;\zeta(0) = z(0).$

(35)

Hence, by choosing $K_c$ that assigns eigenvalues of $A_b+B_bK_c$ with negative real parts it follows that $\zeta$ and $z$ are exponentially stable [5].

Corollary 1. Hypothesis $H_1$ guarantees that $\Theta^{-1}$ is uniformly Lipschitz and, consequently, the exponential convergence to $0$ of $z(t)$ implies the exponential convergence of $x(t)$ to the desired setting point.

---

4. Simulation results

The aim of the simulations reported in this section is to investigate the performance of the proposed algorithm on the ground of the experiments carried out in [14] on mice from which the following parameter values, constraints and initial conditions have been considered.

The desired level for the tumor size is fixed at $r = 135$ mm$^3$ with the end of treatment at $t_{f} = 13$ days. Taking into account the results in [14], the simulations have been performed by fixing the model parameters to the values reported in the Table 1, which are related to a Lewis lung carcinoma. In particular, in [14] the authors fitted the growth parameters ($\lambda$, $b$, $d$) from control data of tumor volume growth, meanwhile, the vascular inactivation rate $c$ and the drug clearance inactivation rate $\eta$ were estimated from data of tumor treated with endostatin and using the fitted Gompertz parameters.

The value of the drug sensitivity $c$ has been fixed to a maximal tolerable daily level of 20 mg/kg/day. Initial conditions are set by $x_1(0) = 200$ mm$^3$ and $x_2(0) = 625$ mm$^3$.

In the previous sections, it has been proved that the control scheme allows to set independently the control and the observer parameters, since the structure of the tumor growth model guarantees the separability of estimation and feedback control algorithms. This separation principle is followed in the setting of the observer-based control parameters.

In particular, we set the observer parameters $(K_0, \beta)$ in (17) as

$K_o = [1.20, \, 0.44, \, 0.05]^{T}, \;\;\;\;\beta = 1.$

(36)

$K_o$ assigns the eigenvalues $[-0.2, \, -0.4, \, -0.6]$ to the matrix $A_b - K_oC_b$. The value of $\beta$ has been empirically chosen in order to have the convergence of the observer for the desired sampling values and gain $K_o$. In the observer initialization we chose $\hat{x}_1(0) = x_1(0) = 200$ mm$^3$, as we can reasonably assume to know the actual tumor volume at the beginning of the simulation, and $\hat{x}_3(0) = x_3(0) = 0$, because there is no drug concentration at $t = 0$. Conversely, the only a priori knowledge on the initial value of the carrying capacity is that it must be higher than the initial tumor volume, but we ignore its actual value. Simulations have been carried out for several initial estimates of the kind $\hat x_2(0) = \alpha x_2(0)$ with $\alpha \in [0.6, 2]$.

The feedback control gain in (16) is $K_c = [20.25, \, 24.75, \, 9.00]$, that assigns the eigenvalues $[-1.5, \, -3.0, \, -4.5]$ to the matrix $A_b + B_bK_c$.

In order to evaluate the performance of the algorithm we define the following goals.

1. $x_1(13) \in [0.9r, 1.1r]$,

2. $m = \int_0^{13}u(t)dt/13 \leq 20$ (mg/kg/day),

where the first goal ensures that the tumor size at the end of the treatment (day $13$) is close to the desired value $r$, within $10\%$ of tolerance, and the second goal requires that the average drug administration ($m$) is maintained below $20$ mg/kg/day. Notice that, from a clinical viewpoint, the item concerning the tumor growth at the end of the treatment could be relaxed to have only an upper bound. The requirement for $x_1(13)$ to be close to the reference level within $10\%$ aims at evaluating the performance of the control law, independently of clinical requirements.

Different simulation campaigns have been carried out according to different sampling times, which have been set at $1$, $2$ and $3$ days, that are the measurement intervals used for the experiments reported in [14]. An example of the system dynamic is reported in Fig. 1, with the model parameters provided by Table 1 and the sampling time fixed at 3 days. Notice that the observer converges within less than $2$ day.

In order to investigate the robustness of the algorithm efficiency for different model parameter settings, we simulated the tumor growth for a "population" of $1000$ different mice, each controlled by the same regulator keeping fixed both the observer parameters $(K_{o}, \beta)$, and the control gain $K_{c}$. In particular, each mouse from the population is associated to:

$\bullet$ the initial condition $x(0)$ uniformly sampled in the interval $[0.4, \, 1.4]*x^{N}(0)$, where $x^{N}(0) = (200, 625, 0)$ is the nominal initial condition given above;

$\bullet$ the initial value of the carrying capacity for the observer, $\hat x_2(0)$, uniformly sampled in the interval $[0.6, \, 2]*x_2(0)$;

$\bullet$ a set of model parameters, with each parameter (namely $p_j$) uniformly sampled in the interval $\Delta_{p, j} = [0.5, \, 1.5]*p^{N}_{j}$, where $p^{N}_{j}$ represent the nominal values of $p_{j}$ reported in Table 1.

Figure 2 shows the percentage of successes for the two proposed criteria for the population. It is apparent that the performances of the closed-loop control do not vary for sampling intervals in $[1, 3]$ days. It is worthwhile to stress that, even for the less-than-$20\%$ cases where the final $x_1$ value exceeds the bound $x_{1}(13)>149 = 1.1r$ mm$^{3}$, we have a strong reduction of the tumor growth, which is at most equal to 167 mm$^{3}$ (the worst case).

---

5. Conclusions

Based on a mathematical model of tumor growth, this work proposes a closed-loop control law aiming at reducing the tumor volume. The model accounts for angiogenic stimulation and inhibition, and is one of the most adopted in the literature to simulate and predict the effects of anti-angiogenic drug delivery. Based on the feedback linearization theory, the control law makes use of an observer for nonlinear systems in order to design the model-based control by means of only available measurements. Theoretical results ensure that the control gain of the regulator can be set independently of the observer gain, thanks to the structural properties of the tumor growth model. Numerical simulations show the effectiveness of the control law in spite of a wide range of variation of the (not measured) carrying capacity as well as a noticeable level of robustness with the respect the uncertainties affecting the model parameters.

---