Nonlinear Adaptive Optimal Controller Design for Anti-Angiogenic Tumor Treatment

Abstract

Angiogenesis is an important process in tumor growth as it represents the regime when the tumor recruits blood vessels from the surrounding tissue to support further tumor growth. Anti-angiogenic treatments aim to shrink the tumor by interrupting the vascularization of the tumor; however, the anti-angiogenic agents are costly and the tumor response to these agents is nonlinear. Simple dosing schemes, e.g., a constant dose, may yield higher cost or lower efficacy than an approach that considers the tumor system dynamics. Hence, in this study, the administration of anti-angiogenic treatment is considered as a nonlinear control problem. The main aim of the controller design is to optimize the anti-angiogenic tumor therapy, specifically, to minimize the tumor volume and drug dose. Toward this aim, two nonlinear optimal controllers are presented. The first controller ensures exponential tracking of a desired, optimal tumor volume profile under the assumption that all parameters in the system model are known. The second controller, on the other hand, assumes all the parameters are unknown and provides asymptotic tracking. Both controllers take pharmacokinetics and pharmacodynamics into account, as well as the carrying capacity of the vascular network. Lyapunov based arguments are used to design the controllers, using stability arguments, and numerical simulation results are presented to demonstrate the effectiveness of the proposed method.

1. Introduction

Breast cancer has the highest incidence rate compared to other types of cancer [1]. An overview of the tumor angiogenesis is provided by the US National Cancer Institute [2] and summarized in Figure 1. Initially, tumor cells get nutrients from the surrounding tissues via diffusion. The diffusion supported growth regime ends when physical limits (1–2 mm of tumor diameter) prevent nutrients and oxygen from reaching the core of the tumor. Then, an angiogenesis triggering event occurs; the tumor cells release factors, such as vascular endothelial growth factor (VEGF) or basic fibroblast growth factor (bFGF) [3], that instigate the recruitment of blood vessels to support further tumor growth. The tumor and support vasculature can then co-develop as the blood vessels supply nutrients and remove waste products from the proliferating tumor cells, and the growing tumor signals additional vascular growth, primarily vascular endothelial cells. The next major event is metastasis when the tumor spreads to the other tissues and organs.

Anti-angiogenic treatment was proposed by Folkman [4] and aimed to control the growth of the vascular network, described above, to indirectly control the tumor. The relationship between the growth of the vascular network and the growth of the tumor is dynamic, complex, and nonlinear. Therefore, applying a fixed or a time-varying dose of anti-angiogenic agents, without observing the response of the tumor, i.e., open-loop control, may fail to shrink the tumor. Closed-loop control may be more effective to control the growth of the vascular network, so as to shrink the tumor, by using feedback signals from the co-developing tumor and vascular network. Since the commercially available anti-angiogenic agents, like angiostatin and endostatin, are expensive, an optimal control approach can be employed to keep the tumor size at a certain (desired) value by using an optimum dose of an anti-angiogenic agent. Ledzewicz and her coworkers have proposed many open-loop [5,6] and closed-loop optimal control schemes [7,8]. Swierniak et al. also presented an optimal controller addressing the same problem in [9]. Nath et al. proposed a novel approach that combines closed-loop tracking control with an optimization to minimize the anti-angiogenic agents that needed to move the system to a setpoint (tumor size) [10]. In 2013, Hasirci et al. proposed enhancements to the core tracking controller: a robust controller in [11] and an adaptive controller in [12]. Researchers have used nonlinear adaptive and robust control techniques in a wide range of applications, success in controlling similar dynamical systems [13,14,15,16,17] inspires the current work.

As a general control design methodology, the salient features of the system to be controlled are first approximated with a mathematical model. In biological systems, especially, it is a good and usual practice to make simplifying assumptions to make the mathematical model tractable, while balancing accuracy and resolution. As experience and evidence build, additional complexities can be added to the model. The controllers reviewed above are based on a simplified dynamic model where tumor volume and vascular carrying capacity are the only state variables. This model assumes drug concentration is equal to drug dose, and also assumes that the effect of the anti-angiogenic agents is instantaneous. A more realistic and, as a result, more complex view of the angiogenesis process dictates the incorporation of pharmacokinetics (PK) and pharmacodynamics (PD) into the system model. As clearly stated in [18], pharmacokinetics equations describe the concentration of a drug and pharmacodynamics equations model the effectiveness of drug dose. Incorporating PK/PD equations to the angiogenic tumor growth model begins to add the clinic reality of drug delivery to the model.

