CONTROL PROBLEMS ARISING IN COMBINED ANTIANGIOGENIC THERAPY AND RADIOTHERAPY Andrzej Swierniak Department of Automatic Control Silesian University of Technology ul. Akademicka 16, 44-100 Gliwice Poland [email protected] 1) The discovery of basic fibroblast growth factor as the first pro-angiogenic molecule 2) The discovery of vascular endothelial growth factor and its receptor tyrosine kinases on activated endothelial cells 3) The discovery of angiopoietins and their tyrosine kinase receptors 4) The discovery of endogenous inhibitors of angiogenesis 5) The discovery of additional molecular markers in newly formed blood vessels 6) The development of quantitative assays for angiogenesis 7) Recognition of the prognostic significance of tumor angiogenesis 8) Lack of acquired resistance to direct acting antiangiogenic drugs 9) The discovery of the impact of angiogenesis on liquid hematologic malignancies 10) The discovery of the accidental antiangiogenic effects of various conventional or new anticancer drugs

ABSTRACT A model of combined antiangiogenic therapy and radiotherapy of tumors is proposed and analyzed. The model is a modified Hahnfeldt model (proposed by d’Onofrio and Gandolfi) and we discuss its properties including its asymptotic behaviour under constant and periodic therapy. Then we propose an optimization problem in the final horizon which leads to the optimal treatment protocols. We find necessary conditions of its solution and discuss their biological meaning. KEY WORDS Modelling, optimization, radiotherapy, tumor angiogenesis, nonlinear systems

1. Introduction During progression of tumor molecular factors called activators (stimulators) and inhibitors (blockers) of angiogenesis are released by tumor to develop its own vascular network which enables its growth and in the next stage determines possibility of cancer metastasis. Since this network is necessary for tumor development, in late sixties of the last century a new anticancer therapy was proposed target of which was not directly the cancer cells but the new born vasculature. This therapy is known as antiangiogenic therapy and the idea is to reduce the tumor volume reducing its vasculature. It has been first time hypothesized by Folkman [1] more than thirty years ago. The most important obstacle against successful chemotherapy is drug resistance acquired by cancer cells while the normal tissues retain sensitive to the drugs. This negative feature of chemotherapy may be used as an advantage in the antiangiogenic therapy which is directed towards special part of normal tissues and only indirectly destroys tumor cells and it is why it has been called by Kerbel [2] a therapy resistant to drug resistance. Tumor angiogenesis belongs to the most inspiring areas of cancer research in oncology. Kerbel presents 10 significant reasons for the explosive growth in tumor angiogenesis research and development of antiangiogenic drugs:

555-015

Nevertheless still the most important constraint in efficient antiangiogenic therapy is the accessibility of antiangiogenic agents. This is why the rational anticancer therapy should contain combination of antiangiogenic therapy with more standard modalities of anticancer treatment for example radiotherapy. We discuss a model of such combined therapy. The complexity of the process of vascularization results in the complicated models (see e.g. [3]) applicable for simulation of the process but completely not useful in synthesis or even analysis of therapy protocols. The exception is a class of models proposed by Hahnfeldt et al [4] who suggested that the tumor growth with incorporated vascularization mechanism can be described by Gompertz type or logistic type equation with variable carrying capacity which defines the dynamics of the vascular network. Roughly speaking the main idea of this class of models is to incorporate the spatial aspects of the diffusion of factors that stimulate and inhibit angiogenesis into a non-spatial two-compartmental model for cancer cells and vascular endothelial cells. This type of model or more precisely its

113

The growth is bounded by:

modification proposed in [5] will be used by us in our study. In [5] it has been proved that using sufficiently high doses of antiangiogenic drugs we are able to annihilate completely the vascular network of the tumor and indirectly eradicate the tumor itself. It can be reached not only using a constant dose of the drug but also by periodic therapy more reasonable from clinical point of view. Nevertheless since the results have an asymptotic character it means that the process of eradication is theoretically infinite and the same the patient once treated by the antiangiogenic therapy should remain under such control to the end of his life. To overcome this difficulty we have proposed [6] to optimize the therapy in finite horizon which is the subject of section 3. We formulated necessary conditions of optimality based on Pontryagin Maximum Principle and explained the results in terms of therapy protocols. Here we extend our results onto the case when two types of therapy are used. The optimization problem for yet another modification of Hahnfeldt model was solved in [7] by Ergun and coworkers. Their model has a drawback that there exists only one cross-coupling between the two compartments which results for example in trivial stability conditions. The rigorous treatment of this model only in the case of antiangiogenic therapy has been recently presented in [8].

