Numerical solution for a problem arising in angiogenic signalling

1.   Introduction

A growing tumor needs a steady supply of oxygen and nutrients for cell duplication, thus the growth of new vessels of a tumour are understood to be stimulated by mainly the principal stimulus, known as the angiogenic switch. In most cases it appears to be because of oxygen deprivation, although other stimuli such as inflammation, oncogenic mutations and mechanical stress may also play a role. Thus, such angiogenic switch leads to tumor expression of pro-angiogenic factors and increased tumor vascularization ^[12].

Initially, during avascular growth, it is provided through the surrounding environment. As the tumor becomes larger these mechanisms become inadequate and tumor cells enter the dormant stage of the cell cycle. As a consequence, vascular endothelial growth factors (VEGF) are released that stimulate the formulation of new blood vessels and capillaries in order to supply the tumor with needed nutrients. This process is called tumor angiogenesis. Hence tumor anti-angiogenesis is a treatment approach for cancer that aims at depriving the tumor of this vasculature.

Ideally, without an adequate support network, the tumor shrinks. Anti-angiogenic treatment was already proposed in the early seventies by J. Folkman ^[13], but became medically possible only with the discovery of the inhibitory mechanisms of the tumor in the nineties ^[9,20]. It brings in external anti-angiogenic agents to disrupt the growth of endothelial cells which form the lining of newly developing blood vessels and capillaries. The intent to directly kill tumour cells or prevent their proliferation has in many cases proved futile as the kinetic understanding of tumour control and sensitivity characteristics reveal that tumour population is far from stable. Therefore, since the tumour vasculature does not exploit tumour cell sensitivities, Hahnfeldt et al. ^[14] realized that it relies on tumour suppression consequent to inhibition of associated vasculature. This has paved the way for antiangiogenic therapy to control an exceptionally heterogeneous, unconstrained tumour population via a relatively homogeneous and constrained endothelial population as it allows one to disregard a vast array of spatial and temporal details of tumour cell expression. As a consequence, no clonal resistance to angiogenic inhibitors has been observed in experimental cancer ^[2].

Since developing drug resistance all too often is the limiting factor in conventional chemotherapy treatments as cancers have a formidable capacity to develop resistance to a large and diverse array of chemical, biologic, and physical anti-neoplastic agents, Kerbel ^[18], claimed that it can be largely traced to the instability of the tumour cell genome, and the resultant ability of tumour cell populations to generate phenotypic variants rapidly. Therefore, anti-cancer strategies should be directed at eliminating those genetically stable normal diploid cells that are required for the progressive growth of tumours. Hence tumor anti-angiogenesis has been called a new hope for the treatment of tumors ^[21]. Although these high hopes have not been realized in practice, there still strong interest and active research on tumor anti-angiogenesis as a method that normalizes the vasculature ^[15,16] and thus, when combined with traditional treatments like chemotherapy or radiotherapy, enhances the efficiency of these procedures.

Apart from formulating a class of mathematical models for tumor anti-angiogenesis as optimal control problems, Ledzewicz and Cardwell ^[30] considered the fact on how to schedule an a priori given amount of anti-angiogenic (e.g., vessel disruptive) agents in order to minimize the tumor volume ^[37,38], they also analyzed these models for a class of mathematical models that include, based on a model that was developed and biologically validated by Hahnfeldt, Panigrahy, Folkman and Hlatky ^[14]. The principal state variables are the primary tumor volume, $p$, and the carrying capacity of the vasculature, $q$, where the latter is a measure for the tumor volume sustainable by the vascular network. The dynamics describes the interactions between these variables and the tumor volume $p$ changes according to some growth function dependent on the variable carrying capacity $q$, where the $q$-dynamics consists of a balance of stimulatory and inhibitory effects. While significant modeling changes are made in the dynamics for the vascular support in this model, the solutions to the optimal control problem are in fact qualitatively identical.