In this paper, we present a novel approach for real-time dosing of the anti-angiogenic treatment by considering the richer tumor growth model, which incorporates PK/PD. The main aim of the control design is to optimize the anti-angiogenic tumor therapy for minimizing the tumor volume by using an optimal amount of the anti-angiogenic agents. Two different controllers are presented to achieve the control objective. To minimize the tumor volume, we first formulate a performance index to be minimized and then a nonlinear optimal exact model knowledge controller is designed to drive the tumor volume to an optimal desired time-varying trajectory. Then, a more realistic control design approach is also presented for the same control objective; a nonlinear optimal adaptive controller assuming all the system parameters is unknown. For the adaptive controller, a least-squares estimation technique is used to identify the system parameters in the performance index and the desired optimal trajectory is generated by using an optimization algorithm, which seeks the minimum of the performance index. As demonstrated by simulation results, this optimization algorithm successfully finds the optimum values of the carrying capacity of vasculature and the tumor volume. The paper is organized as follows: Section 2 contains discussion about the modeling issues and a definition of the control problem. Section 3 contains the details of the design of an exact model knowledge nonlinear optimal controller, and then an adaptive controller. Performance results of the proposed adaptive scheme are obtained via numerical simulation in Section 4. The last section discusses the results and highlights some conclusions.

2. System Model and Control Problem Definition

Early mathematical models for tumor growth comprised one state variable: tumor volume (p(t)) [19]. The carrying capacity of the vasculature network (q(t)) concept was proposed to account for co-growth of the tumor and vasculature, i.e., angiogenesis, and then a killing factor was introduced to model anti-angiogenic therapy, which targets endothelial cells instead of the tumor cells [20]. Many mathematical models have been proposed for this two-state system, some of which attempt to fully describe the complexity of biological process [21,22,23,24,25,26,27]. Since these models are not tractable for mathematical analysis [28], Hahnfeldt et al. proposed and biologically validated a simpler two-state model to describe the interaction between the carrying capacity of vasculature and tumor size [29]. A useful modification of this model was presented by Ergun et al. [30], and more recently by d’Onofrio and Gandolfi [31]. Incorporation of PD/PK equations to this modified model has been presented by Ledzewicz et al. [32] and is given by:

$$\begin{array}{l}\dot{p}=\alpha p\left(1-\frac{p}{q}\right)\\ \dot{q}=bq-d{p}^{2/3}q-Gscq\\ \dot{c}=-mc+hu ; c\left(t_0\right)=0\end{array}$$

(1)

where p(t) is the tumor volume in mm³, q(t) is the carrying capacity of the endothelial cells, also measured in mm³, c(t) is the plasma concentration of the exogenous angiogenic inhibitors, applied as treatment and measured in [conc.], and α is the tumor growth parameter. The term bq accounts for the proliferation kinetics of the endothelial cells and the term dp^2/3q models endogenous inhibition of the tumor. Because the inhibitors release through the surface of the tumor mass, the 2/3 exponent represents the conversion of the tumor volume into a tumor surface area, and the product p^2/3q is due to the interactions of endogenous inhibitors with endothelial cells. The parameters b and d are constants. The constant G denotes the anti-angiogenic killing parameter, and u(t) is the manipulated control input which corresponds to the anti-angiogenic dose. The drug dose u(t) and concentration of c(t) of the inhibitors are linked by a first-order, linear, ordinary differential equation where m and h are constant parameters [32]. The effect of the applied drug is proportional to the concentration of the inhibitor, given as effect = sc where s $\in$ [0,1]. Note that “[conc.]” is generic concentration and that the units of concentration do not affect the states q(t) and p(t) when the parameters G, m, and h are written for a specific unit of concentration.

Assumption 1.

Only the biologically realistic domain where all the state variables are always greater than zero for all time instants, i.e., p(t) > 0, q(t) > 0, and c(t) > 0, will be considered. It is assumed that the tumor volume p(t), the carrying capacity of the vascular network q(t), and the concentration of the angiogenic inhibitors c(t) are measurable.

The control objective is to minimize the tumor volume, p(t), by using an optimum dose of angiogenic inhibitors, u(t). This objective is met if the following performance index is minimized:

$$J=p+{\left(\frac{Su}{d}\right)}^{3/2}$$

(2)

where S = Gs. The performance index J(t) $\in ℝ$ captures the treatment goal by using the summation of tumor size and drug dose. To obtain the same units, mm³, for each term in (2), the drug dose is multiplied by a factor of S/d and then raised to 3/2. The steady-state expression of the performance index is:

$$J_0=p_0+{\left(\frac{S{u}_0}{d}\right)}^{3/2}$$

(3)

where J₀, p₀, and u₀ are the steady-state values of J, p, and u, respectively. The minimum of this performance index gives the minimum tumor volume that can be obtained by using the minimum drug dose; these values are found from the steady-state properties of the system given in (1) by calculating its equilibria corresponding to a constant drug dose, u₀. By setting left-hand sides of the state equations in (1) to zero, an equilibrium at p₀ = q₀ = 0 is obtained, which is not an admissible point [26], and:

$$p_0=q_0={\left(\frac{b-S{c}_0}{d}\right)}^{3/2}, c_0=\frac{h}{m}{u}_0.$$

(4)

This equilibrium point analysis reveals some important properties of the system given in (1). First, it is clear that the tumor volume is equal to the carrying capacity at steady-state. It is also clear that an uninterrupted or a long-term therapy is needed to prevent the growth of the tumor. By selecting the value of the parameter related to PD as s = 0.8, and by setting the PK parameters m = h = 1, which is a case often considered in the literature [33], and finally, by using the system parameter values taken from [31] and given in

Table 1. Parameter Values [31].

