Adv Syst Sci Appl 2020; 04:83–104 Published online at https://ijassa.ipu.ru. On Optimal Therapy Protocols in the Mathematical Model of Prostate Cancer Progression Marie-Christin Litzinger1, Y. Todorov1, Miriam Föller-Nord1, Manitjaiswal K. Chaudhary2, Alexander S. Bratus3,4∗ 1University of Applied Sciences, Paul-Wittsack-Str. 10, 68163, Mannheim, Germany 2Nikolsky Mathematical Institute, Rudn University, Miklukho-Maklaya str.6, Moscow, 117198, Russia 3Russian University of Transport, Obraztsova 15, Moscow, 127994, Russia 4Moscow Center of Fundamental and Applied Mathematics, MSU, GSP-1, Leninskie Gory, Moscow, 119991, Russia Abstract: A mathematical model for the treatment of prostate cancer regarding androgen ablation is considered. The effect of the therapy on healthy cells is taken into account and modelled as a second order state constraint. A mathematical overview about possible therapy strategies is presented. Keywords: prostate cancer, androgen ablation, therapy strategies, optimal control, state constraint 1. INTRODUCTION Regarding prostate cancer progression and development testosterone and its role in the biological mechanisms of epithelial prostate cells are of particular interest for developing therapy strategies. Therefore the main steps of testosterone procession are now outlined briefly. After entering the cells testosterone is either able to bind to androgen receptors (AR) straight away or it is first converted to dihydrotestosterone (DHT). In contrast to testosterone, DHT and AR form a more stable complex. After being phosphorylated and dimerized that complex binds to the DNA, where it causes a higher transcription rate of with proliferation, survival and differentiation associated genes [11]- [10]. Consequently, there are several points in the mechanism where the therapy can intervene, for example through inhibiting the testosterone production or blocking androgen receptors [2]. After some time prostate cancer cells are able to grow even in the absence of androgen. Due to the following therapy resistance an intermittent androgen therapy can be applied [9, 19]. It is important that healthy and cancer cells both produce prostate-specific antigene (PSA). In noncancerous cells only a limited amount of PSA enters into the bloodstream due to natural barriers. In the presence of cancer these barriers break down and more PSA leaks into the blood. The PSA-level in the blood serum therefore acts as an indicator for the stadium and the presence of prostate cancer [11]. Due to the complex mechanisms of cancer the cooperation of different scientific areas is necessary and mentioned in [3] and [6]. Since the last fifty years the application of ∗Corresponding author: alexander.bratus@yandex.ru 84 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS mathematical approaches in order to analyze the dynamics of cancer is therefore common and discussed in [3]. Furthermore mathematical models and methods for the analysis of the dynamics of tumor-immune interaction are considered to evolve research perspectives [4]. However, not only the competition between cancer and healthy cells, but also the toxicity of the therapy, which damages healthy cells, are negative side effects. These are essential for developing effective medicine and therapy strategies. Consequently, the mathematical theory of optimal control processes is needed to develop successful therapy protocols [17]- [22]. Due to the adverse effects caused by the ablation of androgen, for example the decrease of quality of life [5], it is convenient to introduce a limitation of the healthy epithelial prostate cells during androgen ablation therapy. In the present work three approaches – an alternative, suboptimal and optimal one – in the mathematical model for prostate cancer progression in response to androgen ablation therapy are introduced. In all three approaches the side effects of the androgen ablation therapy are considered as a second order pure state constraint on the amount of the healthy epithelial prostate cells. The goal is to minimize the amount of cancer cells while keeping the amount of healthy cells above a given critical limit. The resistance of prostate cancer cells during the androgen ablation therapy is also taken into account. In the following a brief description of each therapy strategy is given. The alternative therapy strategy is based on the dynamical analysis of the system and its critical points. This approach circumvents the classical optimal control theory. In combination with a numerical improvement algorithm it provides feasible results which are similar to the optimal one. Pontraygin’s maximum principle without state constraints is used to develop a therapy protocol at which the conditions of the states on the constraint boundary are analyzed. The considered approach always guarantees a feasible therapy process if there exists such possible one. The results are used to obtain the so called suboptimal therapy strategy. An analytical determination of an optimal therapy strategy in the presence of state constraints of order two is considered by means of analytical and numerical techniques. Therefore the software tool BOCOP is used. The results of the presented approaches are compared and discussed. 2. MATHEMATICAL MODEL OF PROSTATE CANCER In this article a mathematical model of prostate cancer based on [11] is considered. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION 85 dE(t) dt = µEE(t) ( 1− E(t) ηE ) − εE(t)(M(t) +N(t)) ηE − δEE(t) dN(t) dt = (1− αm)µNN(t) ( 1− N(t) +M(t) η ) − δNN(t) dM(t) dt = (µMM(t) + αmµMN(t)) ( 1− N(t) +M(t) η ) − δMM(t) dR(t) dt = αR − λRR(t)− kDf R(t)D(t) + kDr A(t)− kTf R(t)T (t) + kTr A(t) dD(t) dt = βT T (t) kT + T (t) − λDD(t)− kDf R(t)D(t) + kDr A(t) dT (t) dt = f(Ts)− λTT (t)− βT (t) T (t) kT + T (t) − kTf R(t)T (t) + kTr At(t) dA(t) dt = −λAA(t) + kDf R(t)D(t)− kDr A(t) dAt(t) dt = −λAtAt(t) + kTf R(t)T (t)− kTr At(t) dP (t) dt = βp(A(t) + At(t) + ϕ)− αPP (t)Efrac − γpP (t) (M(t) +N(t))2 kP + (M(t) +N(t)) − λPP (t) t ∈ [0, T ], E(0), N(0), M(0), R(0), D(0), T (0), A(0), At(0), P (0) ∈ R+, (2.1) where Efrac = vc E PV ol , PV ol = αvt dv bdvv + tdv . In (2.1) it is assumed that the healthy cells E have the proliferation rate µE , the death rate δE and a limit ηE . The magnitude of the competition between healthy and cancer cells is modelled by means of the parameter ε. The androgen-dependent cancer cells N proliferate at rate µN , decease with rate δN and have the limit ηN . With the probability αm ∈ (0, 1) cancer cells N mutate to castration- resistant cancer cells M where the parameters µM and δM represent the proliferation and dead rate respectively. Because of the character of this article the remaining state variables of the dynamic system (2.1) which represent specific hormone levels in the organism are only mentioned: T - testosterone concentration, P - tissue PSA concentration, R - free androgen receptor concentration, D - dihydrotestosterone concentration, A - dihydrotestosterone-activated androgen receptor concentration,At - testosterone-activated androgen receptor concentration. In this article all cell types are estimated in billions. Detailed information, the biological background of the variables and the associated parameters are given in [?]. Assuming µN = µM = µ, δN = δM = δ, N(t) and M(t) can be substituted in order to represent the total number of the prostate cancer cells. From the addition N(t) +M(t) = Z(t) follows dZ(t) dt = µZ(t) ( 1− Z(t) η ) − δZ(t). (2.2) The following figures illustrate the behaviour of healthy and cancer cells as well as the PSA-level in the absence of therapy. The parameters of the model which were used in this Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 86 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS article describe a scenario in which the amount of healthy cells is not strongly influenced by cancer cells. The set of parameters and their description is given (see A.1). 0 50 100 150 non-dimensional time 0 50 100 150 200 250 300 Fig. 2.1. Time response of healthy E(t) and cancer cells Z(t) 0 50 100 150 non-dimensional time 0 1 2 3 4 5 6 7 8 9 10 Fig. 2.2. Time response of the PSA function P (t) Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION 87 As mentioned before intermittent therapy is one way to slow down the development of therapy resistance. The following figures show the dynamics of the system when such therapy is applied. 0 50 100 150 non-dimensional time 0 1 2 3 4 5 Fig. 2.3. Amount of the administered drug u(t) 0 50 100 150 non-dimensional time 0 50 100 150 200 250 300 Fig. 2.4. Time response of healthy E(t) and cancer cells Z(t) Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 88 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS 0 50 100 150 non-dimensional time 0 1 2 3 4 5 6 7 8 9 10 Fig. 2.5. Time response of the PSA function P (t) In the following sections different therapy strategies are determined and used to minimize cancer cells without violating the constraint boundary on healthy cells. 3. AN ALTERNATIVE THERAPY STRATEGY Consider the optimal control problem after the substitution (2.2) min u J = ∫ T 0 Z(t)dt→ min (3.3) Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION 89 dE(t) dt = µEE(t) ( 1− E(t) ηE ) − εE(t)(M(t) +N(t)) ηE − δEE(t) dZ(t) dt = µ(h)Z(t) ( 1− Z(t) η ) − δZ(t) dR(t) dt = αR − λRR(t)− kDf R(t)D(t) + kDr A(t)− kTf R(t)T (t) + kTr A(t) dD(t) dt = βT T (t) kT + T (t) − λDD(t)− kDf R(t)D(t) + kDr A(t) dT (t) dt = f(Ts)− λTT (t)− βT (t) T (t) kT + T (t) − kTf R(t)T (t) + kTr At(t) dA(t) dt = −λAA(t) + kDf R(t)D(t)− kDr A(t) dAt(t) dt = −λAtAt(t) + kTf R(t)T (t)− kTr At(t) dP (t) dt = βp(A(t) + At(t) + ϕ)− αPP (t)Efrac − γpP (t) (M(t) +N(t))2 kP + (M(t) +N(t)) − λPP (t) dh(t) dt = −γh(t) + u(t) 0 ≤ u(t) ≤ umax, t ∈ [0, T ], T -fixed E(0), Z(0), R(0), D(0), T (0), A(0), At(0), P (0), h(0) ∈ R+, (3.4) and the state constraint on the amount of healthy cells g(t) = E(t)− EC ≥ 0 ∀t ∈ [0, T ]. (3.5) The last equation in (3.4) describes the pharmacokinetic and the functions µ(h) = µ0 ( 1− kh h+ 1 ) and µE(h) = µ0 ( 1− kEh h+ 1 ) the pharamcodynamic of the administered medicine. The parameters γ, k, kE, µ0, µ0 E ∈ R+, where γ stands for the dissipation rate of the administered drug respectively, k and kE describe the effect of the drug on cancer and healthy cells and µ0 and µ0 E are the replication rates of the two cell populations. The equation (3.5) has to guarantee that the amount of the healthy cells does not undercut a given critical limit Ec during the therapy process. Because of the mathematical complexity of the optimal control problem (3.3)-(3.5) an alternative approach, which will further be called alternative therapy strategy, is considered at first. In order to analyze whether small changes of the initial conditions of the considered system lead to a different behaviour of the corresponding trajectories its critical points are analyzed (see A.2). Due to the practical irrelevance of the trivial critical point (Ē, Z̄, h̄) = (0, 0) only Z̄ = ( 1− δ µ(h) ) η, Ē = ηE − ε(µ(h)− δ)η µE(h)µ(h) − δEηE µE(h) , h̄ = u γ (3.6) is taken into account. From the dynamical analysis of the considered system (see A.3) it is obvious, that for realistic parameters of the model the critical point (3.6) is a stable node. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 90 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS Assuming now that u is a constant parameter such that 0 ≤ u ≤ umax we consider the objective function J̃(u) := J̃(Ē, Z̄) = Z̄ + κ(Ē − Ec)2 → min, (3.7) where with the scale of the weighting factor κ ∈ R+ the state constraint on the amount of healthy cells is taken into account. It follows from the implicit function theorem that Ē and Z̄ can be considered as functions of u and consequently the minimization of the functional (3.7) can be approached as a mathematical programming problem. Due to the theorem of Weierstrass a solution of this problem always exists. Assuming that u = ū ∈ [0, umax] is the value for which the minimum of the functional (3.7) is achieved. Now consider the function u(t) = umax, 0 ≤ t ≤ t̂. Here umax and t̂ are two constants such that the solution of the last equation in (3.4) with initial condition h(0) = 0 reaches the value h = ū γ at the moment t = t̂. Since h(t) = umax γh (1− e−γt) (3.8) the condition h(t̂) = ū γ is fulfilled iff t̂ = 1 γ ln ( 1− ū umax ) . (3.9) The last equality provides the required time of the function h(t) to reach the value of ū γ if the control function u(t) = umax in the time interval 0 ≤ t ≤ t̂ , which is called intensive therapy time. In the remaining time interval t̂ < t ≤ T the control function u(t) has to keep the value ū γ which minimizes the objective function J̃(u). This therapy interval time is called relaxation therapy time. The determined strategy is called alternative therapy strategy and can be defined as u(t) = umax, 0 ≤ t ≤ t̂, ū γ , t̂ < t ≤ T. (3.10) According to the considerations above the alternative control strategy consists out of two stages: the stage of the intensive therapy and the stage of relaxation. In order to avoid disadvantages in the approach considered above, the following modification, presented in [22], will be used. This requires the introduction of the shifting variable s ∈ R+ in the penalty term of the objective function (3.7) J̃(u) := J̃(Ē, Z̄) = Z̄ + κ(Ē − (Ec + s))2 → min. (3.11) In the numerical procedure some iterations are usually necessary to obtain a good result and the new value of the shifting variable after every iteration has to be calculated as follows: s(q+1) = s(q) + max{Ec − E(t)| E(t) < Ec, t ∈ (0, T )}, where q is the number of the last iteration. If the maximal number of iteration qmax or the desired accuracy θ ≥ |s(q) − s(q−1)| is reached the procedure has to be terminated. The following numerical results illustrate the time response of healthy and cancer cells as well as the PSA-level and the control function using the alternative therapy strategy. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION 91 0 50 100 150 non-dimensional time 44 46 48 50 52 54 56 58 Fig. 3.6. Time response of healthy cells E(t) using the alternative therapy strategy 0 50 100 150 non-dimensional time 115 120 125 130 135 140 145 150 Fig. 3.7. Time response of cancer cells Z(t) using the alternative therapy strategy Cancer and healthy cells decrease (Fig. 3.6), (Fig. 3.7). Healthy cells reach the constraint boundary without violating it. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 92 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS 0 50 100 150 non-dimensional time 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Fig. 3.8. Amount of the administered drug u(t) using the alternative therapy strategy 0 50 100 150 non-dimensional time 0 1 2 3 4 5 6 7 8 9 10 Fig. 3.9. Time response of the PSA function P (t) using the alternative therapy strategy The maximal dose of therapy is used for a short amount of time and then a constant amount is administered (Fig. 3.8). Due to the reduction of cancer cells less PSA is produced in the tissue and can consequently be found in the blood serum (Fig. 3.9). Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION 93 4. SUBOPTIMAL THERAPY STRATEGY In this section a suboptimal therapy strategy for the problem (3.3)-(3.5) is introduced. In order to minimize the objective function (3.3) it is natural to apply the maximal admissible amount of the drug so long as possible. However, it has to be guaranteed that the amount of healthy cells does not violate the constraint (3.5). I.e. that the time derivation Ė on the constraint boundary has to be non negative. From Ė(tj) ≥ 0 follows µE(h)Ec ( 1− Ec ηE ) − εEcZ(tj) ηE − δEEc ≥ 0 ⇒ µ0 E ( 1− kEh(tj) h(tj) + 1 ) ≥ εZ(tj) + δEηE ηE − Ec ⇒ kEh(tj) ≤ ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) (h(tj) + 1) (4.12) ⇒ h(tj) ( kE − ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E )) ≤ ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) ⇒ h(tj) ≤ ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) kE − ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) , where tj := {t ∈ (0, T )| E(t) = Ec}. Due to the continuity of E the strong inequality in the last expression can not appear. The last expression in (4.12) gives the necessary condition for h on the constraint boundary, which guarantees an admissible process. Therefore the following analysis is needed. From the last equation of the system (3.4) holds for u(t) = umax h(tj) = e−γtj ∫ tj 0 eγtumaxdt = umax γ ( 1− e−γtj ) . (4.13) Thus if the maximal amount of the drug, u(t) = umax, t ∈ [0, tj), is administered until the critical boundary of the healthy cells is reached, the following relation has to be fulfilled ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) kE − ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) = umax γ (1− e−γtj). (4.14) Note, that this is a special case which can only appear with a given set of parameters of the model and initial conditions. In this case the therapy strategy can be changed when the amount of the healthy cells reaches the critical value without violating the therapy process. The therapy strategy stated below is in this case optimal (see A.4). Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 94 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS u(t) =  umax, 0 ≤ t ≤ tj, γ 1− εZ(t) + δEηE (ηE − Ec)µ0 E  kE− 1− εZ(t) + δEηE (ηE − Ec)µ0 E  , tj < t ≤ T. (4.15) In general the change of the administered amount of the drug can not induce the change of the slope of the development of the healthy cells immediately. This is caused by the inertia of the dynamic system. Therefore investigations of the more realistic cases have to be done. The case in which the maximal admissible amount of the drug has to be changed before the critical boundary of the healthy cells is reached, is considered. I.e.,( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) kE − ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) < umax γ (1− e−γtj). (4.16) The admissible process can be guaranteed if the maximal admissible amount of the drug is reduced at the point of time t′ < tj so that h(t′) = ( 1− εZ(tj(umax)) + δEηE (ηE − Ec)µ0 e ) kE − ( 1− εZ(tj(umax)) + δEηE (ηE − Ec)µ0 E ) holds, where Z(tj(umax)) is the amount of the cancer cells on the constraint boundary using the maximal dose of medicine (u(t) = umax, t ∈ [0, tj]). Thus the suboptimal therapy strategy can be given by u(t) =  umax, 0 ≤ t ≤ t′, γ 1− εZ(tj(umax)) + δEηE (ηE − E(t))µ0 E  kE− 1− εZ(tj(umax)) + δEηE (ηE − E(t))µ0 E  , t′ < t ≤ T. (4.17) The following figures present the results of the suboptimal control approach. Cancer and healthy cells decrease and the healthy cells do not reach the constraint boundary (Fig. 4.10), (Fig. 4.11). Since cancer cells are decreasing continuously the quantity of PSA in the tissue and therefore in the blood serum also decreases (Fig. 4.13). 5. OPTIMAL THERAPY STRATEGY In this section the therapy strategy related to the optimal control problem (3.3)-(3.5) is investigated with Pontryagin’s maximum principle, which is the necessary condition for an optimal process. The restriction on the minimal admissible amount of healthy cells (3.5) has to be interpreted as a second order state constraint. Because of the high complexity of optimal Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION 95 0 50 100 150 non-dimensional time 44 46 48 50 52 54 56 58 Fig. 4.10. Time response of healthy cells E(t) using the suboptimal therapy strategy 0 50 100 150 non-dimensional time 120 125 130 135 140 145 150 Fig. 4.11. Time response of cancer cells Z(t) using the suboptimal therapy strategy control problems with higher order state constraints the introduced approach is a composition of analytical investigations and numerical results. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 96 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS 0 50 100 150 non-dimensional time 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Fig. 4.12. Time response of the quantity of the u(t) using the suboptimal therapy strategy 0 50 100 150 non-dimensional time 0 1 2 3 4 5 6 7 8 9 10 Fig. 4.13. Time response of the PSA function P (t) using the suboptimal therapy strategy Because of the nature of the maximum principle it is essential to primarily investigate the optimal therapy strategy without state constraints. The optimality of this process till the state variable E reaches the constraint boundary is determined. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION 97 The investigations given in A.3 prove the following theorem. Theorem 5.1: The optimal therapy strategy in the absence of the state constraint (3.5) is given as u∗(t) = umax ∀t ∈ [0, T ]. From the practical point of view this result is equivalent to the situation in which the test results of a patient are acceptable during the therapy. Therefore the therapy has to be executed with the maximal admissible amount of the drug in order to eradicate the cancer cells. In the opposite case the amount of healthy cells E reaches the critical value Ec. Because of the complexity the strict mathematical investigation of the constraint problem is not considered. Instead of that numerical results using the software tool BOCOP are presented. The numerical results show the development of healthy cells, cancer cells and the optimal control. 0 50 100 150 non-dimensional time 40 42 44 46 48 50 52 54 56 58 60 Fig. 5.14. Time response of healthy cells E(t) using the optimal therapy strategy Healthy and cancer cells decrease (Fig. 5.14), (Fig. 5.15). Healthy cellsE(t) then progress along the boundary. After a period where the maximal amount of drug is administered follows a short relaxation phase. Then the maximal admissible amount of medicine which minimizes the objective function is applied. The following table presents the values of the objective function, the healthy cells and cancer cells at the end of the therapy process (T = 150). Note that in the results achieved by the alternative therapy strategy the modification (3.11) is used. Therefore the results of the alternative strategy are better than the ones obtained with the suboptimal one. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 98 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS 0 50 100 150 non-dimensional time 100 110 120 130 140 150 160 Fig. 5.15. Time response of the cancer cells Z(t) using the optimal therapy strategy 0 50 100 150 non-dimensional time 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Fig. 5.16. Amount of the administered drug u(t) using the optimal therapy strategy 6. SUMMARY A simplified mathematical model for the therapy of prostate cancer regarding androgen ablation is considered. The effect of the therapy on cancer and healthy cells is assumed to Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION 99 Therapy strategy J(T ) E(T ) Z(T ) Alternative therapy strategy 1.988 · 104 46.222 117.88 Suboptimal therapy strategy 2.029 · 104 47.414 123.14 Optimal therapy strategy 1.195 · 104 46.000 114.98 Table 5.1. Comparison of the terminal values of the three considered approaches reduce the birth rate of the populations of the considered cell types. An optimization problem with a constraint on the amount of healthy cells during the therapy process is introduced. Three therapy strategy approaches of the problem are illustrated and solved by combinating analytical and numerical procedures. The presented results of the therapy strategies using the suboptimal and alternative control demonstrate that the mathematical challenges of the classical optimal control problem with state constraints can be circumvented and the approaches can be used for medicine. Moreover the suggested mathematical model produces possibilities for further investigations regarding the physical condition of a patient which can be modelled by means of different parameters of the model. From the mathematical point of view the considered problem covers theoretical investigations on the necessary and sufficient conditions of optimal control problems with state constraints. ACKNOWLEDGEMENTS This research is supported by the Russian Scientific Foundation Grant 19-11-00008 and Grant RFBI 20-04-60157 and Ministry of Science of Higher education grant 075-15-2019- 1621. REFERENCES 1. Agus, D.B., et al. (1990), Prostate cancer cell cycle regulators: Response to androgen withdrawal and development of androgen independence, J Natl Cancer Inst, 91, 1869– 1876. 2. Balk, S.P. (2002), Androgen receptor as a target in androgen-independent prostate cancer, Urology, 60, 132–138. 3. Byrne, H.M. (2010), Dissecting cancer through mathematics: from the cell to the animal model, Nature Reviews Cancer, 10, 221–230. 4. Bellomo, N. & Preziosi, L. (2000), Modelling and mathematical problems related to tumor evolution and its interaction with the immune system, Mathematical and Computer Modelling, 32, 413–452. 5. Donovan, K.A., Walker, L.M., Wassersug, R.J., Thompson, L.M. & Robinson, J.W. (2015), Psychological effects of androgen-deprivation therapy on men with prostate cancer and their partners, Cancer, 121, 4286–4299. 6. Gatenby, R.A. and Maini, P.K. (2003), Mathematical oncology: Cancer summed up, Nature, 421, page 321. 7. Gann, P.H., Hennekens C.H., Ma J., Longcope, C., Stampfer, M.J. (1996), Prospective study of sex hormone levels and risk of prostate cancer, J Natl Cancer Inst, 88, 1118– 1126. 8. Heinlein, C.A., Chang, C. (2004), Androgen receptor in prostate cancer, Endocr Rev, 25, 276–308. 9. Hirata, Y., Tanaka, G., Bruchovsky, N. and Aihara, K. (2012), Mathematically modelling and controlling prostate cancer under intermittent hormone therapy, Asian Journal of Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 100 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS Andrology, 270–277. 10. Hussain, M., et al. (2009), Prostate-specific antigen progression predicts overall survival in patients with metastatic prostate cancer: Data from the Southwest Oncology Group Trials 9346 (Intergroup Study 0162) and 9916, J Clin Oncol, 27, 2450–2456. 11. Jain, H.V., Clinton, S.K., Bhinder, A. & Friedman, A. (2011), Mathematical modelling of prostate cancer progression in response to androgen ablation therapy, Proc. Natl Acad. Sci. USA, 108, 19701–19706. 12. Ledzewicz, U., Schättler, H. (2006), Application of optimal control to a system describing tumor anti-angiogenesis, Proceedings of the 17th International Symposium on Mathematical Theory of Networks and Systems (MTNS), Kyoto, Japan, 478–484. 13. Ledzewicz, U., Schättler, H. (2005), A synthesis of optimal controls for a model of tumor growth under angiogenic inhibitors, Proceedings of the 44th IEEE Conference on Decision and Control, 945–950. 14. Miyamoto, H., Messing, E.M., Chang, C. (2004), Androgen deprivation therapy for prostate cancer: Current status and future prospects, Prostate, 61, 332–353. 15. Murray, J.M. (1990), Optimal control for a cancer chemotherapy problem with general growth and loss functions, Math Biosci., 98(2), 273–287. 16. Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V. & Mishchenko EF (1962), The Mathematical Theory of Optimal Processes, Interscience Publishers. 17. Rubinow S.I., Lebowitz J.L. (1976), A mathematical model of the chemotherapeutic treatment of acute myeloblastic leukemia, Biophys. J., 1257–1271. 18. Seidenfeld J., et al. (2000), Single-therapy androgen suppression in men with advanced prostate cancer: A systematic review and meta-analysis, Ann Intern Med, 132, 566–577. 19. Suzuki, Y., Sakai, D., Nomura, T., Hirata, Y. and Aihara, K. (2014), A new protocol for intermittent androgen suppression therapy of prostate cancer with unstable saddle-point dynamics, Journal of Theoretical Biology, 350, 1–16. 20. Swanson, K.R., True, L.D., Lin, D.W., Buhler, K.R., Vessella, R. and Murray, J.D. (2001), A Quantitative Model for the Dynamics of Serum Prostate-Specific Antigen as a Marker for Cancerous Growth: An Explanation for a Medical Anomaly, The American Journal of Pathology, 158(6), 2195–2199. 21. Swan, G.W. (1990), Role of optimal control theory in cancer chemotherapy, Math Biosci., 101(2), 237–284. 22. Todorov, Y., Nuernberg, F. (2014), Optimal therapy protocols in the mathematical model of acute leukemia with several phase constraints, OCAM, 35(5), 559–574. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION101 A. APPENDIXES A.1. Model parameter values The following values were used for the model parameters (taken from [11]). Parameter Value Parameter Value Parameter Value ηE 110.24 δE 0.0081 δ 0.00405 γ 0.1 µE0 0.22 µ0 0.1 kE 0.889 k 0.99 aTs 1.00 bTs 22.23 cTs 0.23 kTf 0.14 kTr 0.07 αP 0.001515 ε 0.009 λD log(2) 9 λP log(2) 12.3 λAt log(2) 3 βM 0.003 βN 0.002 βT 4.569 pN 10 pM 10 KP 7 γP 0.6 104 kDf 0.018 kDr 0.053 100 KT 0.104 AEt 1.62 λR log(2) 3 λA log(2) 3 λT log(2) 3 aV 18.397 bV 67.753 cV 28.4 dV 12 Vc 5.56 106 r0 E 90 0.002622 β Table A.2. Parameters of the model More detailed information is given in [11]. A.2. Dynamical Analysis The motivation for the construction of the alternative control strategy rests upon the asymptotically behaviour of the considered system (3.4). A relocation of an asymptotic stable critical point in order to minimize the functional (3.3) provides the alternative therapy strategy. Therefore the type of the nontrivial critical point given in (3.6) is analyzed. Define the Jacobian J = ∂f ∂x =  ∂E(t) ∂E ∂E(t) ∂Z ∂Z(t) ∂E ∂Z(t) ∂Z  = µE(h)− 2EµE(h) ηE − εZ ηE − δE −εE ηE 0 µ(h)− 2Zµ(h) η − δ  At the considered critical point (3.6) holds J |(Ē,Z̄) = −µE(h̄) + ĒµE(h̄) ηE + δE −εĒ ηE 0 −µ(h̄) + δ  . Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 102 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS In order to analyze the behaviour of the considered critical point it is necessary to analyze the sign of the eigenvalues of the Jacobian. I.e., to find the roots of the polynomial defined by det(J − λĒ) = ( − µE(h̄) + ĒµE(h̄) ηE + δE − λ )( − µ(h̄) + δ − λ ) . The eigenvalues λ1, λ2 of J are λ1 = −µ(h̄) + δ < 0 and λ2 = −µE(h̄) + ĒµE(h̄) ηE + δE. For the given set of parameters (see A.1) the critical point (Ē, Z̄) is always stable. A.3. Optimal Control without state constraint Constructing the Pontryagin-Hamilton function [16] to the problem (3), (4) H := H(E,Z,R,D, T,A,At, P, h, ψ0, ψ1, ψ2, ψ3, ψ4, ψ5, ψ6, ψ7, ψ8, ψ9, u) = + ψ0(Z) + ψ1 ( µEE ( 1− E ηE ) − εEZ ηE − δEE ) + ψ2 ( µZ(t) ( 1− Z(t) η ) − δZ(t) ) + + ψ3 ( αR − λRR− kDf RD + kDr A− kTf RT + kTr A ) + + ψ4 ( βT T kT + T − λDD − kDf RD + kDr A ) + + ψ5 ( f(Ts)− λTT − βT T kT + T − kTf RT + kTr At ) + + ψ6 ( − λAA+ kDf RD − kDr A) + ψ7(−λtAt + kTf RT − kTr At ) + + ψ7 ( −λAtAt(t) + kTf R(t)T (t)− kTr At(t) ) + + ψ8 ( βp(A+ At + ϕ)− αPPEfrac − γpP (N +M)2 kP +N +M − λPP ) + ψ9 ( − γh+ u ) the optimal control law is given as u∗(t) = umax, ψ9(t) > 0, 0, ψ9(t) < 0, us ∈ (0, umax), ψ9(t) = 0 ∀t ∈ (t1, t2), (t1, t2) ⊂ [0, T ], where ψ9(t) is the so called switching function. Note, that the optimal control function u∗(t) maximizes H in all points t ∈ [0, T ]. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) ON OPTIMAL THERAPY PROTOCOLS IN THE MATH. MODEL OF PROSTATE CANCER PROGRESSION103 From Pontryagin’s maximum principle the system of adjoint variables is given by ψ̇1(t) = −∂H ∂E = −ψ1 ( +µE(h(t))− 2E(t)µE(h(t)) ηE − εZ(t) ηE − δE ) , (A.18) ψ̇2(t) = −∂H ∂Z = −ψ0 + ψ1 E(t)ε ηE + ψ2 ( −µ(h(t)) + 2µZ(t) η + δ ) , (A.19) ψ̇9(t) = −∂H ∂h = −ψ1 dµE(h(t)) dh(t) E(t) ( 1− E(t) ηE ) − ψ2 dµ(h(t)) dh(t) µZ(t) ( 1− Z(t) η ) + +ψ9γ, ψ1(T ) = ψ2(T ) = ψ9(T ) = 0. Lemma A.1: In the optimal control problem (3.4) holds ψ0 = −1. Proof From ψi(T ) = 0, 1 ≤ i ≤ 9, and the condition (ψ0, ψ1, ..., ψ9) 6= 0 follows that ψ0 = −1. From Lemma A.1 follows ψ1(t) = 0 ∀t ∈ [0, T ], ψ2(t) = ew2(t) (∫ t 0 e−w2(s)ds− ∫ T 0 e−w2(t)dt ) , ψ9(t) = eγt[− ∫ t 0 e−γsψ2(s) dµ(h(s)) dh(s) Z(s) ( 1− Z(s) η ) ds+ + ∫ T 0 e−γtψ2(t) dµ(h(t)) dh(t) Z(t) ( 1− Z(t) η ) dt], where w2(t) = ∫ t 0 ( −µ(h(s)) + 2µ(h(s))Z(s) η + δ ) ds. ψ9(t) depends on ψ2(t) and therefore this function is analysed for the mathematical solution of this problem. Note, that ψ9(t) also depends on ψ1(t) but as shown before ψ1(t) ≡ 0. Proof of theorem 5.1 Lemma A.2: ψ2(t) is a negative monotone increasing function on [0, T ]. Proof In order to analyse the behaviour of ψ2(t) the auxiliary function ψ̃2 = e−w2(t)ψ2(t) which has the same sign and zeros as ψ2(t) is used. It holds ψ̃2(t) < 0, t ∈ [0, T ) due to T∫ 0 e−w2(t)dt > t∫ 0 e−w2(s)ds ∀t ∈ [0, T ) and ψ̃2(T ) = 0. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) 104 M.-C. LITZINGER, Y. TODOROV, M. FÖLLER-NORD, M.K. CHAUDHARY, A.S. BRATUS Because t∫ 0 ew2(s)ds is a strictly increasing function ψ̃2(t) has no zeros on [0, T ) and therefore ψ2(t) is a negative monotone increasing function with ψ2(t) ≤ 0 ∀t ∈ [0, T ]. Lemma A.3: ψ9(t) is a positive decreasing function on [0, T ]. Proof From ψ9(t) = eγt ( ψ90 − ∫ t 0 e−γsψ2(s) µ(h(s)) dh(s) Z(s) ( 1− Z(s) η ) ds ) and ψ9(T ) = 0 ψ90 = ∫ T 0 e−γtψ2(t) µ(h(t)) dh(t) Z(t) ( 1− Z(t) η ) dt From the definition of the given problem follows Z(t) ≥ 0 ∀t ∈ [0, T ] and( 1− Z(t) η ) ∈ (0, 1). Moreover it holds ψ2 ≤ 0 and dµ(h) dh < 0 ∀t ∈ [0, T ]. Therefore it follows∫ T 0 e−γtψ2(t) dµ(h(t)) dh(t) Z(t) ( 1− Z(t) η ) dt > ∫ t 0 e−γsψ2(s) dµ(h(s)) dh(s) Z(s) ( 1− Z(s) η ) ds ∀t ∈ [0, T ]. Thus ψ9(t) > 0 ∀t ∈ [0, T ) and u∗(t) = umax. A.4. Proof of the statement Proof From Theorem 5.1 follows u∗(t) = umax for t ∈ [0, tj). Because of the minimization of the objective function (3.3) the maximal admissible value of u on the constraint boundary has to be chosen. The limit of this value is given by u(t) = ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) kE − ( 1− εZ(tj) + δEηE (ηE − Ec)µ0 E ) , t ∈ [tj, T ]. Copyright c© 2020 ASSA. Adv Syst Sci Appl (2020) Introduction Mathematical model of prostate cancer An alternative therapy strategy Suboptimal therapy strategy Optimal therapy strategy Summary Appendixes Model parameter values Dynamical Analysis Optimal Control without state constraint Proof of the statement