Ledzewicz et al, ^[29] considered two mathematical models for tumour anti-angiogenesis in which one model was originally formulated in ^[14] whereas, the other model is a modification of the model by ^[11] considered as optimal control problem with the aim of maximizing the tumour reduction achievable with an a priori given amount of angiogenic agents. They argued that depending on the initial conditions, the optimal controls may contain a segment along which the dosage follows a so-called singular control, a time-varying feedback control. Thus, the efficiency of piecewise constant protocols with a small number of switchings is investigated through comparison with the theoretically optimal solutions. It is shown that these protocols provide generally excellent suboptimal strategies that for many initial conditions come within a fraction of $1\%$ of the theoretically optimal values. When the duration of the dosages are a priori restricted to a daily or semi-daily regimen, still very good approximations of the theoretically optimal solution can be achieved.

Hahnfeldt et al. ^[14] described the growth of a tumour assuming that tumour growth is strictly controlled by the evolution of the vascular network that supplies oxygen and nutrients to tumour cells and noticed that it provides a framework to represent the effects of antiangiogenic therapies. In their paper, some possible modifications of their model are proposed, and conditions that guarantee the eradication of the tumour under a regimen of periodic antiangiogenic therapy are derived. The model variants considered assume the potential doubling time of the vasculature to be constant, and subdivide the endothelial cell pool, which is involved in angiogenesis, in resting and proliferating cells allowing for a more detailed description of drug effects.

In ^[31] considered the problem of minimizing the tumor volume with a priori given amounts of anti-angiogenic and cytotoxic agents. For one underlying mathematical model, optimal and suboptimal solutions are given for four versions of this problem: the case when only anti-angiogenic agents are administered, combination treatment with a cytotoxic agent, and when a standard linear pharmacokinetic equation for the anti-angiogenic agent is added to each of these models. It is shown that the solutions of the more complex models naturally can be built on the simplified versions. This gives credence to a modeling approach that starts with the analysis of simplified models and then adds increasingly more complex and medically relevant features. Furthermore, for each of the problem formulations considered here, there exist excellent simple piecewise constant controls with a small number of switchings that virtually replicate the optimal values for the objective.

Ledzewicz et al. ^[27] analyzed the scheduling of angiogenic inhibitors as an optimal control problem for a mathematical model for tumor anti-angiogenesis proposed by Ergun et al. ^[11] with a logistic growth function modeling tumor growth. It is shown that optimal controls are bang-bang with at most two switchings.

Sebastien ^[36] introduced a phenomenological model for anti-angiogenic therapy in the treatment of metastatic cancers, which is a structured transport equation with a nonlocal boundary condition describing the evolution of the density of metastases, that at first were analyzed at the continuous level. He presented the numerical analysis of a Lagrangian scheme based on the characteristics whose convergence establishes existence of solutions and proved an error estimate that used the model to perform interesting simulations in view of clinical applications.

In ^[7] anti-angiogenic therapy is considered to make a notable difference in every day cancer treatment. While the technique has many advantages the cost of treatments are often expensive due to the non-personalized administration medical protocols. Thus, in their paper, Czako et al. ^[7] considered a model based solution which aims to lower the medical expenses during the treatment by creating personalized administration plans with the help of control engineering.

A contribution to the theory of optimal control can be traced in ^[17], introduction to nonlinear programming, where the numerical methods for optimal control problem are considered in ^[1], whereas, in ^[33] explains how optimal-control problems can be solved with a common spreadsheet such as Microsoft Excel.

Dontchev and Hager ^[10] analyzed the Euler approximation to a state constrained control problem and showed that if the active constraints satisfy an independence condition and the Lagrangian satisfies a coercivity condition, then locally there exists a solution to the Euler discretization. Their error is bounded by a constant times the mesh size. Their analysis utilizes mappings of the discrete variables into continuous spaces where classical finite element estimates can be invoked.