Parameter	Value	Unit
α	1.08	[day]−1
d	3.63 × 10−4	[day]−1 [mm]−2
b	0.243	[day]−1
G	1.3	[day]−1 [conc.]−1

, respective equilibria of the tumor volume and the carrying capacity are found as p₀ = q₀ = 17,320 mm³, if the control input u(t) is equal to zero. Figure 2 shows the variation of J₀ with respect to the steady-state values of u(t), denoted u₀, and Figure 3 shows the variation of J₀ with respect to the steady-state values p(t) and q(t), denoted by p₀ and q₀ for selected parameter values.

The minimum of the curves in Figure 2 and Figure 3, i.e., the optimum values of the variables, are found to be:

$$\begin{array}{l}J*=\mathrm{12,247}{ \left[\mathrm{mm}\right]}^3\\ p*=q*=6118{ \left[\mathrm{mm}\right]}^3\\ u*=0.1169 [\mathrm{conc}.].\end{array}$$

(5)

Considering the results of the analysis given above, the control objective can be more precisely defined as driving q(t) to its optimum value, q*. Since the tumor volume is equal to the carrying capacity of the vascular network at steady-state, this necessarily means that the tumor volume, p(t), will also be driven to q* = p* by using the minimum drug dose, if the control objective is achieved.

3. Control Design

In this section, two nonlinear optimal backstepping controllers are presented. The first controller encapsulates an exact-model knowledge controller that ensures exponential tracking of a desired optimal tumor volume under the assumption that all parameters in the system model are known. The second controller, on the other hand, is an adaptive controller that supports the practical reality that the parameters are unknown and provides asymptotic tracking.

3.1. Exact Model Knowledge Controller

Assuming all parameters in the tumor model are known, a complete control scheme is shown in Figure 3. The tracking controller uses feedback to drive the tumor volume to a desired, optimal point, q_d. That optimal point, q_d, is known a priori from the static optimization summarized in (5). A tracking controller shown in Figure 4 is needed that moves the tumor system towards this optimum value.

To quantify the performance of the exact model tracking knowledge controller to be designed, an error signal can be defined as:

$$e_1=q-q*.$$

(6)

Differentiating and substituting from (1) yields the error dynamics:

$${\dot{e}}_1=bq-d{p}^{2/3}q-Gscq.$$

(7)

This equation lacks a control input to be manipulated to stabilize the error variable; the actual control input appears in the last state equation given in (1). In such cases of coupled subsystems, the backstepping technique [34] provides a tool to design a control input signal that is propagated through the coupled subsystems. To apply this procedure, a virtual control input, c_d(t), is added and subtracted to the right-hand side of (7) as:

$${\dot{e}}_1=bq-d{p}^{2/3}q-Scq-S{c}_{d}q+S{c}_{d}q.$$

(8)

In practice, c_d(t) is the desired trajectory for the concentration of the inhibitor. A new error variable, e₂(t), is defined to quantify the difference between the actual concentration of the inhibitor and the desired concentration of the inhibitor as:

$$e_2=c-c_d$$

(9)

Temporarily ignoring the last two terms in (9), a desired trajectory for the concentration of the inhibitors, that would force e₁(t) towards zero, can be designed as:

$$c_d=\frac{1}{Sq}\left(bq-d{p}^{2/3}q+k_1{e}_1\right)$$

(10)

where k₁ is a positive control gain. Substituting the proposed trajectory yields the closed-loop dynamics for e₁(t) as:

$${\dot{e}}_1=-k_1{e}_1-Sq{e}_2.$$

(11)

Investigating the e₂(t) dynamics yields:

$$\begin{array}{c}{\dot{e}}_2=\dot{c}-{\dot{c}}_d\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}=-mc+hu-\frac{d}{dt}\left\{\frac{1}{Sq}\left(bq-d{p}^{2/3}q+k_1{e}_1\right)\right\}\\ \hspace{1em}\hspace{1em}=W\theta +hu\end{array}$$

(12)

where

$$W=\left[-c\hspace{1em}\frac{2}{3}{p}^{2/3}\left(1-\frac{p}{q}\right)\hspace{1em}{k}_{1}q{e}_2\hspace{1em}{k}_1^2{e}_1\right]$$

(13)

is a 1-by-4 regression vector containing measurable variables and known gains, and

$$\theta ={\left[\begin{array}{cccc}m& \frac{\alpha d}{S}& S& 1\end{array}\right]}^T$$

(14)

is a 4-by-1 vector containing the parameters. The actual control input, the concentration of angiogenic inhibitor, can be designed as:

$$u=\frac{1}{h}\left(-k_2{e}_2-W\theta +Sq{e}_1\right)$$

(15)

where k₂ is a positive control gain, which renders:

$${\dot{e}}_2=-k_2{e}_2+Sq{e}_1$$

(16)

Remark 1.

The above approach means that q(t) will be controlled directly, and that p(t) is indirectly controlled in the sense that it follows q(t) towards an equilibrium. From (1) and Assumption 1, it is evident that $\dot{p}(t)\le 0$ when $p(t)\ge q(t)$, thus $p(t)$ decreases. If $p(t)<q(t)$, then $p(t)$ will increase until $p(t)=q(t)$. It will start to decrease again if $p(t)\ge q(t)$. Therefore, $p(t)$ is bounded as long as $q(t)$ is bounded.

