Optimizing the Timing and Composition of Therapeutic Phage Cocktails: A Control-Theoretic Approach

Abstract

Viruses that infect bacteria, i.e., bacteriophage or ‘phage,’ are increasingly considered as treatment options for the control and clearance of bacterial infections, particularly as compassionate use therapy for multi-drug-resistant infections. In practice, clinical use of phage often involves the application of multiple therapeutic phage, either together or sequentially. However, the selection and timing of therapeutic phage delivery remains largely ad hoc. In this study, we evaluate principles underlying why careful application of multiple phage (i.e., a ‘cocktail’) might lead to therapeutic success in contrast to the failure of single-strain phage therapy to control an infection. First, we use a nonlinear dynamics model of within-host interactions to show that a combination of fast intra-host phage decay, evolution of phage resistance amongst bacteria, and/or compromised immune response might limit the effectiveness of single-strain phage therapy. To resolve these problems, we combine dynamical modeling of phage, bacteria, and host immune cell populations with control-theoretic principles (via optimal control theory) to devise evolutionarily robust phage cocktails and delivery schedules to control the bacterial populations. Our numerical results suggest that optimal administration of single-strain phage therapy may be sufficient for curative outcomes in immunocompetent patients, but may fail in immunodeficient hosts due to phage resistance. We show that optimized treatment with a two-phage cocktail that includes a counter-resistant phage can restore therapeutic efficacy in immunodeficient hosts.

Appendices

Appendix A: Implementation of Projection Operator $\mathcal {P}_{U}$

Here, we present a closed form of projection operator $\mathcal {P}_{U}$ via a geometric approach (Boyd and Vandenberghe 2004), recall that $u^{*} = \mathcal {P}_{U}(\hat{u})$ in Theorem 1, then we have the following:

$$\begin{aligned} u^{*}&= \hat{u},\ \text {if}\ \hat{u} \in \{ u\ |\ e_{1}^\mathrm{T}u \ge 0,\ e_{2}^\mathrm{T}u \ge 0,\ \mathbb {1}^\mathrm{T}u - 1 \le 0\},\\ u^{*}&= \mathcal {P}_{\mathcal {A}}(\hat{u}),\ \text {if}\ \hat{u} \in \{ u\ |\ \mathbb {1}^\mathrm{T}u - 1\ge 0,\ \widetilde{\mathbb {1}}^\mathrm{T}u + 1 \ge 0,\ \widetilde{\mathbb {1}}^\mathrm{T}u - 1 \le 0 \},\\ u^{*}&= [0,1]^\mathrm{T},\ \text {if}\ \hat{u} \in \{ u\ |\ e_{2}^\mathrm{T}u - 1 \ge 0,\ \widetilde{\mathbb {1}}^\mathrm{T}u - 1 \ge 0 \},\\ u^{*}&= [0,\hat{u}_{2}]^\mathrm{T},\ \text {if}\ \hat{u} \in \{ u\ |\ e_{1}^\mathrm{T}u \le 0,\ e_{2}^\mathrm{T}u \ge 0,\ e_{2}^\mathrm{T}u - 1 \le 0 \},\\ u^{*}&= [0,0]^\mathrm{T},\ \text {if}\ \hat{u} \in \{ u\ |\ e_{1}^\mathrm{T}u \le 0,\ e_{2}^\mathrm{T}u \le 0 \},\\ u^{*}&= [\hat{u}_{1},0]^\mathrm{T},\ \text {if}\ \hat{u} \in \{ u\ |\ e_{1}^\mathrm{T}u \ge 0,\ e_{1}^\mathrm{T}u - 1 \le 0,\ e_{2}^\mathrm{T}u \le 0 \},\\ u^{*}&= [1,0]^\mathrm{T},\ \text {if}\ \hat{u} \in \{ u\ |\ \widetilde{\mathbb {1}}^\mathrm{T}u + 1 \le 0,\ e_{1}^\mathrm{T}u \ge 1 \}, \end{aligned}$$

