Optimal control to reach eco-evolutionary stability in metastatic castrate-resistant prostate cancer

Abstract

In the absence of curative therapies, treatment of metastatic castrate-resistant prostate cancer (mCRPC) using currently available drugs can be improved by integrating evolutionary principles that govern proliferation of resistant subpopulations into current treatment protocols. Here we develop what is coined as an 'evolutionary stable therapy', within the context of the mathematical model that has been used to inform the first adaptive therapy clinical trial of mCRPC. The objective of this therapy is to maintain a stable polymorphic tumor heterogeneity of sensitive and resistant cells to therapy in order to prolong treatment efficacy and progression free survival. Optimal control analysis shows that an increasing dose titration protocol, a very common clinical dosing process, can achieve tumor stabilization for a wide range of potential initial tumor compositions and volumes. Furthermore, larger tumor volumes may counter intuitively be more likely to be stabilized if sensitive cells dominate the tumor composition at time of initial treatment, suggesting a delay of initial treatment could prove beneficial. While it remains uncertain if metastatic disease in humans has the properties that allow it to be truly stabilized, the benefits of a dose titration protocol warrant additional pre-clinical and clinical investigations.

Fig 1. Population densities for stable equilibria.

Population densities $x_{T^{+}}^{*}$, $x_{T^P}^{*}$ and $x_{T^{-}}^{*},$ corresponding to stable equilibria for Λ ∈ [0, 1] and for α₃₂ = 2.0. The gray highlighted regions show the stable equilibria that are within the patient viability constraint (5). The yellow highlighted points represent two specific stable equilibria chosen for further analysis.

Fig 2. Initial tumor compositions for Forwards Backwards Sweep analysis.

100 randomly selected initial tumor compositions used in the Forwards Backwards Sweep optimal control analysis. All initial total tumor volumes satisfy the patient viability constraint $\sum _{i\in \mathcal{T}}{x}_i\left(t_0\right)\le 9000$.

Fig 3. Forwards Backwards Sweep optimized dosing schedules.

Forward Backwards Sweep results for optimal dosing schedule to arrive at two-species stability point (left panel) and three species stability point (right panel). The mean of all 100 paths is shown with symmetric one standard deviation error bars (dosage values <0 are not possible). Standard error of the mean is on the order of 10⁻³.

Fig 4. System state trajectories under optimal dose schedules.

State trajectories from each of the 100 initial tumor compositions to the two species equilibrium point (left panel) and the three-species equilibrium point (right panel). Paths highlighted in red breach the patient viability constraint before reaching the equilibrium point.

Fig 5. Dose adjustment schematic.

Schematic description of the dose adjustment rules based on the measured total tumor volume shown here, attempting to maintain a total tumor burden at V_a, though the same rules apply to V_b.

Fig 6. Initial tumor compositions for clinical feasible protocols.

10, 000 randomly selected initial tumor compositions used to analyze the clinically feasible protocols. All total tumor volumes satisfy the patient viability constraint ${\displaystyle \sum _{i\in \mathcal{T}}}{x}_i\left(t_0\right)\le 9000$.

Fig 7. Kaplan-Meier survival analysis for clinically feasible treatment strategies.

10, 000 patients were given each of the six clinically feasible treatment strategies. In this way this shows the outcome of 60, 000 simulated patients. Patients that had not yet breached the patient viability constraint by the end of the simulation are labeled as censored.

Fig 8. Maximum tolerable dose state dynamics.

State trajectories for 100 initial tumor compositions (left panel), used for optimization of patients under the maximum tolerated dose protocol. All state trajectories end when the total tumor burden violates the patient viability constraint (5). The two blue dots show the location of the two equilibria (two- and three- species). The right panel shows the population densities of the three cell types in a representative case.

Fig 9. Adaptive therapy state dynamics.

The state trajectory (left panel) and population densities (right panel) of an example patient under the adaptive therapy protocol. The two blue dots show the location of the two stable equilibria (two- and three- species).

Fig 10. Titration protocols resulting in patient survival.

Average titration protocols of patients that did not breach the patient viability constraint within the simulation time. The standard error of the mean (SEM) is on the order of 10⁻³ for all cases, therefore here the error bars show one standard deviation. (A) Λ(t₀) = 1 stabilizing at V_a. (B) Λ(t₀) = 0 stabilizing at V_a. (C) Λ(t₀) = 1 stabilizing at V_b. (D) Λ(t₀) = 0 stabilizing at V_b.

Fig 11. State dynamics of patient undergoing titration protocol.

The dynamics here show an example patient under the initial dose of Λ(t₀) = 0 and attempting to stabilize at a tumor volume equal to V_b = 7000. The state trajectory in the left panel shows the population arriving at the two-species equilibrium. The population densities and abiraterone dose are shown in the right panel.

Fig 12. Initial tumor compositions of surviving patients.

The initial tumor compositions of the patients that did not breach the patient viability constraint within the simulation time for each of the six clinically feasible protocols. Their two dimensional projections are available in S9 in S1 File. (A) Maximum tolerable dose. (B) Adaptive therapy. (C) Λ(t₀) = 1 stabilizing at V_a. (D) Λ(t₀) = 0 stabilizing at V_a. (E) Λ(t₀) = 1 stabilizing at V_b. (F) Λ(t₀) = 0 stabilizing at V_b.

Fig 13. Tumor composition at time of viability constraint breach.

Ternary plots where each red dot indicates the tumor composition of T⁺, T^P, and T⁻ cells at the time a patient reached the viability constraint. (Figures made using [57].) The top highlighted triangle in each figure encompasses the tumor compositions with >80% T⁻ cells. Patients with tumor compositions located in this upper triangle suffered from competitive release of the T⁻ cells. Outside of this upper triangle, treatable cells were still present at the time of viability constraint breach. (A) 100% of patients are located in the top triangle: n = 10000. (B) Adaptive Therapy. 99.90% of patients are located in the top triangle: n = 8861. (C) Λ(t₀) = 1 stabilizing at V_a. 87.13% of patients are located in the top triangle: n = 9088. (D) Λ(t₀) = 0 stabilizing at V_a. 60.30% of patients are located in the top triangle: n = 7580. (E) Λ(t₀) = 1 stabilizing at V_b. 99.97% of patients are located in the top triangle: n = 7961. (F) Λ(t₀) = 0 stabilizing at V_b. 81.86% of patients are located in the top triangle: n = 3445.