Theorem 1.

For the tumor growth model given in (1), the dose of the angiogenic inhibitor applied according to (15) ensures that the concentration of the inhibitor exponentially follows its trajectory, c_d(t), in order for the carrying capacity, q(t), of the vascular network to reach the constant value, q*, exponentially fast. In other words, e₁ goes to zero exponentially fast. It also guarantees that all signals in the closed-loop system are bounded.

Proof 1.

Consider the following nonnegative Lyapunov function:

$$V_1=\frac{1}{2}{e}_1^2+\frac{1}{2}{e}_2^2.$$

(17)

Taking the time derivative and then substituting the closed-loop dynamics of e₁ and e₂, given by (11) and (16), respectively, yields:

$${\dot{V}}_1=-k_1{e}_1^2-k_2{e}_2^2.$$

(18)

By considering (17) and (18) together, one can write:

$$V_1=V_1(t_0)\exp \left\{-\eta t\right\},\eta \in {ℝ}_{+}$$

(19)

where V₁(t₀) is the initial condition of V₁(t) and exp{·} is the base of the natural logarithm. This means that V₁(t) and its components, e₁(t) and e₂(t), go to zero exponentially. According to (6), if e₁(t) goes to zero, then the carrying capacity, q(t), is bounded, i.e., $q(t)\in {\mathcal{L}}_{\infty }.$ From Remark 1, $p(t)\in {\mathcal{L}}_{\infty }$; therefore, $c_d(t)\in {\mathcal{L}}_{\infty }.$ From (9), it can be concluded, $c(t)\in {\mathcal{L}}_{\infty }.$ Finally, from (15), $u(t)\in {\mathcal{L}}_{\infty }.$ □

The following simulation is presented to illustrate the performance and feasibility of the proposed controller applied to the tumor volume problem defined in this section. Note that the approach used to design this exact model controller will be expanded to account for parameter uncertainty in the following section. The simulation parameter values were provided in Section 2. Figure 5 shows the time variation of the error between the actual and optimum carrying capacity, defined in (6)—it goes to zero exponentially, as claimed. Practically, this means that the actual value of the carrying capacity of the vascular network goes to its optimum value, 6118 mm³, if the dose of the angiogenic inhibitors is applied to the tumor, as in Figure 6. This figure also demonstrates that the required dose of the anti-angiogenic agents is bounded. Finally, Figure 7 illustrates that the tumor volume follows the carrying capacity to reach its optimum value at steady-state. Values of the control gains were selected as k₁ = 0.1 and k₂ = 10 for this simulation. Initial values were assigned as p(t₀) = 12,000 mm³ and q(t₀) = 14,000 mm³. This approach is now extended to account for parameter uncertainty in the vector in (14).

3.2. Adaptive Controller

Determination of the “exact” values of the system parameters b, d, and S, which appear in the carrying capacity dynamics, for a specific patient is a very difficult and time-consuming task. The biological experiments that would be required to measure the exact values of the parameters for each patient would be costly, may delay the start of treatment, and would likely ultimately fail because of the uncertainty in measuring such a complex biological system. A more realistic and clinically feasible approach is to achieve the control objective, to minimize the tumor volume, p(t), by using an optimum dose of angiogenic inhibitors, u(t), despite measurement uncertainties in tumor growth model parameters. To end this, an estimate of the performance index given in (2) is defined, denoted by $\widehat{J}(t)\in ℝ$, as:

$$\widehat{J}=p+{\left(\frac{\widehat{S}c}{\widehat{d}}\right)}^{3/2}$$

(20)

where $\widehat{S}(t),\widehat{d}(t)\in ℝ$ are the estimates of the parameters S and d, respectively. Such a definition reveals an important issue with the control problem definition; the control problem will no longer be a set point problem, i.e., the desired trajectory for the carrying capacity will not be a fixed point, q*. Instead, the desired carrying capacity will be a time-varying trajectory as knowledge about the model parameters evolves.

Figure 8 shows an outline of the solution and indicates the three key components. First, the calculation of the performance index given in (20) requires the individually estimated values of the system parameters; a convergent estimation scheme is proposed to find the estimates of these parameters. Specifically, a least-squares estimation technique is adopted to generate the estimates of the unknown constant parameters. Second, an adaptive controller is proposed so that the carrying capacity of the vascular network follows a time-varying desired optimum trajectory $q_d(t)$, that is $q(t)\to q_d(t)$ as $t\to \infty$. Third, the desired trajectory, $q_d(t)$, is generated dynamically and online, using a numerical optimization method to minimize the performance index given by (20), such that $q_d(t)\to q*$ where $q*$ is the optimum value of $q(t)$ at steady-state.

3.2.1. Least-Squares Parameter Estimator

The least-squares estimation technique is used in many types of parameter estimation schemes. The main idea behind the technique is to extract the maximum information about the unknown parameters by using the measurable quantities in the system. Since the tumor volume, carrying capacity, and the concentration of the inhibitors are assumed to be measurable, information about the unknown parameters can be extracted as explained in the following.

