Optimal control in combination therapy for heterogeneous cell populations with drug synergies

TL;DR

The paper develops a general optimal control framework for modeling treatment responses in heterogeneous cell populations under multi-drug regimens, explicitly incorporating multiplicative drug effects and drug–drug interactions through a coupled ODE system for pharmacodynamics. It formulates a quadratic cost balancing tumor burden and drug toxicity, derives adjoint equations and optimality conditions, and demonstrates the approach on three canonical cancer-led examples, including a two-population cervical cancer model, a three-population neuroblastoma model, and a proportion-control extension. The results show that time-varying, interaction-aware dosing can outperform naive constant dosing in certain regimes and highlight how parameter context (e.g., proliferation rate, interconversion dynamics, relative drug toxicity) shapes optimal pharmacodynamics, suggesting routes to personalized and risk-aware therapies. The framework is intended to be broadly applicable to multi-drug, multi-population treatments and is complemented by publicly available code for numerical exploration and replication of findings.

Abstract

Cell heterogeneity plays an important role in patient responses to drug treatments. In many cancers, it is associated with poor treatment outcomes. Many modern drug combination therapies aim to exploit cell heterogeneity, but determining how to optimise responses from heterogeneous cell populations while accounting for multi-drug synergies remains a challenge. In this work, we introduce and analyse a general optimal control framework that can be used to model the treatment response of multiple cell populations that are treated with multiple drugs that mutually interact. In this framework, we model the effect of multiple drugs on the cell populations using a system of coupled semi-linear ordinary differential equations and derive general results for the optimal solutions. We then apply this framework to three canonical examples and discuss the wider question of how to relate mathematical optimality to clinically observable outcomes, introducing a systematic approach to propose qualitatively different classes of drug dosing inspired by optimal control.

Figures (6)

• Figure 1: Schematic of the two-population ODE model. Cervical cancer cells are divided into two populations, A and B, according to their phase of the cell cycle. G1 phase (non-proliferative) cells are denoted A, whilst S/G2 phase (proliferative) cells are denoted B. Cells enter the cell cycle at rate $k$. Cisplatin (concentration denoted by $u_c$), reduces the rate of transfer of population B to population A and kills cells in the G1 phase at rate $\delta_A$. Paclitaxel (concentration denoted by $u_p$), reduces the rate of proliferation in population B and leads to cell death at rate $\delta_B$.
• Figure 2: Efficacy of optimal control treatment in two-population model. A: Heat map of relative treatment efficacy, $\eta$ (as defined in Equation \ref{['eq:eta']}) for different system parameters, $\alpha$, and $\beta$. B: Total integrated cost of the drug administered. C: model solutions using optimal control at $\alpha = 0.05$, $\beta = 0.5$ (dashed line) and $\alpha = 0.5$, $\beta = 0.05$ (solid line). D: Optimal controls in the same parameter regimes as panel C. E-G: Optimal controls for $\alpha = \beta = 0.5$ and varying $R_0$ (the entry along the leading diagonal of the cost matrix $R$).
• Figure 3: Efficacy of treatment upon varying drug toxicity. A: Heat map showing the cost incurred by paclitaxel administration with constant cost matrix as in Figure \ref{['fig:ex1EqualToxicity']}. B: Heat map showing cost incurred by cisplatin administration. Black circle indicates the endpoint of the trajectory in parameter space visualised in D. C: Parameterised trajectory through parameter space along curve with significant cisplatin cost variation. D: Heat map showing treatment efficacy along the path in panel C as relative cisplatin toxicity is varied.
• Figure 4: Schematic of the three-population neuroblastoma ODE model. Cells are divided into I-type sympathoblasts, N-type adrenergic cells, and S-type mesenchymal cells. The model includes spontaneous interconversion between all cell types, treatment with RA acting on the I-to-N differentiation pathway, chemotherapy affecting I-type cells, combined trkA-NGF signaling affecting I and N-type cells, and trkB-BDNF signaling affecting N and S-type cells.
• Figure 5: Effect of varying key system parameters on drug cost. A: Diagonal plots showing maximum marginal cost of drug treatment as system parameters are varied. Off-diagonal plots showing pairwise maximum marginal costs as system parameters are varied. B: Phase plane showing regions of parameter space where variation in $\delta$ leads to most change in drug response. The two sensitive drugs are chemotherapy and RA, and they form two distinct regions of the ($\delta_{\text{apop}}$, $\lambda$) space.