N ∞ = N 0 eα / β

called in population dynamics the carrying capacity. The same solution is obtained when we use non-linear Gompertz equation in the form: N& = α N − β N ln N / N 0 , N (0) = N 0 , or N& = − β N ln N / N ∞

N& / N = − β ln N / N ∞ ≈ 1/ PDT

The simplest model of population kinetics for cancer tissues is given by Malthusian growth which assumes exponential relationship between a size of the population and time. The dynamics is described by the equation:

N = N 0 eat , a = ln 2 / PDT

(2)

N = N 0 eα / β (1−e

−βt

)

(8)

Although equation (8) looks similarly to the second of equations (6) but now the carrying capacity is not constant but varies with changes of the volume of the vessels. The dynamics of the growth of this volume represented by its PDT depends on the stimulators of angiogenesis (SF), inhibitory factors secreted by tumor cells (IF) and natural mortality of the endothelial cells (MF):

PDTk = f ( MF , SF , IF )

with N denoting the size of the population and a Malthusian parameter defined by the inverse of the potential doubling time (PDT). The unlimited growth in this model can be avoided if we introduce a varying coefficient a (t ) as it is in Gompertz model:

N& = a(t ) N , N (0) = N 0 , a& = − β a, a (0) = α ⇒

(7)

Hannfeldt [4] proposed to treat the carrying capacity which constrains the tumor growth as a varying tumor volume sustainable by the vessels and roughly proportional to the vessel volume: N ∞ = K , N& / N = − β ln N / K

(1)

(6)

It enables defining PDT by:

2. A Model of Cancer Growth Including Vascularization, Radiotherapy and Antiangiogenic Therapy

N& = aN , N (0) = N 0

(5)

(9)

In [5] it has been assumed that the inverse of PDT is the sum of these three factors i.e.

1 / PDTk = MF + SF + IF

(3)

(10)

The spontaneous loss of functional vasculature represented by MF (e.g. through natural mortality of the endothelial cells) is supposed to be negative constant, the stimulatory capacity of the tumor upon inducible vasculature represented by SF (e.g. through angiogenic factors like vascular endothelial factor) is found to grow at rate K b N c slower than the endogenous inhibition of

(4)

114

previously generated vasculature represented by IF (e.g. through endothelial cell death or disaggregation) where: b+c ~2/3

x = ln N / N * , y = ln K / K * , x* = y* = 0, τ = β t , ϑ = (γ − µ ) / β , x′ = dx / dτ , y ′ = dy / dτ ,

(11)

It could be transformed into the following quasi-linear system:

It results from the assertion that tumor driven inhibitors from all sites act more systematically whereas tumorderived stimulators act more locally to the individual secreting tumor site From the other hand analysing a diffusion-consumption equation for the concentration of stimulator or inhibitor inside and outside the tumor, Hahnfeld et al concluded that inhibitor will impact on target endothelial cells in the tumor in a way that grows ultimately as the area of the active surface between the tumor and the vascular network which in turn is proportional to the square of the tumor diameter. It leads to the conclusion that IF is proportional to the tumor volume in power 2/3 since volume is proportional to the cube of the diameter. The expression for K suggested in [4] has therefore the following form:

K& / K = γN / K − (λN 2 / 3 + µ )

x′ = y − x, y ′ = ϑ (1 − e2 / 3 x )

z = y − x, x ′ = z , z ′ = − z − ϑ (e 2 / 3 x − 1)

Application of antiangiogenic therapy can be incorporated to the model by a factor increasing multiplicatively the mortal loss rate of the vessels:

of stimulation, inhibition and natural mortality, respectively. The modification of this model proposed in [5] which also satisfies Hahnfeldt’s suggestions given by (11) assumes that the effect of SF and MF on the inverse of PDT K is constant while the IF is proportional to the active surface of the area of tumor being in contact with the vascular network and the same to the square of the tumor radius:

(13)

K& / K = γ − (λN 2 / 3 + µ )

(14)

(18)

Using the first method of Lyapunov it could be found that the equilibrium point is locally asymptotically stable. Moreover it is globally asymptotically stable that could be proved [5] using the Lurie type Liapunov function and the second method of Lyapunov [9].