In ^[19] considered the reduction of the effects of modeling impr$\acute{e}$cisions, that is, the actually measured state variable is used as the starting point in the next cycle within a horizon-length cycle, where a cost function is minimized under a constraint that mathematically represents the dynamic properties of the system under control. Thus, the nonlinear programming approach, the state variables as well as the control signals are considered over a discrete time-resolution grid, and the solution is computed by the use of Lagrange's reduced gradient method. They have suggested that instead of exerting the estimated control signals, the estimated optimized trajectory is adaptively tracked within the given horizon and they found out that the transients of the adaptive controller that appear at the boundaries of the finite-length horizons reduce the available improvement in the tracking precision. In contrast to the traditional Receding Horizon Control, in which decreasing horizon length improves the tracking precision.

The shortage of limited resources in every undertaking is a very serious concern to the survival of human kind. Thus, in this paper, we would like to provide an adequate analysis of the optimal problems which arise as a result of agiogenic signalling of tumor cells. To this end, it is evident that in the literature more work required to be done as far as qualitative and quantitative features of these type of problem are concerned. In turn, this can ensure that the implementation of such models in real life are indeed cost effective across all stake holders. Therefore, instead of defining admissible singular arcs as in ^[14] without presenting models' solutions, thus, our focus in this paper is to analyse the equilibrium state of the models, use their derived singular arcs to implement a robust numerical method based on the qualitative behaviors of the the models.

The rest of the paper is arranged as follow, Section 2 states the problem description, whereas Section 3 highlights the Hamiltonian and Lagrange multipliers. We analyse the equilibrium state of the models in Section 4 and state the singular controls for the models in Section 5. Numerical method and the stability of the method are presented in Section 6 and 7, respectively. We discuss our numerical results in Section 8 and conclude the paper with Section 9.

2.   Problem description

Let $p, \xi, q, \gamma, u, a, A, b, d, \mu q$ denote the primary tumor volume, tumor growth parameter, endothelial support, anti-angiogenic killing parameter, treatment with an anti-angiogenic agent, a priori set maximum dosage, positive constant, birth rate, death rate, net balance between endothelial cell proliferation and loss to the endothelial cells through natural causes such as death and the parameter $\theta\in[0, 1]$. Then, to reduce the volume ($p$) of a tumour efficiently results into the maximization of the tumour volume reduction achievable with an apriori amount of angiogenic inhibitors ^{[11,22,24,26,28]}

for a free terminal $T$, minimizes the value $p(T)$ subjects to the dynamics

with initial conditions $p(0) = p_0>0, q(0) = q_0>0, y(0) = 0$ ^[14]. Equation (2.1) together with (2.2) is an optimal control problem. Therefore, in the next section we determine the Hamiltonian and Lagrange multipliers of the optimal control problem.

3.   Hamiltonian and Lagrange multipliers

The Pontryagin maximum principle ^[3,4,34] enables us to determine the necessary conditions for optimality of a control $u$. Thus, for a row-vector ${\bf \lambda} = (\lambda_1, \lambda_2, \lambda_3)^t\in\mathbb{R}^3$ the Hamiltonian $H: = H({\bf \lambda}, p, q, u)$ is

where, $S$ and $I$ denote endogenous inhibition, stimulation terms. Therefore, the individual Hamiltonians ^[14] corresponding to equation (2.2) are

over all Lebesgue measurable functions $u:[0, T]\rightarrow[0, a]$, for which the corresponding trajectory satisfies $y(T)\leq A$ and the transversality conditions are

Let $\bar{x}: = (p, q, y)$, by Samaee et al. ^[35], we obtain

where,

obtained through equations in (2.2). Applying equation (3.3) to model $H_1$ we obtain

which implies that

Equation (3.7) is also obtained for other models. Solving equation (3.7), we obtain

where $C$, is a constant of integration. Using the transversality conditions we obtain

4.   Equilibrium state

In order to develop the robust numerical methods it is necessary to analyse the steady state behaviour of these models. Therefore, in the next subsections we deduce the stability conditions of the models.

4.1.   Model $I_{\theta}$

For this model, we let

then

We see that $H_p = H_q = H_u = 0$, where the subscripts imply the partial derivatives with respect to a subscript $p$, $q$ and $u$ in that order.