To facilitate the estimator development, the q-dynamics given in (1) are parameterized as:

$$Q=W_e{\theta }_e$$

(21)

where Q denotes $\dot{q}(t)$,

$$W_e(t)\triangleq \left[q\hspace{1em}-p^{2/3}q\hspace{1em}-cq\right]\in {ℝ}^{1\times 3}$$

(22)

is a measurable regression vector, and

$${\theta }_e={\left[\begin{array}{ccc}b& d& S\end{array}\right]}^T\in {ℝ}^3$$

(23)

is a vector of unknown parameters. Such a parameterization allows the design of an adaptive estimation rule for the unknown parameters, given in (23), through the known quantities found in (22). To further facilitate the estimator design, a prediction error $\epsilon (t)\in ℝ$ is defined as:

$$\epsilon \triangleq Q_f-{\widehat{Q}}_f$$

(24)

in which $Q_f(t)\in ℝ$ is a filtered signal defined in the form of:

$${\dot{Q}}_f\triangleq -\beta Q_f+\beta Q ; Q_f(t_0)=0$$

(25)

where $\beta \in ℝ$ is a positive constant. It should be noted that (25) cannot be implemented directly since Q(t) is not measurable; however, an implementable form of this signal is provided in Appendix A. ${\widehat{Q}}_f(t)\in ℝ$ is the estimate of Q_f(t) and is defined as:

$${\widehat{Q}}_f\triangleq W_f{\widehat{\theta }}_e$$

(26)

where $W_f(t)\in {ℝ}^{1\times 3}$ is a filtered regression vector expressed as:

$${\dot{W}}_f\triangleq -\beta W_f+\beta W_e ; W_f\left(t_0\right)=0_{1\times 3}$$

(27)

and 0_1×3 denotes a 1-by-3 vector of zeros. In (26), ${\widehat{\theta }}_e={\left[\begin{array}{ccc}\widehat{b}& \widehat{d}& \widehat{S}\end{array}\right]}^T\in {ℝ}^3$ is the estimate vector of unknown parameters. Substituting (21) into (25) yields:

$${\dot{Q}}_f+\beta Q_f=\beta W_e{\theta }_e.$$

(28)

The last expression can be rewritten as:

$${\dot{Q}}_f+\beta Q_f={\dot{W}}_f{\theta }_e+\beta W_f{\theta }_e$$

(29)

where (27) was used. After taking the time derivative of (26), and then adding and subtracting the term ${\dot{W}}_f{\widehat{\theta }}_e$ to the right-hand side of the resulted expression, the following expression will be obtained:

$${\dot{\widehat{Q}}}_f+\beta {\widehat{Q}}_f=\frac{d}{dt}\left(W_f{\widehat{\theta }}_e\right)+\beta W_f{\widehat{\theta }}_e$$

(30)

where (26) and (27) were used. After subtracting (30) from (29), and utilizing (24) and (27), the resulting expression can be written as:

$$\dot{\epsilon }+\beta \epsilon =\frac{d}{dt}\left(W_f{\tilde{\theta }}_e\right)+\beta W_f{\tilde{\theta }}_e$$

(31)

where ${\tilde{\theta }}_e\in {ℝ}^3$ is the estimation error signal, defined as:

$${\tilde{\theta }}_e\triangleq {\theta }_e-{\widehat{\theta }}_e.$$

(32)

From (31), it can be shown that a mathematically useful, but unrealizable form of the prediction error $\epsilon (t)$ given in (24) can be written as [35]:

$$\epsilon =W_f{\tilde{\theta }}_e.$$

(33)

Based on the subsequent analysis, the following continuous least-squares update law ${\dot{\widehat{\theta }}}_e\in {ℝ}^3$ is employed for estimating the unknown parameters:

$${\dot{\widehat{\theta }}}_e\triangleq \mathrm{\Gamma }{W}_f^T\epsilon$$

(34)

where $\Gamma (t)\in {ℝ}^{3\times 3}$ is the least-squares estimation gain matrix which is designed as:

$${\dot{\mathrm{\Gamma }}}^{-1}\triangleq W_f^T{W}_f.$$

(35)

Remark 2.

If Q(t) is bounded, from (25), it can be seen that$Q_f(t),{\dot{Q}}_f(t)\in {\mathcal{L}}_{\infty }.$Similarly, if$W_e(t)\in {\mathcal{L}}_{\infty },$is bounded, from (27), it can be seen that$W_f(t),{\dot{W}}_f(t)\in {\mathcal{L}}_{\infty }$ [35].

Remark 3.

It should be noted that if${\mathrm{\Gamma }}^{-1}\left(t_0\right)$is selected to be positive definite and symmetric, then$\mathrm{\Gamma }\left(t_0\right)$is also positive definite and symmetric. Therefore, it follows that both${\mathrm{\Gamma }}^{-1}\left(t\right)$and$\mathrm{\Gamma }\left(t\right)$are positive definite and symmetric. The following expression can be obtained from (35):