where $\widetilde{\mathbb {1}} = [-1,1]^\mathrm{T}$ and $\mathcal {A} = \{u\ |\ \mathbb {1}^\mathrm{T}u - 1 = 0\}$. Computing the projection of $\hat{u}$ onto $\mathcal {A}$ is straightforward. The orthogonality principle yields $\hat{u} - \mathcal {P}_{\mathcal {A}}(\hat{u})$ which must be colinear with $\mathbb {1}$, also $\mathcal {P}_{\mathcal {A}}(\hat{u}) \in \mathcal {A}$, i.e., $\hat{u} - \mathcal {P}_{\mathcal {A}}(\hat{u}) = z \mathbb {1}$ and $\mathbb {1}^\mathrm{T}\mathcal {P}_{\mathcal {A}}(\hat{u}) - 1 = 0$. This yields $z = (\mathbb {1}^\mathrm{T}\hat{u} - 1)/(\mathbb {1}^\mathrm{T}\mathbb {1})$ and thus $\mathcal {P}_{\mathcal {A}}(\hat{u}) = \hat{u} - z \mathbb {1}$.

Appendix B: The Optimality System of Optimal Control in Monophage Therapy

Here, we derive the necessary conditions for the optimal control problem (28) via Pontryagin’s maximum principle (PMP) (Pontryagin 2018).

Theorem 2

If $u_{1}^{*}$ is an optimal control that solves problem (28), and $x^{*}(t)$ is the corresponding state trajectory of the initial value problem (24)–(27), and $\lambda ^{*}(t)$ is the costate trajectory of the following terminal value problem:

$$\begin{aligned} \dot{\lambda }^{*}&= -\left( \frac{\partial f}{\partial x}(x^{*})\right) ^\mathrm{T}\lambda ^{*} - \left( \frac{\partial \mathcal {L}}{\partial x}(x^{*}, u^{*})\right) ^\mathrm{T},\quad \lambda ^{*}(t_{f}) = [\theta _{f}, \theta _{f}, 0, 0]^\mathrm{T} \end{aligned}$$

(29)

where $\frac{\partial f}{\partial x}$ is the Jacobian of state Eqs.(24)–(27) and $\frac{\partial \mathcal {L}}{\partial x} = [\theta _{B}, \theta _{B}, 0, 0]$, then