Solving for critical point $q^*$ in equation (4.1), we find that

But we know that $q^*\geq p^*$, then this enables us to write

then, $p^*>0$ as $p^* = 0$ is not admissible. Therefore,

which we solve and obtain

From equation (4.2), we obtain the non-zero entries of the Jacobian matrix $J_{I_{\theta}}: = J_{ij}$ for $i = j = 1:3$ as

Using the concept of numeric-analytic dissipativity condition^[6], we obtain the characteristic equation

from $\frac{1}{2}(J+J^t)$. This implies that the model is stable if

which implies that

4.2.   Adjoint for model $I_{\theta}$

Since the Hamiltonian $H = H(\lambda, p, q, u)$, where the $\lambda$'s are constants multipliers then

because $d\lambda/dt = 0$ and by the stationary condition we have $\partial H_{I_{\theta}}/\partial u = 0$. Therefore, for the steady state equation in (4.10) becomes

Using equation (4.11) we find the corresponding critical points by linearizing the Jacobian matrices as follow

and

Solving for the critical point $q^*$ in (4.13) we find

The Jacobian matrix is

Therefore, the adjoint is stable if and only if and the eigenvalues are

which implies that

4.3.   Model $H_e$

Applying the same procedures as in the above section we have,

then

and we also see that H_p=H_q=H_u=0, where the subscripts denote the partial derivatives with respect to p, q and u, respectively.

Then from the second equation in (4.15) we see that

which implies that $bq^{-1/3}-dq^{1/3}-\mu = 0$, as $q^*\neq 0$. This implies that

because $u^*=0$. Take $q^*\geq p^*$ and the non-zero entries of the Jacobian matrix $J_{E}:=J_{ij}$ where $i=1:3$ and $j=1:3$ are

which implies that the model is stable if and only if

4.4.   Adjoint for model $H_e$

Let $H=H(\lambda,p,q,u)$, then

as $d\lambda/dt = 0$ and by the stationary condition we see that $\partial H/\partial u = 0$. Thus, for the steady state we have

In view of equation (4.19) we have

and from

which implies that

The corresponding Jacobian matrix $J_{H_{E}}:=J_{ij}$ for $i=j=1:3$ is

where, the non-zero entries are

Therefore, the adjoint of this model is stable if

4.5.   Model $H_1$

We let

so that

where, we see that $H_p=H_q=H_u=0$. The subscripts imply the partial derivatives with respect to $p$, $q$ and $u$, respectively. Therefore,

and the Jacobian matrix $J_{H_1}:=J_{ij}$ for $i=j=1:3$ is

where, the non-zero entries are

In view of equation (4.22), we see that,

The first equation in (4.26) requires that $-\xi p=0$ or $\ln\left(p/q\right)=0$. However, based on the construction of this model, neither $\xi\neq 0$ nor $p\neq 0$, then the only choice is

However, further basic requirement on this model is such that $\ln\left(p/q\right)$ should be a decreasing function and this is only possible if $q^*\geq p^*$. Solving for $q$ in the second equation in (4.26) we obtain

But $q^*\geq p^*$, then in view of equation (4.28), we see that

which enables us to rewrite equation (4.28) as

Hence, the model is stable if and only if

4.6.   Adjoint for model $H_1$

For this model we have

as $d\lambda/dt=0$ and by the stationary condition we see that $(\partial H/\partial u)=0$. Thus, for the steady state equation (4.32) becomes

Using equation (4.32) we have,

In view of equation (4.32) we see that

which implies that

whereas

which implies that $p^*/q^* = 0$, which is possible if $p^* = 0$ and $q^*\neq 0$. The Jacobian matrix is

in which we see that

and

Thus, the adjoint of this model is unconditionally stable.

5.   Singular controls for the models

Since the Hamiltonian ($H$) is linear in $u$, then minimizing the control requires that $u = 0$ or $u = a$ ^[3]. This is known as the bang controls. In view of equations in (3.1), we obtain the switching function ($\Phi$) as