$$\dot{\mathrm{\Gamma }}=-\mathrm{\Gamma }{W}_f^T{W}_f\mathrm{\Gamma }.$$

(36)

It can be easily seen from (36) that$\dot{\mathrm{\Gamma }}\left(t\right)$is negative semidefinite; therefore,$\mathrm{\Gamma }\left(t\right)$is always constant or decreasing. Hence, it follows that$\mathrm{\Gamma }\left(t\right)$is bounded [34,36].

Theorem 2.

The update law defined in (34) ensures that$‖{\tilde{\theta }}_e‖\to 0$as$t\to \infty$, provided that the following Persistency of Excitation condition [36] holds:

$${\kappa }_1{I}_3\le {\displaystyle \underset{t_0}{\overset{t_o+\delta }{\int }}{W}_f^T\left(\tau \right)W_f\left(\tau \right)d\tau \le {\kappa }_2{I}_3}$$

(37)

where ${\kappa }_1,{\kappa }_2,\delta \in ℝ$ are positive constants and I₃ is a standard 3-by-3 identity matrix.

Proof 2.

See Appendix B. □

3.2.2. Adaptive Controller Design

The parameter estimates, ${\widehat{\theta }}_e(t)$, are ultimately used to generate the desired trajectories; that is, the tracking objective for the control system. The current challenge is to now create a tracking controller for a system with uncertain parameters—this suggests an adaptive controller solution. To proceed with the development of an adaptive control law to achieve the control objective, we again consider c(t) as the virtual control input and define the error variables to rewrite the error dynamics as:

$$e_{a1}=q-q_d$$

(38)

$$e_{a2}=c-c_d.$$

(39)

As stated before, c_d(t) is the desired optimal trajectory for the concentration of the inhibitors. Note that q-dynamics can be rewritten as:

$$\begin{array}{cc}\hfill A_1\dot{q}& =A_{1}bq-A_{1}d{p}^{2/3}q-cq\hfill \\ & =W_1{\theta }_1-cq\hfill \end{array}$$

(40)

where $A_1\triangleq 1/S,$ $W_1=\left[\begin{array}{cc}q& -p^{2/3}q\end{array}\right]\in {ℝ}^{1\times 2}$ is a 1-by-2 regression vector containing measurable variables and ${\theta }_1={\left[\begin{array}{cc}{A}_{1}b& A_{1}d\end{array}\right]}^T\in {ℝ}^2$ is a 2-by-1 vector containing the uncertain parameters. After substituting (39) into (40), the following expression is obtained:

$$A_1\dot{q}=W_1{\theta }_1-e_{a2}q-c_{d}q.$$

(41)

After taking the time-derivative of (38) and multiplying both sides of the resulting expression by A₁, the following expression will be obtained:

$$A_1{\dot{e}}_{a1}=W_2{\theta }_2-e_{a2}q-c_{d}q$$

(42)

where (41) was utilized. In (42), $W_2=\left[\begin{array}{cc}{W}_1& -{\dot{q}}_d\end{array}\right]\in {ℝ}^{1\times 3}$ is a 1-by-3 regression vector containing measurable variables and ${\theta }_2={\left[\begin{array}{cc}{\theta }_1^T& A_1\end{array}\right]}^T\in {ℝ}^3$ is a 3-by-1 vector containing the uncertain parameters. Then, the desired trajectory for the concentration of the inhibitors can be designed as:

$$c_d=\frac{1}{q}\left(W_2{\widehat{\theta }}_2+k_{a1}{e}_{a1}\right)$$

(43)

in which k_a1 is a positive control gain. In (43), ${\widehat{\theta }}_2(t)\in {ℝ}^3$ is an estimate of the constant vector, ${\theta }_2.$ A parameter estimation error vector can be defined to quantify estimation performance as:

$${\tilde{\theta }}_2={\theta }_2-{\widehat{\theta }}_2.$$

(44)

Then the final dynamics for e_a1 will be:

$$A_1{\dot{e}}_{a1}=-k_{a1}{e}_{a1}-e_{a2}q+W_2{\tilde{\theta }}_2.$$

(45)

Following the same procedure, i.e., dividing both sides of the c-dynamics given in (1) by h and then substituting the results into the time derivative of (39), one can write:

$$A_2{\dot{e}}_{a2}=W_3{\theta }_3+u$$

(46)

where $A_2\triangleq h^{-1}\in ℝ,$ $W_3=\left[\begin{array}{cc}-c& -{\dot{c}}_d\end{array}\right]\in {ℝ}^{1\times 2}$ is a 1-by-2 regression vector containing measurable variables and ${\theta }_3={\left[\begin{array}{cc}m{A}_2& A_2\end{array}\right]}^T\in {ℝ}^2$ is a 2-by-1 vector containing the uncertain parameters. Based on the subsequent stability analysis and the control input, the concentration of angiogenic inhibitors, can be designed as:

$$u=-W_3{\widehat{\theta }}_3-k_{a2}{e}_{a2}+e_{a1}q$$

(47)

where k_a2 is a positive control gain and ${\widehat{\theta }}_3(t)\in {ℝ}^2$ is an estimate of the constant vector ${\theta }_3,$ which renders:

$$A_2{\dot{e}}_{a2}=-k_{a2}{e}_{a2}+q{e}_{a1}+W_3{\tilde{\theta }}_3$$

(48)

where

$${\tilde{\theta }}_3={\theta }_3-{\widehat{\theta }}_3.$$

(49)

Theorem 3.

For the tumor growth model given in (1), the concentration of the angiogenic inhibitors proposed in (47) ensures that the carrying capacity of the vascular network follows a desired trajectory, q_d(t), despite uncertainty about the values of parameters in the tumor dynamics with update laws:

$$\begin{array}{l}{\dot{\widehat{\theta }}}_2\triangleq {\gamma }_1{W}_2^T{e}_{a1}\\ {\dot{\widehat{\theta }}}_3\triangleq {\gamma }_2{W}_3^T{e}_{a2}\end{array}$$

(50)

where${\gamma }_1,{\gamma }_2\in ℝ$are positive constants. In other words, e_a1 and e_a2 go to zero. It also guarantees that all signals in the closed-loop system are bounded.

Proof 3.

See Appendix C. □

3.2.3. Optimum Trajectory Generation

The desired carrying capacity input to the adaptive controller will be a time-varying trajectory due to lack of knowledge about the model parameters, i.e., the cost function is changed based on parameter updates from the least-squares estimator. The gradient descent algorithm, also known as steepest descent, is a first-order optimization algorithm to minimize a function. For a function defined by a set of parameters, the gradient descent algorithm starts with an initial set of parameter values and then, by taking steps proportional to the negative of the gradient of the function, reaches the set of parameter values that minimizes the function. In general, it is considered as one of the simplest optimization algorithms. It is used in this overall approach to determine an optimum trajectory for the carrying capacity.

The desired trajectory, q_d(t), is needed as the input to the adaptive controller. This continuous optimum trajectory is designed to minimize the performance index given in (20), where the estimates of parameters S and d, obtained as described in Section 3.2.1, are utilized. For the minimization of the estimate of the performance index, a gradient descent algorithm is employed, which guesses the optimum value ${\overline{q}}_d[n]$ at each time step of the algorithm [10]. The output of the algorithm, ${\overline{q}}_d[n]$, is passed through a second-order, stable, and proper low-pass filter to generate the continuous and bounded signals $q_d(t)$ and ${\dot{q}}_d(t).$ The following filters are utilized:

$$q_d(t)=\frac{{\varsigma }_1}{{\varsigma }_2{s}^2+{\varsigma }_{3}s+{\varsigma }_4}{\overline{q}}_d[n]$$

(51)

$${\dot{q}}_d(t)=\frac{s{\varsigma }_1}{{\varsigma }_2{s}^2+{\varsigma }_{3}s+{\varsigma }_4}{\overline{q}}_d[n]$$

(52)

where ${\varsigma }_1, {\varsigma }_2, {\varsigma }_3, {\varsigma }_4\in ℝ$ are positive filter constants and $n\in \mathbb{Z}$ is a positive integer that represents the iteration step of the algorithm. At step n, the optimum trajectory holds the output, ${\overline{q}}_d[n]$, constant until the response of the closed-loop system, $q(t),$ has reached a steady-state near $q_d(t).$ A new target optimum, ${\overline{q}}_d[n+1]$, is then issued. In other words, the algorithm waits for certain thresholds to be satisfied before it proceeds to the next iteration. For instance, if $\left|q(t)-q_d(t)\right|\le {\overline{e}}_1,$ $\left|{\overline{q}}_d[n]-q_d(t)\right|\le {\overline{e}}_2,$ and $\left|q(t)-p(t)\right|\le {\overline{e}}_3,$ then n = n + 1, where ${\overline{e}}_1, {\overline{e}}_2, {\overline{e}}_3\in ℝ$ are threshold constants. Furthermore, the desired trajectory can be concluded to have converged when the gradient of $\widehat{J}(t)$, with respect to $q(t)$, is within a certain threshold. Once $q_d(t)$ has reached the termination threshold, the optimization algorithm stops updating ${\overline{q}}_d[n].$ As the performance index approaches its minimum value, the desired trajectory $q_d(t)$ and $u(t)$ approach $\widehat{q}*$ and $u*,$ respectively. $\widehat{q}*, u*\in ℝ$ are the estimates of the optimum values of $q(t)$ and $u(t),$ respectively, that result from the optimum seeking algorithm. As mentioned earlier, at steady-state $p(t)=q(t);$ therefore, $p(t)\to p_d(t)\to p*$ as $q(t)\to q_d(t)\to q*$, where $p_d(t)=q_d(t)$ and $\widehat{p}*=\widehat{q}*.$ Then, $\widehat{p}*$ is the estimated optimum tumor volume that can be realized by applying the estimated optimum drug dose, $\widehat{u}*.$

4. Simulation Results