γ , λ , µ being constant parameters representing the effect

N& / N = − β ln N / K

(17)

or:

(12)

K& / K = γ − (λ N 2 / 3 + µ + η u (t )),

(19)

where u (t ) denotes the dose of the agent scaled to its effect on vascular network and η is a constant parameter. For the constant dose U, the equilibrium points take the form:

N * = K * = ((γ − µ − ηU ) / λ )3 / 2

(20)

which according to the conditions of stability given in [5] leads to the conclusion that for:

ηU ≈ γ − µ ⇒ K * → 0

This system has a nontrivial equilibrium point (N*,K*):

N& / N = K& / K = 0 ⇒ N * = K * = ((γ − µ ) / λ )3 / 2

(16)

(15)

(21)

The form of condition (21) results from the suggestion that even if the dose is not exactly equal to the value found from the equilibrium condition the convergence to 0 takes place. In other words the vascular network and in turn the tumor can be eradicated. In [5] it has been proved that the same effect could be reached for periodic therapy with mean value satisfying condition (21) or greater.

The model is strongly nonlinear but by logarithmic change of variable and some scaling transformation:

115

The effect of radiotherapy should be included in both compartment because radiation destroys both cancer and normal tissues. The classical LQ model (e.g. [10]) assumes that the damage to DNA, which is the principal target for the radiation has two components: a linear one that is a consequence of a simultaneous break in both DNA strands caused by a single radiation particle, and a quadratic one that is the result of two separate but adjacent breaks in different strands caused by two different particles. In our model we omit this second term and introduce only linear (in dose) effect into both equations. It leads to the following model: N& / N = − β ln N / K

− ψv

with known parameters: U m ,Vm , Ξ, Φ.

constraining The integral

constant constraints

imposed on control variables although similar in form have different meaning. For radiotherapy it measures the feasible cumulated negative affect of the radiation while the one for antiangiogenic agent represents mostly the shortage in the availability of the agent and only in part the possible side effects of the drugs (not sufficiently recognized yet). Due to isoperimetric form of the problem it could be transformed into the problem with the integral part of the performance index instead of the integral constraint on the control:

(22)

Tk

Tk

0

0

J = N (Tk ) + r ∫ u (t )dt + s ∫ v(t )dt

(25)

0 ≤ u ≤ U m ,0 ≤ v ≤ V m K& / K = γ − (λN

2/3

+ µ ) − ηu − ξv

(23) This problem in turn can be approximated by the following quasi-linear problem in logarithmic variables:

• Where v(t) denotes the dose of the radiotherapy scaled to its effect on tumor and normal tissues and ξ and ψ are constant scaling parameters. Of course the additional radiotherapy supports the effect of antiangiogenic therapy. Moreover the effect of tumor eradication may be achieved easier and faster although still the theoretical results based on the theory of stability have asymptotic form.

Tf

I = gx(T f ) + hy (T f ) + r ∫ u(τ )dτ + 0

Tf

(26)

+ s ∫ v(τ )dτ ,0 ≤ u ≤ 1, T f = Tk β 0

3. Optimization of Therapy in Finite Horizon x ′ = y − x − εv, y ′ = ϑ (1 − e 2 / 3 x ) + νu + ςv,ν = = −η / β , ε = ψ / β , ς = −ξ / β

The results presented in the previous section are asymptotic which in turn leads to the conclusion that therapy should be permanent to be successful. One way to overcome this difficulty is to find the best (in some sense) therapy protocol which could be realized in finite time. The reasonable solution is to formulate optimal control problem for the system given by the proposed model and the control objective which adequately represents the primary goal of the therapy. In [7] and [8] the optimal control problem for the simplified model in which there is only one way cross-coupling between the two components of the model was presented for a free terminal time. We propose to optimize the protocol in the fixed finite time of therapy with the primary goal which is to find the control maximizing TCP (treatment cure probability) that leads to the following equivalent form of an optimal control problem: Tk

Tk

0

0

J = N (Tk ), ∫ u (t )dt ≤ Ξ, ∫ v(t )dt ≤ Φ

(27)