such that the singular control is ^[3]

where, for the three models we have the optimal singular arcs

6.   Numerical method

Not withstanding the associated optimal synthesis of the models considered in this paper, but it is evident from the stability structures of the continuous models that reliable numerical method should be developed. Thus, in order to accomplish the development of a robust numerical method for optimal problems arising as a result of angiogenic signalling, we believe we first have to consider the existing numerical methods for these types of models. However, in this paper, we consider only one type of the numerical method for the models. Thus, we sub-divide the interval [0, T] into equal pieces with specific points of interest

where $N$ is a positive integer denoting the number of sub-intervals. Since the total-enumeration methods or linear programming techniques can be used to solve optimal control problems such the one in ^[1], because such methods fail to capture the associated optimality, adjoint equation and the transversality condition. Therefore the only applicable methods are Runge-Kutta or adaptive schemes and the boundary value problems such as shooting method ^[5,8]. Hence, following the Forward-backward sweep method ^[32] then the optimal control problem is implemented as we have shown here below.

6.1.   Forward-backward sweep method (FBSM) for $I_{\theta}$

Set the $flag = -1$, then define the step-size $h = 1/N$, and initialize the controls, states and the adjoints with their initial conditions, we have

and Step:. WHILE ($flag<0$) do the following steps.

Step 1a. oldu$= u$; oldp$= p$; oldq$= q$; oldy$= y$; oldlambda1 = $\lambda_1$;

oldlambda2$= \lambda_2$; oldlambda3$= \lambda_3$;

Step 2a

FOR $i = 1, 2, \cdots, N$ set

STOP

Step 3a

FOR $i = 1, 2, \cdots, N$ and $j = N+1-i$ set

7.   Stability analysis of FBSM

Basically the FBSM first solves the state equation with a forward in time Runge-Kutta method, then solves the costate equation backwards in time with the Runge-Kutta method and then updates the control. Then, stability analysis should follow the procedures carried out when one determine the condition of the Runge-Kutta method. Since we have impose the numeric-analytic dissipativity condition ^[6] to the models eigenvalues, then FBSM is A-stable.

8.   Numerical result

Based on the initial conditions $p_0 = 8.00, q_0 = 4.00, u_0 = 0.1$, parameter values $\xi = 0.084, b = 5.85, d = 0.00873, \mu = 0.02; \gamma = 0.01, \theta = 0.1, \delta = 0.1$ (^[26]), we implemented the Forward-backward sweep method (FBSM) for the systems in (2.2) and (3.1) as shown up for the case of $I_{\theta}$, where the numerical approximations are presented in Figure 1 and for the remaining two models the results are presented in Figure 2 and Figure 3. Our aim in this paper is to present the numerical solutions of the three models, we have considered. Thus, we see that the control ($u$) and angiogenesis ($q$) increases monotonically but remain bounded, except for the $H_1$ model. We also see that the tumour volume ($p$) decreases and increases eventually. This is due to the ever growing agiogenesis system of the tumor. Such phenomena is also evident for model $H_1$. The above-mentioned behaviours remain the same, when we perturb the initial values and for an increased values of $T$.

9.   Conlusion

In view of the problem description, Hamiltonian and Lagrange multipliers, we were able to deduce the multipliers for these models. We have also established the stability conditions for each model which in turn guaranteed the stability of the Forward-backward sweep method. In doing so, we believe that this can enable us to attain most features of each model which can give deeper insight of the properties of the models. Since the authours in ^[23,27,28] were mainly interested in attaining the singular arc of the models, it is important to combine the defining element and all the syntheses of optimally controlled trajectories qualitatively and quantitative with the associated solution to a problem. Therefore, this paper should be viewed as a first attempt to combine singular arc with their associated solutions of the optimal problems. Hence, our future research direction is to extend the paper to higher dimensional space, with the inclusion of the spatial effects.

Acknowledgments

We would like to thank the University of the Western Cape for the NRF support.

Conflict of interest

All authors declare no conflict of interest.