$$\begin{aligned} u^{*}_{1}(t) = \left\{ \begin{matrix} 0, &{}\quad \text {if}\quad \lambda _{3}^{*}(t) \ge 0\\ -\frac{q \lambda _{3}^{*}(t)}{\theta _{u}},&{} \quad \ \text {if}\quad -\frac{\theta _{u}}{q}\le \lambda _{3}^{*}(t)< 0\\ 1,\ &{} \quad \text {if}\quad \lambda _{3}^{*}(t) < -\frac{\theta _{u}}{q}. \end{matrix}\right. \end{aligned}$$

(30)

Proof

According to PMP, if $u_{1}^{*}$ is an optimal control that solves problem (28), and if $x^{*}(t)$ and $\lambda ^{*}(t)$ are the corresponding state trajectory and costate trajectory, then the following equations are satisfied:

$$\begin{aligned}&\text {State equation}:\ \dot{x}^{*} = f(x^{*}, u_{1}^{*}) \end{aligned}$$

(31)

$$\begin{aligned}&\text {Costate equation}:\ \dot{\lambda }^{*} = -\left( \frac{\partial f}{\partial x}(x^{*})\right) ^\mathrm{T}\lambda ^{*} - \left( \frac{\partial \mathcal {L}}{\partial x}(x^{*}, u_{1}^{*})\right) ^\mathrm{T} \end{aligned}$$

(32)

$$\begin{aligned}&\text {Maximum principle}: \forall t\in [0,t_{f}],\ \mathcal {H}(x^{*}(t), \lambda ^{*}(t), u_{1}^{*}(t))\nonumber \\&\quad = \text {min}\ \{\mathcal {H}(x^{*}(t),\lambda ^{*}(t), u_{1})~|~u_{1}\in U_{1} \} \end{aligned}$$

(33)

$$\begin{aligned}&\text {Terminal condition}:\ \lambda ^{*}(t_{f}) =\left( \frac{\partial g(x^{*}(t_{f}))}{\partial x}\right) ^\mathrm{T}, \end{aligned}$$

(34)

where $\mathcal {H}$ is the Hamiltonian with the form of $\mathcal {H}(x, \lambda , u_{1}) = \lambda ^\mathrm{T} f(x,u_{1}) + \mathcal {L}(x,u_{1})$. We find that $\mathcal {H}(x, \lambda , u_{1}) = \widetilde{Q}\ +\ q\lambda _{3}u_{1} + (\theta _{u}/2) u_{1}^{2}$, where $\widetilde{Q}$ is the collection of terms that has no argument in $u_{1}$. Minimizing $\mathcal {H}(x, \lambda , u_{1})$ over $u_{1} \in U_{1}$ yields Eq. (30). The costate equation with terminal condition is

$$\begin{aligned} \dot{\lambda }^{*} = -\left( \frac{\partial f}{\partial x}(x^{*})\right) ^\mathrm{T}\lambda ^{*} - \left( \frac{\partial \mathcal {L}}{\partial x}(x^{*}, u_{1}^{*})\right) ^\mathrm{T},\ \lambda ^{*}(t_{f}) =\left( \frac{\partial g(x^{*}(t_{f}))}{\partial x}\right) ^\mathrm{T}, \end{aligned}$$

where $\partial \mathcal {L}/\partial x = [\theta _{B}, \theta _{B}, 0, 0]$, $\lambda ^{*}(t_{f}) = [\theta _{f}, \theta _{f}, 0, 0]^\mathrm{T}$. The Jacobian is

$$\begin{aligned} \frac{\partial f}{\partial x} = \begin{bmatrix} J_{11}&{}\quad J_{12} &{}\quad J_{13} &{}\quad J_{14}\\ J_{21}&{}\quad J_{22} &{}\quad 0 &{}\quad J_{24}\\ J_{31}&{}\quad 0 &{}\quad J_{33} &{}\quad 0\\ J_{41}\quad &{} J_{42} &{}\quad 0 &{}\quad J_{44} \end{bmatrix}, \end{aligned}$$

where

$$\begin{aligned} J_{11}&= \frac{r(1 - \mu )(k_{CD} - 2x_{1}-x_{2})}{k_{CD}} - k_{PD}\mathcal {I}(x_{3}) - \frac{\widetilde{\epsilon }(1 + x_{2})x_{4}}{(1 + x_{1} + x_{2})^{2}}\ ,\\ J_{12}&= \frac{\widetilde{\epsilon }x_{1}x_{4}}{(1 + x_{1} + x_{2})^{2}} - \frac{r(1 - \mu )x_{1}}{k_{CD}} \\ J_{13}&= - k_{PD}\frac{\psi x_{1}}{(1 + x_{3})^{2}}\ ,\ J_{14} = - \frac{\widetilde{\epsilon }x_{1}}{1 + x_{1} +x_{2}} \\ J_{21}&= \frac{\mu r (k_{CD} - 2x_{1} - x_{2})}{k_{CD}} - \frac{r'x_{2}}{k_{CD}} + \frac{\widetilde{\epsilon }x_{2}x_{4}}{(1 + x_{1} + x_{2})^{2}}\ ,\\ J_{22}&= \frac{r'(k_{CD} - x_{1} - 2x_{2})}{k_{CD}} - \frac{\mu r x_{1}}{k_{CD}} - \frac{\widetilde{\epsilon }(1 + x_{1})x_{4}}{(1 + x_{1} + x_{2})^{2}}\\ J_{24}&= - \frac{\widetilde{\epsilon }x_{2}}{1 + x_{1} +x_{2}}\ ,\ J_{25} = - k_{PD}\frac{\psi x_{2}}{(1 + x_{5})^{2}} \\ J_{31}&= \beta \mathcal {I}(x_{3}) - \psi x_{3}\ ,\ J_{33} = \frac{\beta \psi x_{1}}{(1 + x_{3})^{2}} - \omega - \psi x_{1}\\ J_{41}&= \frac{\alpha k_{ND}x_{4}(1 - x_{4})}{(k_{ND} + x_{1} + x_{2})^{2}}\ ,\ J_{42} = \frac{\alpha k_{ND}x_{4}(1 - x_{4})}{(k_{ND} + x_{1} + x_{2})^{2}}\ ,\\ J_{44}&= \alpha (\frac{x_{1} + x_{2}}{x_{1} + x_{2} + k_{ND}})(1 - 2x_{4}). \end{aligned}$$

We have explained Eqs. (31)–(34). $\square$

Fig. 7

(Colour figure online) Comparison of time series of population densities with different treatments in the high level of baseline immune response, $I_{0} = 8.5 \times 10^{6}$ cell/g. a Optimal injection rate, $\rho _{S}(t)$, is obtained by solving control problem (28) with tuned regulator weight $\theta _{u} = 10^{11}$ (the largest regulator weight on treatment costs). The Hamiltonian-based algorithm is terminated after k iterations and output control $u_{1}^{k}$, where $k = 11$. The numerical value of objective cost in problem (28) with control $u_{1}^{k}$ is 0.1024, and the convergence indicator $|\Theta (u_{1}^{k})| \approx 1.68 \times 10^{-8}$. Bacteria are eliminated around 30 h post-infection. b Optimal injection rate, $\rho _{S}(t)$ and $\rho _{R}(t)$, is obtained by solving control problem (17) with tuned regulator weight $\theta _{u} = 10^{11}$. The Hamiltonian-based algorithm is terminated after k iterations and output control $u^{k}$, where $k = 10$. The numerical value of objective cost in problem (17) with control $u^{k}$ is 0.1024, and the convergence indicator $|\Theta (u^{k})| \approx 9.34 \times 10^{-7}$. Note that the optimal injection rate of phage $P_{R}$ is nearly zero, i.e., $\rho _{R}(t) \approx 0\ \forall t \in [t_{0}, t_{f}]$. Thus, the optimal injection rates solved from 2D-OC and 1D-OC are nearly identical, i.e., $\rho _{S}(t)$ is a single-pulse signal centered at $t = 2$ h. Bacteria are eliminated around 30 h post-infection. c The practical therapeutic treatment is obtained from optimal injection rate in (b): Single dose, $P_{S}$ phage dose, is injected at 2 h post-infection with amount of $5 \times 10^{2}$PFU. Bacteria are eliminated around 30 h post-infection. See model parameters and simulation details in Sect. 4.2

Appendix C: Effective Single-Dose Treatment in Immunodeficient Hosts (Baseline Immune Response is Sufficiently High)

When the baseline immune response is sufficiently high in immunodeficient hosts, all the treatment strategies (1D-OC, 2D-OC, and practical treatments) can eliminate bacteria with a low dose of phage $P_{S}$ injected at very beginning of treatment (see Fig. 7 in Appendix 1).

Appendix D: Robustness analysis of optimal timing and dose to variations in therapy duration

From our numerical simulations, e.g., Figs. 3, 5, 6, and 7, we observe that successful phage therapy relies on early injection (for both single-dose and multi-dose cases). Here, we show that early ‘hit-hard’ approaches remain robust to variations in treatment duration when the final treatment time $t_{f}$ is changed. Please refer to Fig. 8 for the optimal dose and timing of the practical therapy for different final treatment time.

Fig. 8

(Colour figure online) The optimal timing and dose in practical therapeutic treatment with variation of final time from 2 to 4 days (i.e., $t_{f} \in [48, 96]$ h). The baseline immune response is fixed at $I_{0} = 6 \times 10^{6}$ cell/g. (Left) Minimal phage amount for eliminating bacterial cells with variation in final time $t_{f}$. Optimal dosages of phage $P_{S}$ (red) and phage $P_{R}$ (blue) are maintained at approximately $10^{8}$ (PFU/g) and $10^{7}$ (PFU/g), respectively. (Right) Optimal timing (defined by the peak of the optimal phage injection profile) of two types of phage injection with variation in final time $t_{f}$. The timings of injecting two types of phage dose are both about 2 h post-infection (i.e., $T_{P_{S}} = T_{P_{R}} \approx 2$ h). See model parameters and simulation details in Sect. 4.2