The weight coefficients h, g, r may change in broad ranges depending on the type of therapy used and the strength of the integral constraint. The additional term related to the volume of vascular network may be regarded as yet another constraint imposed on the possible dynamics of the system. On the other hand by the choice of the weighting coefficients we obtain a new possibility of analysis of the mutual dependence between the tumor growth and the volume of the vascular network. Thus it is reasonable to provide an extensive analysis of their effect on the solution of the optimal control problem. Necessary conditions of optimality can be found using Pontryagin maximum principle [11] for Hamiltonian and adjoin variables p, q defined as:

(24)

H = ru + sv + νqu + ζqv + p ( y − x − εv) + qϑ (1 − e 2 / 3 x )

0 ≤ u (t ) ≤ U m ,0 ≤ v(t ) ≤ Vm

p ′ = p + 2 / 3qϑ e2 / 3 x , p (T f ) = g , q ′ = − p, q(T f ) = h

116

(28)

It leads to the following switching functions and bangbang control law:

q = − r /ν , u = {

1 0

⇐ min H

Acknowledgment The author was partly supported by grants internal grant of Silesian University of Technology n.BK208/Rau1/2006.

(29)

References 1 pε − qζ = s, v = { ⇐ min H 0

[1] J. Folkman: Tumor angiogenesis: therapeutic implications, N. Engl. J. Med. 295 (1971), 1182-1186. (30)

[2] R.S. Kerbel: A cancer therapy resistant to resistance, Nature 390 (1997), 335-340. [3] A.R.A. Anderson and M.A.J. Chaplain: Continuous and discrete mathematical models of tumor induced angiogenesis, Bull. Math. Biol. 60 (1998), 857-864.

Rewriting the adjoint equation in the form of scalar second order ODE we have: q ′′ − q ′ + 2 / 3qϑ e2 / 3 x = 0, q (T f ) = h, q ′(T f ) = − g

[4] P. Hahnfeldt, D. Panigraphy, J. Folkman J., and L. Hlatky: Tumor development under angiogenic signaling: A dynamic theory of tumor growth, treatment response and postvascular dormacy, Cancer Res. 59 (1999), 47704778.

(31)

The important finding is that singular arcs are not feasible since there are no finite intervals of constant solutions to the adjoint equation. This leads to the conclusion that intermediate doses of the drug are not optimal and that the optimal protocol contains only switches between maximal dose and no drug intervals. It differs the obtained solution from those obtained in [7] and [8]. The form of conditions allows to find recurrently the solution of the TPBVP composed of the state and co-state equations with bangbang control found from the switching condition by using for example shooting algorithm. Depending on the values of system parameters there may be only one switch of the optimal control or the solution may be quasi periodic.

[5] A. d’Onofrio and A. Gandolfi: Tumour eradication by antiangiogenic therapy analysis and extensions of the model by Hahnfeldt et al (1999) Math. Biosci. 191 (2004), 159-184. [6] A. Swierniak, A. D’Onofrio, A. Gandolfi: Antiangiogenic therapy as a control problem, Proc. Int. IASTED Conference on Biomechanics, Bio Mech 2006, paper n. 542-801. [7] A. Ergun, K. Camphausen, and L.M. Wein: Optimal scheduling of radiotherapy and angiogenic inhibitors, Bull. Math. Biol. 65 (2003), 407-424.

4. Conclusion

[8] U. Ledzewicz and H. Schattler: A synthesis of optimal control for a model of tumour growth, Proc. 44th IEEE CDC and ECC 2005, Seville 2005, pp. 934-939.

In this study we have shown how using quite simple models we can analyze and design therapy protocols of combined antiangiogenic and radiation therapy. This type of cancer treatment is still in experimental and clinical trials. The results are promising however still the knowledge of the processes behind the evolution of cancer vascular network, the equilibrium between stimulated and inhibitory factors, different forms of antiangiogenic therapy, its side effects and the results of combined used of different treatment modalities is far from being complete. We hope that our results may help in the progress in this field. The interesting question is if and how the asymptotic strategy can be compared with the optimal one in the finite horizon in terms of morbidity and cost. Till now we have no answer to this question.

[9]..J. La Salle and S. Lefschetz: Stability by Liapunov’s Direct Method, Academic Press, New York, 1961. [10] H.D. Thames,and J.H. Hendry: Fractionation in radiotherapy, Taylor&Francis, London, 1987. [11] L.S. Pontryagin, V.G. Boltyanski, R.V. Gamkrelidze, and Y.F. Mishchenko: The Mathematical Theory of Optimal Processes, Mac Millan, New York, 1964.

117