A numerical simulation study was conducted to evaluate the proposed tumor treatment optimization technique using the MATLAB/Simulink environment. It should be noted that the least-squares estimator developed in Section 3.2.1, adaptive controller designed in Section 3.2.2, and the optimization algorithm described in Section 3.2.3 were run simultaneously, as shown in Figure 8. The parameter values were taken from [26] and are summarized in Table 1. Initial values were selected as p(t₀) = 10,000 mm³ and q(t₀) = 12,000 mm³ so that $q\left(t_0\right)>p\left(t_0\right),$ i.e., the carrying capacity of the endothelial cells is greater than the tumor volume; thus, the tumor volume is susceptible to an increase in the absence of a proper therapy.

The least-square estimator given in (34) was initialized as ${\widehat{\theta }}_e\left(t_0\right)=0.1{\theta }_e$ and the estimation gain matrix was initialized as ${\Gamma }^{-1}\left(t_0\right)=10I_3.$ The positive constant $\beta$ introduced in (25) was set to $\beta =0.05$. The adaptive update law was initialized as ${\widehat{\theta }}_a\left(t_0\right)=0_3$, where 0₃ denotes a 3-by-1 vector of zeros. The positive control gains were set to $k_{a1}=10, k_{a2}=5,$ and ${\gamma }_1={\gamma }_2=1\times {10}^{-12}.$ The initial guess, ${\overline{q}}_d[n],$ was selected to be 7500 mm³. Figure 9 shows the least-squares estimate of the unknown vector, ${\theta }_e$; it can be seen from this figure that all the parameters are accurately identified.

Figure 10 shows the time evolution of the drug dose, $u(t)$. It can be seen from Figure 10 that $u(t)$ converges to its optimum value $u*=0.1169 [\mathrm{conc}.]$.

The time evolution of the tumor volume, $p(t)$, and the carrying capacity of the endothelial cell volume, $q(t)$, are shown in Figure 11 and Figure 12, respectively. It can be seen from these figures that $p(t)$ and $q(t)$ converge to their respective optimum values $p*=6118{ \left[\mathrm{mm}\right]}^3$ and $q*=6118{ \left[\mathrm{mm}\right]}^3$ given in (5), respectively. It should be noted that if the initial conditions, $p(t_0)$ and $q(t_0)$, were chosen to be smaller values, then the convergence would have been faster.

The time evolution of the estimate of the performance index, $\widehat{J}(t)$, is shown in Figure 13. It can be seen from Figure 13 that $\widehat{J}(t)$ converges to $J*=12,247{ \left[\mathrm{mm}\right]}^3$.

Figure 14 shows the tracking error, $e_{a1}(t)$. It can be seen from the figure that the tracking error, $e_{a1}(t)$, is driven to zero.

Recently, an optimization method has also been proposed to analyze the mathematical model of cancer chemotherapy, which includes anti-angiogenic effects of cytotoxic agents in [37]. Similar to the study in [38], it provides a very useful discussion on modeling and optimization aspects of the problem, but they are open-loop controllers, i.e., both studies try to minimize to tumor volume by using the anti-angiogenic agents without feedback. Some useful methods with state feedback are given in [39], but they lack any optimization algorithm; on the other hand, lack of optimization in [39] leads to a very important advantage of continuous protocols. A very useful controller is given in [40], but it assumes all the system parameters are known. The main benefits of the controller given in existing papers are that it offers a control of the tumor volume and optimization algorithm to the drug dose. However, the parameter estimation scheme creates a disadvantage for the implementation of the control input signal, which is the concentration of the antiangiogenic agents. As seen from Figure 10, time variation of the drug dose has very high chattering. Even if this variation is continuous, a real-world clinical application will be discrete in two ways: drugs are not administrated continuously and the tumor volume is not measured continuously. Therefore, the chattering may affect the performance of the proposed tumor shrinking procedure, according to the selected drug-delivery frequency, especially in the first 6 months of treatment. Effects of continuous and discrete applications of the drug dosing controller has also been discussed in detail in [11].

5. Conclusions

A novel approach for the administration of anti-angiogenic tumor treatment by considering the nonlinear tumor dynamics and drug delivery pharmacokinetics and pharmacodynamics was presented. A performance index was formulated and then minimized to obtain an optimum tumor volume that would be realized; using an optimum drug dose and a nonlinear tracking controller was proposed to drive the system states to this optimum point. It was shown that if all the system parameters are known, the tumor volume can be driven to its optimum value exponentially fast.

The problem is more complicated if the tumor model parameters are unknown; it was shown that the performance index can be estimated, using a least-squares strategy, and the optimum tumor volume trajectory can be generated numerically. The proposed least-squares estimation scheme proves that the estimation errors can be driven to zero, provided that the persistency of excitation condition is satisfied. An adaptive tracking controller then drives the system to this optimum. The proposed adaptive scheme ensures that the tumor volume and the carrying capacity of the vascular network go to their desired values globally and asymptotically. Numerical simulation results were presented to show the performance and feasibility of the proposed estimation and control approaches. The simulations demonstrate that the developed estimation and control strategies successfully minimize the tumor volume, along with the carrying capacity of endothelial cells, with an optimum drug dose, despite of the lack of knowledge about system parameters. It can be concluded from the simulation study that the proposed tumor minimization technique can efficiently reduce the tumor volume with an optimum drug dose.