# Multi-strain phage induced clearance of bacterial infections

Jacopo Marchi, Chau Nguyen Ngoc Minh, Laurent Debarbieux, Joshua S Weitz

► To cite this version:

Jacopo Marchi, Chau Nguyen Ngoc Minh, Laurent Debarbieux, Joshua S Weitz. Multi-strain phage induced clearance of bacterial infections. PLoS Computational Biology, 2025, 21 (2), pp.e1012793. 10.1371/journal.pcbi.1012793. hal-05215158

HAL Id: hal-05215158

https://hal.science/hal-05215158v1

Submitted on 19 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Multi-strain phage induced clearance of bacterial infections

Jacopo MarchiID¹*, Chau Nguyen Ngoc Minh²,³, Laurent Debarbieux ID²,

Joshua S. Weitz ID¹,⁴,⁵*

¹ Department of Biology, University of Maryland, College Park, Maryland, United States of America,

² Institut Pasteur, Université Paris Cité, CNRS UMR6047, Bacteriophage Bacterium Host, Paris, France,

³ Sorbonne Université, Collège Doctoral, Paris, France, ⁴ Department of Physics, University of Maryland, College Park, Maryland, United States of America, ⁵ University of Maryland Institute for Health Computing, North Bethesda, Maryland, United States of America

* Former affiliation: School of Biological Sciences, Georgia Institute of Technology, Atlanta, Georgia, United States of America

* jsweitz@umd.edu

<figure><img src="image_2.png" /></figure>

## Abstract

Bacteriophage (or 'phage' – viruses that infect and kill bacteria) are increasingly considered as a therapeutic alternative to treat antibiotic-resistant bacterial infections. However, bacteria can evolve resistance to phage, presenting a significant challenge to the near- and long-term success of phage therapeutics. Application of mixtures of multiple phages (i.e., 'cocktails') has been proposed to limit the emergence of phage-resistant bacterial mutants that could lead to therapeutic failure. Here, we combine theory and computational models of in vivo phage therapy to study the efficacy of a phage cocktail, composed of two complementary phages motivated by the example of *Pseudomonas aeruginosa* facing two phages that exploit different surface receptors, LUZ19v and PAK_P1. As confirmed in a Luria-Delbrück fluctuation test, this motivating example serves as a model for instances where bacteria are extremely unlikely to develop simultaneous resistance mutations against both phages. We then quantify therapeutic outcomes given single- or double-phage treatment models, as a function of phage traits and host immune strength. Building upon prior work showing monophage therapy efficacy in immunocompetent hosts, here we show that phage cocktails comprised of phage targeting independent bacterial receptors can improve treatment outcome in immunocompromised hosts and reduce the chance that pathogens simultaneously evolve resistance against phage combinations. The finding of phage cocktail efficacy is qualitatively robust to differences in virus-bacteria interactions and host immune dynamics. Altogether, the combined use of theory and computational analysis highlights the influence of viral life history traits and receptor complementarity when designing and deploying phage cocktails in immunocompetent and immunocompromised hosts.

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Marchi J, Minh CNN, Debarbieux L, Weitz JS (2025) Multi-strain phage induced clearance of bacterial infections. PLoS Comput Biol 21(2): e1012793. https://doi.org/10.1371/journal.pcbi.1012793

**Editor:** Rob J De Boer, Utrecht University, NETHERLANDS, KINGDOM OF THE

**Received:** September 30, 2024

**Accepted:** January 14, 2025

**Published:** February 4, 2025

**Copyright:** © 2025 Marchi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data availability statement:** The software and data are available at https://github.com/Jacopo-Marchi/multi-strain-infection-clearance

**Funding:** The work was supported by the National Institutes of Health (1R01AI46592-01 to JSW, LD) and by the Simons Foundation (930382 to JSW). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
**Competing interests:** The authors have declared that no competing interests exist.

## Author summary

Phage represent a therapeutic alternative to treat antibiotic-resistant pathogenic infections. However, bacteria can evolve resistance to phage, presenting a significant challenge to phage therapy. Application of mixtures of multiple phages (i.e., phage cocktails) has been proposed to prevent treatment failure. Here, we combine theory, experiments, and computational models as a means to develop a pathway towards actionable principles to assess the efficacy of phage cocktails within immunomodulated hosts. We develop and analyze a flexible modeling framework to represent tripartite interactions between multiple phages, susceptible and resistant bacteria, along with the response of the innate immune system. We show that cocktails comprised of phage targeting independent bacterial receptors extend the range of immunocompromised conditions in which phage therapy can clear resistant pathogens, bolstering the public health relevance of phage cocktail development. We quantify the robustness of therapy success given variation in immune state, phage adsorption rates, and asymmetry in phage effectiveness within cocktails. We show that moderate adsorption rate variations can have drastic effects on therapy outcomes, potentially making the difference between a successful and a failing treatment. Such sensitivity highlights the importance of selecting phages that are highly effective against the target pathogens and potential resistant mutants.

## Introduction

Bacteriophage (phage, i.e. viruses that infect bacteria) are the most abundant organisms on the planet[1,2]. Bacteria have evolved a myriad of defense mechanisms to tolerate, counter, and resist infections [3,4]. Likewise, phages have co-evolved with bacteria for billions of years, catalyzing the emergence of a myriad of phage variants [5-8]. Diverse phage represent a largely untapped, therapeutic reservoir for treatment of antibiotic-resistant bacteria [9-12]. The therapeutic application of lytic phage is gaining interest as a viable treatment in alternative to antibiotics [13], capitalizing on the evolutionary dynamics of phages to combat bacterial pathogens that have outpaced traditional treatment options. Unlike antibiotics, which exert selective pressure driving the emergence of resistance, phage possess the inherent ability to evolve alongside bacteria, and can potentially adapt to overcome resistance mechanisms and maintain their efficacy over time [9]. Although evolution of phage resistance in bacteria is feasible, the widespread application of phage therapy is unlikely to catalyze broad spectrum resistance, given the high host specificity of phage [14].

In practice, phage therapy has shown significant potential in *in vivo* studies [15,16] and in human treatment in compassionate cases [17,18] and clinical trials [19]. The study [18] demonstrated that a personalized phage cocktail successfully treated a life-threatening, multidrug-resistant *Acinetobacter baumannii* infection in a critically ill patient, where conventional antibiotics had failed. However, there have been cases where therapeutic phage were not shown to be significantly more effective in controlling bacterial infections than standard procedures of care [20,21]. In [20] the authors found that a phage cocktail was safe and well-tolerated in treating *Pseudomonas aeruginosa* burn wound infections, but it did not demonstrate a significant improvement in bacterial load reduction or wound healing compared to standard treatment. The mixed results of clinical trials call for a thorough identification of the factors leading to success or failure of phage therapy [22].
To address this gap, several studies have integrated *in vitro* and *in vivo* experiments, as well as *in silico* computational models, into the design and assessment of phage therapy [23]. These studies have advanced the quantitative understanding of parameters impacting phage therapy outcomes, such as the pathogen strain, infected tissue, and chosen delivery method for the therapeutic phage [24]. A key factor modulating the outcome of *in vivo* phage therapy is the mammalian host immune response. For example, combined use of theory and *in vivo* application of phage in immunomodulated murine hosts revealed that phage and neutrophils work synergistically to eliminate *P. aeruginosa* and prevent the onset of fatal acute pneumonia [16]. In immunophage synergy, phage rapidly infect and lyse susceptible *P. aeruginosa* cells while neutrophils clear both susceptible and subpopulations of phage-resistance cells. However, not all phage and immune interactions may lead to positive impacts on therapy. For instance, alveolar macrophages can reduce the density of circulating phage, potentially jeopardizing phage therapeutic treatment efficacy [25]. When immune systems are compromised [26] or limit phage-induced clearance of target bacteria [25], the proliferation of phage-resistance mutants can lead to therapeutic failure [16].

Indeed phage resistance represents one of the major challenges to phage therapy success. The emergence of bacteria resistant to virulent phage reduces the ability of phage to contain an infection [27–30]. There are multiple approaches to overcome the evolution and proliferation of phage resistant bacteria. First, phage can be trained via asymmetric evolutionary training to target evolved bacteria that are resistance to the original, therapeutic phage [31]. Subsequent application of trained phage can limit the emergence of phage-resistant bacteria, improving therapeutic efficacy. Alternatively, phage may be used in combination with small molecules (e.g., antibiotics) that limit the potential for bacterial evolution. In one well-studied case, phage OMKO1 targets antibiotic efflux pumps within *P. aeruginosa*; hence joint use of phage OMKO1 and antibiotics can lead to therapeutic success [32]. Finally, there may be cases where phage target receptor sites such that evolution of phage resistance comes with significant fitness costs in an *in vivo* context [9,33–38]. All of these examples serve to illustrate the general rule that decreasing the scope of phage-escape bacterial mutants and/or increasing the costs of phage-resistant mutations increase the efficacy of phage therapy. In the same spirit, the combined use of multiple phages in a cocktail that infect distinct receptors makes it harder for bacteria to acquire resistance against all phages at the same time [39–42]. Cocktails with phages targeting different cell receptors have higher efficacy against *P. aeruginosa* [43]. These evolutionary principles can be combined in cocktails that comprise trained phage, which can efficiently control *P. aeruginosa* populations both *in vitro* and in mice lung infections [44].

In this work we focus on evaluating the therapeutic efficacy of phage cocktails, which leverage complementary adsorption paths to infect target bacteria. We combine computational and analytical treatment of single- and double-phage therapy models, informed by Luria-Delbrück (LD) fluctuation tests exposing *P. aeruginosa* to phages LUZ19v and PAK_P1, to address whether phage cocktails can help restore therapy success in immunocompromised hosts in face of phage resistance. We then extend our analysis to include *in vivo* immune system dynamics and structured phage-bacteria interactions, mapping the quantitative conditions for treatment success as a function of the immune state and therapeutic phages effectiveness. Comparing our exploration with the parameters extracted in [16] we note that the quantitative adsorption rate of any of the therapeutic phages may have drastic effects on the treatment outcomes jeopardizing the benefits of phage cocktails, highlighting the importance of selecting a combination of efficient phages against the pathogenic strains especially in hosts with a compromised immune system.
# Methods

## Fluctuation test to infer the probability of resistance mutations of *P. aeruginosa* against two phages

We perform a LD fluctuation test to assess the likelihood that *P. aeruginosa* colonies randomly develop a resistance mutation against either phage PAK_P1 or LUZ19v or against both phages at the same time. For each phage, we grow overnight 360 independent populations of strain PAK-lumi. 300 of these are plated against either PAK_P1 or LUZ19v to count the number of resistant mutants, whereas the remaining 60 are used as control to count the CFUs in the absence of phage. To measure resistance against both phages simultaneously we grow 150 independent colonies and plate them against a mixture of PAK_P1 and LUZ19v. More details on the fluctuation test protocol are in Sec. 2 in S1 Text. Then we use the resistant mutants counts to infer the probability at which cells develop resistance to PAK_P1, LUZ19v or both, within one duplication. We run the inference through the web-tool bz-rates [45], which learns the model parameters via the Generating Function estimator from [46] (see Sec. 3 in S1 Text for more details on the inference).

## Mathematical models of phage therapy

**Summary.** We study two different models of single- and double-phage therapy of bacterial infections in immunomodulated hosts leading to four different scenarios. Fig 1 sketches the models ingredients both in terms of ecological dynamics (panel A) as well as evolutionary multi-strain phage-bacteria interactions (panel B). We start from a simple model of a lung infection, mimicking a signaling deficient immune system as in the case of myeloid differentiation primary response gene 88-deficient mice (MyD88-/-), where the immune system is static as neutrophils are not recruited in the lungs [16] (Fig 1, A₁). We first analyze a single-phage treatment model (Fig 1, B₁), then we add a second phage to the cocktail to address its impact on treatment outcome (Fig 1, B₂ top). Third, switching back to a single-phage therapy, we include more complex ecological dynamics to account for a responsive immune system, still modulated by its carrying capacity (limiting the availability of immune cells within the infected tissue). In this second single-phage treatment model we include intermediate stages of infected bacteria and add a non-linearity in phage infections to mimic a Michaelis-Menten phage-bacteria binding reaction (Fig 1, A₂). Finally we study the effect of a phage cocktail with the inclusion of these new model ingredients, and explore the role of evolutionary relations between the combined therapeutic phages (Fig 1, B₂ bottom). We map the phage therapy outcomes as a function of host immune strength and phage efficacy, depending on the different components of these four models.

**Infection model of interacting phage, multi-strain bacteria and non-responsive host immune system.** The following ordinary differential equation model represents how phage *P* interacts with susceptible bacteria *S* and phage-resistant bacteria *R*, under the action of immune cells *I* during infection of an immunodeficient host:

$$ 
\begin{aligned}
\frac{dS(t)}{dt} &= rS - \tilde{d}S(S + R) + \mu r(R - S) - \frac{\kappa IS}{1 + \frac{S+R}{K_D}} - \phi SP, \\
\frac{dR(t)}{dt} &= rR - \tilde{d}R(S + R) + \mu r(S - R) - \frac{\kappa IR}{1 + \frac{S+R}{K_D}}, \\
\frac{dP(t)}{dt} &= \beta \phi SP - \omega P,
\end{aligned}
\quad (1) $$
<figure><img src="image_4.png" /><figcaption>**Fig 1. Schematic of mathematical model of phage therapy.** A) Ecological interactions between bacteria, phage and immune system (graphical and mathematical symbols in the legend box). Subpanel A₁ shows the interactions in the simple phage therapy model in a signaling deficient host introduced in Section “Infection model of interacting phage, multi-strain bacteria and non-responsive host immune system”, whereas Subpanel A₂ sketches the addition of more complex dynamics in the model introduced in Section “Structured model for single-phage treatment of *in vivo* infections with a modulated responsive immune system”, namely a structured bacteria population including stages of phage infection and the active recruitment of the host immune cells. Panel B) gathers the cross-infection and mutation networks between specific phage and bacteria strains (listed in the legend box) in the different treatment models. Subpanel B₁ sketches the single-phage treatment structure used in Sections “Infection model of interacting phage, multi-strain bacteria and non-responsive host immune system” and “Structured model for single-phage treatment of *in vivo* infections with a modulated responsive immune system”, whereas B₂ summarizes the cross-infection structure in the phage-cocktail treatment models introduced in Sections “Model of a phage combination therapy against multiple pathogenic bacteria strains” and “Multi-phage treatment model, modulating the evolutionary interactions between bacteria and phage strains during infections” (top and bottom respectively). The main model parameters studied in this work, the immune strength *I* and the phage adsorption rate *φ*, are colored in red. Arrow types encode the different kind of dynamical interactions (growth, killing, decay, mutations), with dashed-dotted lines indicating saturating nonlinear rates of the form $f(\cdot) \propto \frac{\phi}{1+p}$. Created in BioRender. NGUYEN NGOC MINH, C. (2025) https://BioRender.com/y51s627</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012793.g001

This model includes the main ingredients pinpointed in early theoretical studies of phage-
bacteria interactions *in vitro* [47] and in the context of phage therapy [48]. Bacteria undergo
logistic growth at maximum rate *r* and with density driven death rate *d*, equivalent to a
carrying capacity term $K_C = r/\bar{d}$. They mutate between the susceptible and resistant types, and the susceptible type can be infected by phage with adsorption rate $\phi$. At the same time phage lyse susceptible bacteria producing $\beta$ virions per lysis event, and decay at a constant rate $\omega$. We build on previous models focusing on acute infections [16,26], therefore we only include the action of innate immune effector cells $I$, mainly representing the action of neutrophils. In this first model we take $I$ as a constant mimicking an immunodeficient host unable to recruit more immune cells other than a fixed baseline as in the case of myeloid differentiation primary response gene 88-deficient mice (MyD88$^{-/}$), as proposed in [16]. Immune cells target both bacteria types through a saturating function, as bacteria evade the immune action at high population densities compared to an evasion threshold $K_D$, consistent with experiments showing that pathogens at high density can activate defenses against the immune system [49,50] through physical shielding [51] or expression of virulence factors [52]. This modeling ingredient, included in [53] to explain the saturation in granulocytes action on *P. aeruginosa* during mice thigh infections and then adapted to the context of phage therapy [26], was shown to be essential to explain the synergy between phage and immune system in clearing bacteria infections [16,26]. Therefore, model (1) represents a baseline model of immunophage synergy extending [26] to enable assessment of phage-resistant bacteria on *in vivo* dynamics given modulation of immune strength ($I$) and phage life history traits (here, primarily explored through variation in $\phi$). Table 1 reports the parameters meaning and values.

**Model of a phage combination therapy against multiple pathogenic bacteria strains.**

In the second model we add a second phage to the infection treatment, now composed of a combination of two phages $P_1$ and $P_2$ (Fig 1, B$_2$ top). The pathogen types include the

<table><caption>**Table 1. Model parameters.** Most values are taken from [26]. The mutation rate is inferred from the fluctuation test (see Sec. 3 in S1 Text). For the viral decay $\omega$ we elect a mid-value between the parameter reported in [26,54] (1 h$^{-1}$) and the one inferred in [16] (0.07 h$^{-1}$).</caption><thead><tr><th>Variable</th><th>Meaning</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>r</td><td>Maximum bacterial growth rate</td><td>1</td><td>h<sup>-1</sup></td></tr><tr><td>d̄</td><td>Bacteria density-driven death rate</td><td>10<sup>-10</sup></td><td>g/(h CFU)</td></tr><tr><td>μ</td><td>Mutation probability per replication</td><td>4.3 · 10<sup>-8</sup></td><td></td></tr><tr><td>φ</td><td>Phage adsorption rate</td><td>[10<sup>-12</sup>, 10<sup>-6</sup>]</td><td>g/(h PFU)</td></tr><tr><td>β</td><td>Burst size</td><td>100</td><td>(PFU/cell)</td></tr><tr><td>ω</td><td>Viral decay</td><td>0.4</td><td>h<sup>-1</sup></td></tr><tr><td>κ</td><td>Immune killing rate</td><td>8.2 · 10<sup>-8</sup></td><td>g/(h CFU)</td></tr><tr><td>I</td><td>Immune cells density</td><td>[10<sup>4</sup>, 10<sup>9</sup>]</td><td>(immune cell)/g</td></tr><tr><td>K<sub>D</sub></td><td>Cell shielding density</td><td>10<sup>7</sup></td><td>(CFU)/g</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1012793.t001
wild-type infecting bacteria S, susceptible to both phages, which can in turn mutate to a type $R_1$ resistant to $P_1$ or to a type $R_2$ resistant to $P_2$:

$$ 
\begin{aligned}
\frac{dS(t)}{dt} &= rS - \tilde{d}SB_{tot} + \mu r(R_1 + R_2 - S) - \frac{\kappa IS}{1 + \frac{B_{tot}}{K_D}} - \phi SP_1 - \phi SP_2, \\
\frac{dR_1(t)}{dt} &= rR_1 - \tilde{d}R_1 B_{tot} + \mu r(S - R_1) - \frac{\kappa IR_1}{1 + \frac{B_{tot}}{K_D}} - \phi R_1 P_2, \\
\frac{dR_2(t)}{dt} &= rR_2 - \tilde{d}R_2 B_{tot} + \mu r(S - R_2) - \frac{\kappa IR_2}{1 + \frac{B_{tot}}{K_D}} - \phi R_2 P_1, \\
\frac{dP_1(t)}{dt} &= \beta \phi (S + R_2) P_1 - \omega P_1, \\
\frac{dP_2(t)}{dt} &= \beta \phi (S + R_1) P_2 - \omega P_2.
\end{aligned}
\quad (2) $$

We initially assume that the two phages have the same adsorption rate $\phi$, which we vary together with the immune strength $I$ to evaluate therapeutic outcomes when using a phage cocktail. Importantly, model (2) assumes that bacteria cannot develop resistance to both phages at the same time within an infection timescale, which was suggested in several studies proposing cocktail of phages that infect bacteria through different routes [42,43,55]. We also assume that there is no cross-infectivity between the two phages with respect to the resistant types. The susceptible bacteria mutate equally to either resistant type (and vice versa), but resistant types do not mutate directly among themselves – hence we ignore the near-term chance of the emergence of a double-phage resistant mutant.

**Structured model for single-phage treatment of in vivo infections with a modulated responsive immune system.** We extend the single phage therapy model Eq. (1) to include a realistic in vivo innate immune response and more complex phage infection dynamics [16] (Fig 1, A2):

$$ 
\begin{aligned}
\frac{dS(t)}{dt} &= rS - \tilde{d}SB_{tot} + \mu r(R - S) - \frac{\kappa i S}{1 + \frac{B_{tot}}{K_D}} - S\phi F(P_1), \\
\frac{dR(t)}{dt} &= rR - \tilde{d}RB_{tot} + \mu r(S - R) - \frac{\kappa i R}{1 + \frac{B_{tot}}{K_D}}, \\
\frac{dE_1^1(t)}{dt} &= \phi SF(P_1) - L\eta E_1^1, \\
\frac{dE_1^i \in \{2...L\}(t)}{dt} &= L\eta E_1^{i-1} - L\eta E_1^i, \\
\frac{dP_1(t)}{dt} &= \beta L\eta E_1^L - \phi SF(P_1) - \omega P_1, \\
\frac{di(t)}{dt} &= \alpha i \left(1 - \frac{i}{I}\right) \frac{B_{tot}}{B_{tot} + K_N}.
\end{aligned}
\quad (3) $$

The innate immune cells $i$ are recruited in the lungs at a maximum rate per capita $\alpha$ until they reach a saturation density $I$ (immune capacity), which is one of the parameters we vary in this study to modulate hosts immune responses. Immune cells are recruited proportionally to a saturating function of the total density of bacteria $B_{tot}$ with half saturation constant $K_N$. These saturation profiles in immune recruitment within infected tissues are supported by empirical evidence that neutrophils can only be produced up to a limit [56], and that innate immune responses saturate as a function of bacterial loads [53]. Further evidence supports that killing
of bacteria saturates at high neutrophils counts [57]. The immune dynamics in Eq. 3 build on previous phage therapy models [16,26], which have shown that limitation on immune clearance at high bacterial densities and self-limitation of the immune response are necessary ingredients to recapitulate the empirical evidence of bacterial persistence within the host in specific treatment contexts [16].

Additionally, we structure the bacteria population introducing L infected bacteria stages, $E_1^{i \in \{1...L\}}$, between phage adsorption and lysis to model a finite phage infection time distributed as an Erlang distribution with average $\frac{1}{\eta}$, an approach applied in several models of epidemiological dynamics and phage-bacteria interactions [58–60]. We set $L = 10$ as in previous studies to ensure a coefficient of variation for the lysis time distribution of $1/\sqrt{10} \sim 0.3$, in line with empirical estimates [60–62]. Finally, as proposed in [16], we add a saturation term in phage adsorption $F(P_1) = \frac{P_1}{1+P_1/P_c}$, as if phage adsorbed into bacteria through Michaelis-Menten reaction kinetics [63,64], producing a phage-adsorption profile similar to Monod growth [65]. Table 2 provides a list of parameter definitions and their values.

**Multi-phage treatment model, modulating the evolutionary interactions between bacteria and phage strains during infections.** Finally, we develop a multi-phage therapy model composed of two phages $P_1$ and $P_2$, in the presence of the new modeling ingredients presented in the previous section. As in eq. (2), bacteria susceptible to both phages, S, can mutate to a type $R_1$ resistant to $P_1$ or to a type $R_2$ that can be infected by $P_1$. Compared to Section “Model of a phage combination therapy against multiple pathogenic bacteria strains”, now we introduce a parameter $p$ continuously modulating the evolutionary relations between the two phages and the target bacteria, in an asymmetric fashion. $P_2$ can either infect $R_1$ with probability $p$, or $R_2$ with probability $1-p$. So when $p=1$ the two phages are complementary in infecting each other’s resistant type with rate $\phi$, and $R_i$ truly has the meaning of bacteria resistant to phage $i$. When $p=0$ the two phages are phenotypically equivalent in the sense that they target the same set of bacteria, as if resistance would arise to both at the same time. Intermediate values of $p$ mimic a generalist phage $P_2$ capable of infecting either bacterium, but with

<table><caption>**Table 2. Model parameters.** Most values are taken from [16]. The mutation rate is inferred via a fluctuation test (see Sec. 3 in S1 Text). The lysis rate $2h^{-1}$ is compatible with reported lysis times of phages infecting *P. aeruginosa* [66,67]. The number of infected stages $L$ is the same as several previous models of phage-bacteria dynamics [60,68].</caption><thead><tr><th>Variable</th><th>Meaning</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>r</td><td>Maximum bacterial growth rate</td><td>0.75</td><td>h<sup>-1</sup></td></tr><tr><td>d̄</td><td>Bacteria density-driven death rate</td><td>7.5 · 10<sup>-9</sup></td><td>g/(h CFU)</td></tr><tr><td>μ</td><td>Mutation probability per replication</td><td>4.3 · 10<sup>-8</sup></td><td></td></tr><tr><td>φ</td><td>Phage adsorption rate</td><td>[10<sup>-9</sup>, 10<sup>-6</sup>]</td><td>g/(h PFU)</td></tr><tr><td>P<sub>c</sub></td><td>Phage adsorption saturation density</td><td>1.5 · 10<sup>7</sup></td><td>(PFU)/g</td></tr><tr><td>η</td><td>Lysis rate</td><td>2</td><td>h<sup>-1</sup></td></tr><tr><td>L</td><td>Number of infected stages</td><td>10</td><td></td></tr><tr><td>β</td><td>Burst size</td><td>100</td><td>(PFU/cell)</td></tr><tr><td>ω</td><td>Viral decay</td><td>0.07</td><td>h<sup>-1</sup></td></tr><tr><td>κ</td><td>Immune killing rate</td><td>8.2 · 10<sup>-8</sup></td><td>g/(h CFU)</td></tr><tr><td>I</td><td>Immune cells capacity</td><td>[10<sup>4</sup>, 10<sup>8</sup>]</td><td>(immune cell)/g</td></tr><tr><td>K<sub>D</sub></td><td>Bacteria shielding density from immune action</td><td>4.1 · 10<sup>7</sup></td><td>(CFU)/g</td></tr><tr><td>K<sub>N</sub></td><td>Bacteria density when immune recruitment is half its maximum</td><td>10<sup>7</sup></td><td>(CFU)/g</td></tr><tr><td>α</td><td>Maximum immune cells recruitment rate</td><td>0.97</td><td>h<sup>-1</sup></td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1012793.t002
lower specificity to either type since the sum of the effective adsorption rates to all bacteria is always $\phi$. The equations for this model are:

$$ 
\begin{aligned}
\frac{dS(t)}{dt} &= rS - \tilde{d}SB_{tot} + \mu r(R_1 + R_2 - S) - \frac{\kappa i S}{1 + \frac{B_{tot}}{K_D}} - S\phi F(P_1) - S\phi F(P_2), \\
\frac{dR_1(t)}{dt} &= rR_1 - \tilde{d}R_1B_{tot} + \mu r(S - R_1) - \frac{\kappa i R_1}{1 + \frac{B_{tot}}{K_D}} - R_1p\phi F(P_2), \\
\frac{dR_2(t)}{dt} &= rR_2 - \tilde{d}R_2B_{tot} + \mu r(S - R_2) - \frac{\kappa i R_2}{1 + \frac{B_{tot}}{K_D}} - R_2\phi F(P_1) - R_2(1-p)\phi F(P_2), \\
\frac{dE_1^1(t)}{dt} &= \phi SF(P_1) + \phi R_2F(P_1) - L\eta E_1^1, \\
\frac{dE_1^{i \in \{2...L\}}(t)}{dt} &= L\eta E_1^{i-1} - L\eta E_1^i, \\
\frac{dE_2^1(t)}{dt} &= \phi SF(P_2) + (1-p)\phi R_2F(P_2) + p\phi R_1F(P_2) - L\eta E_2^1, \\
\frac{dE_2^{i \in \{2...L\}}(t)}{dt} &= L\eta E_2^{i-1} - L\eta E_2^i, \\
\frac{dP_1(t)}{dt} &= \beta L\eta E_1^L - \phi SF(P_1) - \phi R_2F(P_1) - \omega P_1, \\
\frac{dP_2(t)}{dt} &= \beta L\eta E_2^L - \phi SF(P_2) - (1-p)\phi R_2F(P_2) - p\phi R_1F(P_2) - \omega P_2, \\
\frac{di(t)}{dt} &= \alpha i \left(1 - \frac{i}{I}\right) \frac{B_{tot}}{B_{tot} + K_N},
\end{aligned}
\quad (4) $$

where $E_1^i$ denotes a bacterium infected by $P_1$ (at infection stage $i$), and $E_2^i$ denotes a bacterium infected by $P_2$. We start by considering phages with the same life history traits, and then relax this assumption studying a case where $P_1$ adsorbs at a rate $\phi_1$ while $P_2$ adsorbs at a rate $\phi_2$. In this case we study a scenario where $F(P_1) = \frac{P_1}{1+\frac{P_1+P_2}{P_c}}$, mimicking apparent competition between the two phages, as if they bind to the same “substrate” or surface receptor, or if they shared common molecular machinery in order to inject genetic material into the cell.

## Results

### Condition for single-phage therapy success in the face of phage resistance and non-responsive immunity

We start by analyzing the single-phage model presented in Section “Infection model of interacting phage, multi-strain bacteria and non-responsive host immune system”. In the absence of phage, the model admits a stable fix point $B_1^S$ and an unstable one $B_1^U$ for bacteria densities. After the addition of phage, whenever bacteria are driven below the infectious dose $B_1^U$ then the immune system is able to clear the infection, producing a successful treatment with phages and immune system working in synergy to control the pathogens [16,26]. Neglecting the emergence of phage-resistant bacteria, the phage induced bacteria steady state becomes unstable when $B_p = \frac{\omega}{\beta\phi} < K_D (\sqrt{\frac{\kappa i K_c}{rK_D}} - 1)^2$ (derived in [26]). In terms of immune strength,
this means that for immunophage synergy to function then the host innate immunity needs to be bigger than a threshold $I_0$:

$$I > I_0 = \frac{\tilde{d}K_D}{\kappa} \left( 1 + \frac{\omega}{\beta\phi K_D} \right)^2, \qquad (5)$$

which in the absence of phage resistance would be a sufficient condition for therapeutic success.

In contrast, if phage-resistant bacteria are present in the population, then the phage-induced dynamical instability is not enough to drive bacteria down as the phage resistant bacteria grow to an immune-dependent stable density (see Sec. 5 in S1 Text). The mathematical condition $B_I^U > 0$ becomes necessary for therapy success together with Eq. (5), so that the bacteria extinction steady state $B^* = 0$ becomes stable. In this case, phage can drive the susceptible bacteria below $B_I^U$ before the resistant bacteria proliferate from rare to densities above this threshold. When bacteria densities are depleted below $B_I^U$, the immune system can clear the remainder of susceptible and resistant bacteria in the host. This mechanism restores the ability of phage and the host immune system to synergistically clear an infection in the presence of phage resistant bacteria, reinforcing phage-immune synergy against multiple strains of bacteria. Fig 2 shows how, with $B_I^U > 0$, crossing the threshold $I > I_0$ drives a transition from therapeutic failure to success as phage drive the susceptible bacteria below $B_I^U$ before the resistant type grows enough, at which point the immune system can clear both phage-susceptible and phage-resistant bacteria.

Next, we assume that the immune evasion threshold $K_D$ is lower than bacteria carrying capacity $K_C$ and that the immune system is well below the strength where it would clear the infection on its own (and hence no therapy would be necessary). In this case we can derive the

<figure><img src="image_5.png" /><figcaption>**Fig 2. Phage-immune synergy in face of phage resistance.** Simulated dynamics for phage, susceptible and resistant bacteria, when the immune system strength is below $I_0$ (left panel) or above (center, right panels), with phage (left, center) or without (right). Increasing the immune system strength above $I_0$, with $B_I^U > 0$, phage can drive susceptible bacteria below $B_I^U$ before the resistant type grows, driving a transition from therapy failure to success (left to center panel). This is a signature of phage-immune synergy as the infection would persist without phage (right). Simulation parameters are reported in Table 1 with $\phi = 10^{-11}$ g/(h PFU), $I = 1.5 \cdot 10^7$ and $2.9 \cdot 10^7$ cells/g respectively in the left and the other two panels.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012793.g002
immune strength attaining $B_I^U > 0$, necessary for immunophage synergy (see Sec. 5 in S1 Text for the derivation):

$$I > I_b = \frac{r}{\kappa}. \qquad (6)$$

If phage drive the susceptible bacteria below $B_I^U$ fast enough, the simultaneous satisfaction of the conditions in Eqs. (5) and (6), or more compactly

$$I > I_S = \max(I_0, I_b), \qquad (7)$$

yields phage therapy success. This condition summarizes the necessary host immune strength as a function of phage life history traits, so that the application of a single phage clears the infection in an immunomodulated host, overcoming the evolutionary challenge posed by bacteria developing resistance to phage. The intuitive explanation of this generalization of immunophage synergy is as follows. First, the immune system needs to be strong enough so phage therapy can control susceptible bacteria, represented by the condition in Eq. (5). Second, resistant bacteria must be controllable at relatively low densities by the immune system, Eq. (6). Combining the two we recover the complete success condition in Eq. (7).

We test our analytic results through several simulations of the model in Eq. (1) varying the immune system strength $I$ and the phage adsorption rate $\phi$. Fig 3 shows the phase diagram for the density of bacteria after 200 hrs post-infection, comparable to in vivo experiments, averaged over the last day (details in Sec. 4 in S1 Text). In this case bacteria are either eliminated in a successful treatment, or phage-resistant types grow to an immune system-dependent stable density and therapy fails. As predicted by Eq. (7), and confirmed by the numerical results, the immune system needs to be strong enough for a single-phage therapy to be successful and phage with a high adsorption rate can lead to therapy success in weaker hosts (lower $I$).

## Phage cocktails can improve therapy efficacy and overcome phage resistance in immunocompromised hosts

Immunophage synergy can eliminate a population of phage-resistant and phage-susceptible bacteria insofar as the immune system is sufficiently strong given the life history traits of the therapeutic phage (as summarized in Eqs.(5) and (6)). Likewise, monophage therapy can fail due to the proliferation of phage-resistant bacteria. Hence, we next explore how phage cocktails can potentially expand the range of immunocompromised hosts in which phage therapy is effective. To do so, we simulate model (2), which includes two phages that can target the susceptible bacteria as well as one additional bacterium (i.e., each of the resistant bacteria is only resistant to one of the two phages in the cocktails). This complementarity with respect to resistance reflects that, in our in vitro experiment growing 150 separate populations of *P. aeruginosa*, we found no mutant conferring simultaneous resistance against a cocktail of phages PAK_P1 and LUZ19v. Hence we assume that the emergence of simultaneous resistance against both phages from a background of susceptible bacteria is negligible.

Fig 4 A shows the expanded region of immunophage synergy given a phage cocktail via numerical simulations of model (2). The addition of a second phage can restore therapeutic success in immunocompromised hosts, provided phages have sufficiently effective life history traits (here, we focus on variation in the adsorption rate $\phi$). This is expected as the added phage will lyse the bacteria resistant to the other one, as shown in the population dynamics in Fig 4C,D, which show how bacteria strains evolve under the selection imposed by phages eventually leading to extinction. Even so, if phages have low adsorption rate then they will not clear the infection (e.g. see dynamics in Fig 4B).
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Bacteria density as a function of immune strength and phage adsorption rate. Numerical simulations of the model in Eq. (1) varying *I* and *φ*. The color map represents the density of bacteria in the last part of the numerical simulations. Single-phage therapy works such that bacteria are driven to elimination with a strong enough immune system (denoted in the white region). Efficient phage (higher *φ*) broaden the therapy success conditions. The analytic condition in Eq. (7) (blue dashed line) accurately predicts the therapy outcome transition. The two black triangles correspond to the left and central panel in Fig 2. Simulation parameters are reported in Table 1.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012793.g003

We propose a heuristic derivation of the phase space of immunophage synergy, as detailed in Sec. 6 in S1 Text. The result predicts that without the immune system ($I=0$) phages drive bacteria to the extinction threshold $T$ (1 bacterial cell in the lungs) for $\phi > \phi_c$, where $\phi_c$ satisfies the transcendental equation:

$$ \frac{2\phi_c P_0}{r} = 2 + \ln\left(\frac{\phi_c P_0}{r}\right) + \ln\left(\frac{\omega}{\phi_c \Omega \beta T}\right) \frac{\omega}{r} \quad (8) $$

Here $P_0$ is the initial phage density. $\Omega$ is a numerical factor ensuring bacteria exponential decay before extinction (see Sec. 6 in S1 Text), which needs to satisfy $\Omega \gg 1$ and $\Omega < \frac{\omega}{\phi_c \beta T}$. $\Omega$ affects the equation only logarithmically, and $\phi_c$ decreases with $\Omega$. The condition in Eq. (8) gives the lowest $\phi$ above which phage can drive bacteria below the extinction threshold $T$ before bacteria resume growth (due to phage decay at rate $\omega$). The leading factors determining the phage cocktail success are the initial rate of phage killing $\phi P_0$ and bacteria growth $r$, which drive the transition to phage-bacteria coexistence. The dashed vertical green line in Fig 4A shows $\phi_c$ for $\Omega = 10$. We observe that the numerical simulations agree with the heuristic prediction Eq. (8) when $I \ll I_b$. Note that when phage are not very efficient ($\phi \ll \phi_c$) the presence of resistance is not the major treatment failure driver, therefore Eq. (5) is an upper bound
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Phage cocktail therapy succeeds for a wide parameters range. A) Numerical simulations of the model in Eq. (2) varying *I* and *φ*. The color map represents the density of bacteria in the last part of the numerical simulations. Adding a second phage clears the infection for a much wider parameter range compared to single-phage treatment, if phages have high enough *φ*, as highlighted by the comparison with the transition for single-phage therapy success in Eq. (7) (blue dashed line). The vertical dashed green line represents the value of *φ*<sub>c</sub> obtained solving Eq. (8) with Ω = 10 The black annotations correspond to the parameters yielding the population dynamics in panels B),C),D) respectively that show the evolution of phage and bacteria strains during the simulated infection. Simulation parameters are reported in Table 1.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012793.g004

approximating the success condition. Together with the baseline single-phage treatment
model results, our analysis of a simple phage cocktail model reveals that phage can restore
therapeutic success even in immunodeficient hosts and even when phage-resistant bacte-
rial mutants are present, insofar as sufficiently effective phage – according to the quantitative
condition $\phi > \phi_c$ – are utilized.

Synergy between single-phage treatment and modulated responsive immune system against *in vivo* infections by multiple bacteria

The previous sections addressed the efficacy of monophage and phage cocktail therapy assum-
ing that interactions between phage and bacteria occur in a well-mixed environment. Here
we expand these results to a more realistic *in vivo* context by incorporating a nonlinear phage
adsorption profile (reported to fit experimental adsorption curves [63,64]), incorporat-
ing a structured phage infection dynamics, and including a dynamic response of the innate
immune system. With the inclusion of these ingredients the model structure becomes similar to previous *in vivo* models that quantitatively reproduced experimental dynamics measured when treating lung infections with phage in immunomodulated mice [16]. Therefore we can directly compare our theoretical results with the parameters previously inferred, obtaining quantitative insights on the dynamics produced by single-phage treatments of *in vivo* infections.

Fig 5 shows the final density of bacteria produced by numerical simulations of Eq. 3 after 200 hrs, for a broad range of immune system capacities I and phage adsorption rates $\phi$. We find a qualitatively similar pattern as in Section “Condition for single-phage therapy success in the face of phage resistance and non-responsive immunity”, as the immune system needs to be strong enough to clear the infection in synergy with phage, which in turn can improve treatment outcomes with higher adsorption rates. We can compare infection outcomes for the parameter values explored in Fig 5 to the parameters inferred in [16] for wild-type mice with a fully competent immune system (upward light-blue triangle) and for MyD88⁻/⁻ mice with a deficient immune condition preventing neutrophils recruitment above a low baseline (downward dark-green triangle). Fig 5 extrapolates the treatment outcome for biologically relevant parameters around these experimentally inferred values. Notably, we confirm that the current model predicts phage therapeutic success in immunocompetent mice, however the model also predicts that this outcome is contingent on the use of a phage with a sufficiently rapid adsorption rate.

To understand these numerical results quantitatively, we proceed as in Section “Condition for single-phage therapy success in the face of phage resistance and non-responsive immunity” by considering the system without phage $P = 0$, which yields the same fixed points for bacteria density while the immune cells relax to $i^* = I$. Hence Eq. (6) is again a necessary condition for therapy success, so that phage resistant bacteria can be controlled while still below the infectious dose $B_I^U$. To progress further in the analysis we assume that phage infection dynamics are fast compared to other timescales, disregarding the dynamics of infected bacteria $E_1^{i \in \{1...L\}}$, which corresponds exactly to the model in [16]. Assuming also a scenario in which phage are so abundant that $P \gg P_c$, and that the immune cells dynamics relax fast to I, we find (details of the derivation in Sec. 7 in S1 Text) a new immune threshold $I_S$ that constrains the host immune capacity for which the single-phage therapy succeeds:

$$I > I_S = \max(I_1, I_b), \qquad (9)$$

with

$$I_1 = \frac{r - \phi P_c}{\kappa} \left( 1 + \frac{B_0}{K_D} \right). \qquad (10)$$

As in Section “Condition for single-phage therapy success in the face of phage resistance and non-responsive immunity” the condition for therapeutic success arises from combining two necessary immune constraints, so that phage and immune system can synergistically control both the susceptible bacteria (Eq. (10)), and the phage resistant mutants (Eq. (6)). Notably, the condition to overcome phage resistance in a timely fashion is determined by the balance between bacteria growth $r$ and immune killing $\kappa$ regardless of the new model ingredients. Eq. (9) depends on the infecting bacteria density at the time of therapy administration $B_0$, highlighting the importance of prompt intervention when treating the infection. The dash blue line in Fig 5 shows Eq. 9 as a function of $\phi$. The analytical result agrees with the numerical results of the general model in Eq. 3 and with the findings of [16] where therapy worked
<figure><img src="image_8.png" /><figcaption>**Fig 5.** *Single-phage therapy succeeds against phage resistance with a strong enough immune system.* Numerical simulations of the model in Eq. (3) varying *I* and *φ*. The color map represents the density of bacteria in the last part of the numerical simulations, where the white region denotes bacteria clearance. Single-phage therapy works with a strong enough immune system, which needs to be even stronger with inefficient phages (lower *φ*). The analytic condition in Eq. (9) (blue dashed line) predicts well the therapy outcome transition. The two triangle markers correspond to the *I* and *φ* inferred in [16] on a model without infected bacteria classes from *in vivo* experiments, for wild-type (immunocompetent, upward light-blue triangle) and MyD88⁻/⁻ (immunodeficient, downward dark-green triangle) mice. The green (grey) dashed line represents Eq. (9) with a 3-fold increase (decrease) in the concentration of bacteria at the beginning of the therapy. The bigger the bacteria inoculum (worse infection), the harder it is to clear the infection using phage with moderate to low *φ*. Simulation parameters are reported in Table 2.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012793.g005

in an immunocompetent cohort but failed in immunocompromised mice. The green and grey dashed lines in Fig 5 show the quantitative impact of bacteria concentration at the treatment start, as they correspond to Eq. (9) with respectively a 3-fold increase and decrease in $B_0$. We note that the highest $B_0$ scenario brings the success to failure transition very close to the parameters that saw successful treatment in immunocompetent mice in [16].

## Robust benefits of therapeutic phage cocktails within an *in vivo* infection model, modulating immune responses and phage-bacteria interactions

Finally, we address the sensitivity of phage cocktail efficacy given variation in innate immune efficiency and phage life history traits. When adding a second phage to an *in vivo* infection model with a responsive modulated immune system and structured nonlinear phage-bacteria interactions, as described in Section “Multi-phage treatment model, modulating the
<figure><img src="image_9.png" /><figcaption>**Fig 6. Efficient phage cocktails improve therapy success in an in vivo model of immunocompromised hosts.** Numerical simulations of the model in Eq. (4) varying *I* and *φ*, with *p* = 1. The color map represents the density of bacteria in the last part of the numerical simulations, where the white region denotes bacteria clearance. Phage cocktails can drive bacteria to extinction in immunocompromised hosts provided that phages have a high enough adsorption rate *φ*. The dashed lines show Eq. (11) for *p* = 0, 0.5 and 1 (blue, yellow, green), which agrees well with the simulation results. The black dashed line shows Eq. (10), representing therapy success when bacteria do not develop phage resistance. The two triangle markers correspond to the *I* and *φ* inferred in [16] from *in vivo* experiments, for wild-type (immunocompetent, upward light-blue triangle) and MyD88<sup>-/-</sup> (immunodeficient, downward dark-green triangle) mice. Simulation parameters are reported in Table 2.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1012793.g006

evolutionary interactions between bacteria and phage strains during infections”, we predict that treatment would succeed if (see Sec. 8 in S1 Text)

$$I > \max\left(\frac{r - \phi P_c}{\kappa}\left(1 + \frac{B_0}{K_D}\right), \frac{r - p\phi P_c/2}{\kappa}\right). \qquad (11)$$

The first term in the maximum, the same as in Eq. (9), ensures that susceptible bacteria are killed, whereas the second term is the condition under which resistant types can be kept in check. Fig 6 shows the infection clearance pattern as a function of *I* and *φ* for simulations of model (4) with *p* = 1. The dashed lines show Eq. (11) for *p* = 0, 0.5 and 1 (blue, yellow, green). The black line represents a scenario where bacteria do not develop phage resistance, in which case the success condition is given just by Eq. (10). Our simulations agree with the analytical prediction in Eq. (11), as shown in *p* = 1 in Fig 6 (complemented by S2 Fig showing the numerical population dynamics with *I* = 10⁵, for three values of *φ* close to the sharp outcome transition) and in S1 Fig for the other cases.
The simulations findings in Fig 6 show that the treatment with two phages improves the therapeutic outcome in immunocompromised hosts if the adsorption rate is high enough, confirming the qualitative picture presented in Section “Model of a phage combination therapy against multiple pathogenic bacteria strains”. Eq. (11) quantifies how the therapy success depends on phage-bacteria evolutionary interactions through variation in the cross-phage resistance parameter *p*. There is a critical $\phi_c = \frac{2r}{pP_c}$ above which therapy succeeds when $I = 0$. We interpret this finding to mean that higher *p*, i.e. using therapeutic phages with more complementary phenotypes, improves the therapeutic outcome. In contrast, when *p* = 0 bacteria can evolve simultaneous resistance to both phages leading to the same result as in the single-phage treatment. Hence it is crucial to pre-select phages against which it is less likely that bacteria could evolve combined resistance, for instance making sure that they do not share a common receptor target, or that receptors targets are not pleiotropically linked [43], like in our LD experiment scoring a potential cocktail of PAK_P1 and LUZ19v against *P. aeruginosa*. When *p* = 1 and the phage combination covers all bacteria resistant mutants, there is a region of parameters, between the black and the green lines in Fig 6, where therapy fails due to phage resistance despite the phage cocktail treatment. Here the only way to improve the treatment efficacy would be to use phages with stronger lytic features.

Interestingly, when $p < 2/3$ and $P_2$ is a generalist phage used alone, the maximum phage killing term in Eq. (4) against the most resistant bacteria type would be larger than when using two phages for the same value of *p* (see Sec. 8 in S1 Text for the derivation). Therefore, it would be better to deliver $P_2$ alone than together with another phage. S3 Fig shows that, in the special case $p = 0.5$, a single $P_2$ produces the same numerical results as a phage cocktail with $p = 1$, and indeed gives better outcomes than two phages with $p = 0.5$, as expected from the theory. This result is conditioned on the specific assumption we made in this model that phages compete to infect bacteria. Nevertheless it showcases an extreme situation where using a phage cocktail could be worse than a single phage, suggesting that it is crucial to quantify experimentally the effect of phage combinations in the system of interest. Even in this scenario, sub-optimal treatment with two phages with intermediate *p* still yields much better therapy outcomes in immunocompromised hosts than a treatment composed of a single or multiple phages that do not target all bacteria mutants in the pathogen population, as evident when comparing the dashed blue and yellow lines in Fig 6. These results can be generalized to cases when the life history traits of the two phage differ from one another (see Sec. 9 in S1 Text). Altogether, this analysis reinforces the importance of selecting a combination of efficient phage against target bacteria when designing phage cocktails for treating infections in immunocompromised hosts.

## Discussion

In this work we analyzed tripartite population dynamics models of single-phage and phage combination treatments to clear infections by multiple strains of bacteria capable of evolving phage resistance in immunomodulated hosts. In doing so we extended previous work that considered the impacts of phage therapy when bacteria were exclusively susceptible to phage [26] by considering the potential combined use of phages with complementary modes of infection. This assumption is supported by new fluctuation test experiments in which *P. aeruginosa* was unlikely to randomly acquire double resistance mutations to phage LUZ19v and PAK_P1. Beginning with monophage treatment and extending this to phage cocktails, we find that immunophage synergy underlies the curative treatment of bacterial infections given sufficiently efficient phage. Notably, the use of phage cocktails can extend the range of immunocompromised conditions in which phage therapy can clear a pathogen. Finally, we
extended core theoretical findings to a realistic *in vivo* modeling contexts, showing the robustness of immunophage synergy given variation in immune state, phage adsorption rates, and asymmetry in phage effectiveness within cocktails. Our analytical results quantify the importance of a prompt infection treatment and of selecting phages that are highly effective against the target pathogens including potential resistant mutants [69], especially when dealing with immunocompromised hosts. We extrapolate therapy outcome predictions around parameters inferred in *in vivo* lung infections by *P. aeruginosa* in immunomodulated mice [16], showing that moderate adsorption rate variations in any employed phage can have drastic effects on therapy outcomes, potentially making the difference between a successful and a failing treatment in experimental applications.

Our theoretical exploration of immunophage synergy builds on a set of assumptions that come with caveats to be addressed in future work - broadly speaking we categorize this in terms of simplifications in our representation of the immune system, phage infection and administration dynamics, and the evolutionary relationship between phage and bacteria. First of all, we consider a relatively simple impact of the host immune system on pathogens, through quantitative features that have been proposed in past infection models [16,53,56]. The first part of this study, focusing on non-responsive immunity, is consistent with a signaling deficient immune system as in the case of myeloid differentiation primary response gene 88-deficient mice (MyD88⁻/⁻) [16]. The predictions in this limit may also be tested in *ex vivo* experiments mimicking infections in the lungs [70] by inoculating therapeutic phage and a fixed amount of immune cells. The innate immune responses considered in the second part of our work focuses on neutrophils, the first immune barrier against invading pathogens [71], while neglecting the other components of the innate immune response and the adaptive immune response altogether. Notably, phage therapy can be hindered by inhibition of phage via immune responses [54] and by the reduction of circulating infectious phage by macrophages [25]. While our analysis focuses on the density of pathogens in the lungs as a proxy for phage therapy outcomes, other factors can contribute to the clinical condition of the host. For instance, infected patients can die from sepsis - i.e., organ failure due to the dysregulation of the host response to an infection [72]. A model accounting for clinical factors affecting host health would require a hazard function relating several immune system components and other indicators to the likelihood of patient survival, which is an ideal topic for future studies focusing on alternative *in vivo* contexts with dysregulated immune activity. It will be important to increase our quantitative understanding of the impact of the complex immune dynamics arising when treating infections with phage, including those triggering sepsis.

The overall modeling structure used here adopts an implicit view of complex spatial processes taking place during phage treatment of respiratory lung infections [73]. Although it is possible to include effective, nonlinear interaction terms to mimic spatial complexity [16], moving forward it will be paramount to evaluate the explicit impact of spatial structure on quantitative phage-pathogens dynamics, for instance leveraging *ex vivo* technologies [74,75], so that future models can incorporate spatial components [76]. In this study we also do not address the impact of treatment timing on therapy outcomes. A theoretical work applying control theory on a phage cocktail model suggested that additional treatment improvements may be possible by optimizing the timing and distribution of phage titers [77]. Such optimum would depend on the pathogen population composition and likelihood of resistance mutations, as well as on the quantitative features of phage-bacteria interactions during an infection, such as the functional shape of phage adsorption profiles. Simultaneous administration has been shown to outperform sequential treatments as controlling the pathogen population size right away reduces the chances of multi-resistance [41,78], even though the generality of these
results is unclear. It will be important to integrate empirically constrained population dynam-
ics models of local infections with pharmacokinetics parameters describing the likelihood and
delay of delivering phage in the desired infected tissue [24].

Finally, analysis of phage cocktail impacts here assumes relatively simple evolutionary
interactions between phage and bacteria such that bacteria cannot evolve resistance to both
phages at the same time. Our choice is motivated by previous works that suggested combi-
nations of phages that target different bacteria receptors in order to improve treatment effi-
cacy [42,43,55]. This hypothesis is supported by fluctuation test findings in which *P. aerugi-
nosa* did not randomly develop simultaneous resistance to LUZ19v and PAK_P1 (see Sec. 3
in S1 Text). It is important to note that although a wild type population susceptible to both
phages need not necessarily evolve a double resistant mutant immediately, this does not pre-
clude the potential that a population already selected to persist given exposure to one of the
two phages may evolve to become doubly resistant. The current study shows that phage cock-
tails can perform robustly even when modulating rules regarding phage complementarity.
It will be crucial to further explore more complex eco-evolutionary processes that can arise
between pathogenic bacteria and therapeutic phage during the course of an infection, whether
in an acute or chronic infection context.

Despite these caveats, the combined use of experiments, simulations and theory provides
guidance on the expected range of phage therapeutic efficacy whether using monophage
or phage cocktail treatments. Building upon earlier findings [16], our framework provides
testable predictions on the quantitative impact of different modes of tripartite phage-bacteria-
immune interactions on therapeutic outcomes over a range of host immune conditions and
phage life history traits, highlighting the success of single-phage therapy in synergy with a
strong enough immune system and the benefit of phage combination therapy in immuno-
compromised hosts. Importantly, the parameters inferred in [16] fall right at the boundary
between treatment failure and success in immunodeficient hosts, which could result in dras-
tic differences in the infection outcomes given small variations in phage and immune features.
Such sensitivity makes it crucial to inform model development with *in vitro* and *in vivo* data
to improve therapeutic design.

More broadly, this work represents a further step towards quantitatively addressing the
impact of evolutionary considerations on phage therapy outcomes. Here we showcased a spe-
cific example of how we can harness the evolutionary potential of phages to develop phage
cocktails that target specific strains of multi-drug resistant bacteria, facing the challenge posed
by the evolution of phage resistance [43,44]. Whether we seek to exploit *in vitro* phage train-
ing [31,44] or evolutionary trade-offs [32–37], in the future it will be essential to integrate
further experimental evidence into quantitative models [27], tackling different aspects of
phage-bacteria evolutionary interactions during therapy. Predictive models that integrate
the tripartite interactions among pathogenic bacteria, therapeutic phages, and the eukaryotic
host, along with evolutionary processes, will be crucial for designing effective phage treatment
strategies both in the near and long term.

**Supporting information**

**S1 Text. Supplementary information.**

(PDF)

**S1 Fig. Numerical simulations confirm the cocktail success condition.** Numerical simula-
tions of the model in Eq. (4) varying *I* and *φ*, with *p* = 0 (A), *p* = 0.5 (B) and without allowing
for phage resistance (C). The colormap represents the density of bacteria in the last part of
the numerical simulations. The dashed line represents the prediction in Eq. (11) for *p* = 0 in
A) and 0.5 in B) (blue and yellow respectively). The black dashed line in C) shows Eq. (10), representing therapy success when bacteria do not develop phage resistance. Simulation parameters are reported in Table 2.
(TIF)

**S2 Fig. Phage cocktail treatment dynamics.** Numerical simulations of the model in Eq. (4), with $p = 1$ and $I = 10^5$. Each panel shows the population dynamics for a different value of $\phi$ close to the sharp therapy success transition in Eq. (11). A) $\phi = 8 \cdot 10^{-8}$ falls within the red region of the phase diagram in Fig 6, with bacteria saturating to carrying capacity $K_C$. B) With $\phi = 10^{-7}$, right at the theoretical transition, bacteria coexist with phage at an intermediate density. C) Just above the transition at $\phi = 1.1 \cdot 10^{-7}$ bacteria are cleared, corresponding to the white region in the phase diagrams. The other simulation parameters are reported in Table 2.
(TIF)

**S3 Fig. A generalist phage works better alone than in combination.** Numerical simulations of the model in Eq. (4) varying $I$ and $\phi$, with a single phage $P_2$ and $p = 0.5$. The colormap represents the density of bacteria in the last part of the numerical simulations. In this extreme case, a single-phage treatment works better than two phages. The dashed lines compare Eq. (11) for $p = 0.5$ with two phages (yellow) with the transition for a single generalist phage (lightblue, equivalent to Eq. (11) with $p = 1$). Simulation parameters are reported in Table 2.
(TIF)

**S4 Fig. Phage cocktail treatments with different viral traits produce different outcomes.** Numerical simulations of the model in Eq. (4) with $\phi_2 = \sigma\phi_1$, varying $I$ and $\phi_1$, for $\sigma = 1/5$ (A) and $\sigma = 5$ (B). The color map represents the density of bacteria in the last part of the numerical simulations. A) While phage cocktails can still drive bacteria to extinction in immunocompromised hosts, a low $\sigma$ can jeopardize the treatment outcome. B) A high $\sigma$ improves therapy, especially in immunocompetent hosts. The dashed lines show Eq. (S19) for $\sigma = 1/5$, 1 and 5 (lighter green, green, teal). The simulations agree reasonably well with the analytical approximation. The two triangle markers correspond to the $I$ and $\phi$ inferred in [16] from in vivo experiments, for wild-type (immunocompetent, upward light-blue triangle) and MyD88$^{-/}$ (immunodeficient, downward dark-green triangle) mice. Simulation parameters are reported in Table 2.
(TIF)

**S5 Fig. Simulation dynamics with different $\phi$s.** Numerical simulations of the model in Eq. (4) varying $I$ and $\phi_1$, for $\sigma = 1/5$ $\phi_1 = 5 \cdot 10^{-7}$ and $I = 10^4$. The rest of the simulation parameters are reported in Table 2.
(TIF)

**Acknowledgments**

We thank Jeremy Seurat and Rogelio Rodriguez-Gonzalez for insights in the development of the in vivo phage therapy model.

**Author contributions**

**Conceptualization:** Jacopo Marchi, Joshua S. Weitz.

**Data curation:** Jacopo Marchi, Chau Nguyen Ngoc Minh.

**Formal analysis:** Jacopo Marchi.
**Funding acquisition:** Laurent Debarbieux, Joshua S. Weitz.

**Investigation:** Jacopo Marchi, Chau Nguyen Ngoc Minh.

**Methodology:** Jacopo Marchi.

**Project administration:** Laurent Debarbieux, Joshua S. Weitz.

**Software:** Jacopo Marchi.

**Supervision:** Laurent Debarbieux, Joshua S. Weitz.

**Validation:** Jacopo Marchi, Joshua S. Weitz.

**Visualization:** Jacopo Marchi, Chau Nguyen Ngoc Minh.

**Writing – original draft:** Jacopo Marchi.

**Writing – review & editing:** Jacopo Marchi, Laurent Debarbieux, Joshua S. Weitz.

## References

1. Suttle CA. Marine viruses–major players in the global ecosystem. Nat Rev Microbiol. 2007;5(10):801–12. https://doi.org/10.1038/nrmicro1750 PMID: 17853907

2. Koskella B, Brockhurst MA. Bacteria-phage coevolution as a driver of ecological and evolutionary processes in microbial communities. FEMS Microbiol Rev. 2014;38(5):916–31. https://doi.org/10.1111/1574-6976.12072 PMID: 24617569

3. Luria SE, Delbrück M. Mutations of Bacteria from Virus Sensitivity to Virus Resistance. Genetics. 1943;28(6):491–511. https://doi.org/10.1093/genetics/28.6.491 PMID: 17247100

4. Labrie SJ, Samson JE, Moineau S. Bacteriophage resistance mechanisms. Nat Rev Microbiol. 2010;8(5):317–27. https://doi.org/10.1038/nrmicro2315 PMID: 20348932

5. Hampton HG, Watson BNJ, Fineran PC. The arms race between bacteria and their phage foes. Nature. 2020;577(7790):327–36. https://doi.org/10.1038/s41586-019-1894-8 PMID: 31942051

6. Bratbak G, Heldal M, Norland S, Thingstad TF. Viruses as partners in spring bloom microbial trophodynamics. Appl Environ Microbiol. 1990;56(5):1400–5. https://doi.org/10.1128/aem.56.5.1400-1405.1990 PMID: 16348190

7. Wichels A, Biel SS, Gelderblom HR, Brinkhoff T, Muyzer G, Schütt C. Bacteriophage diversity in the North Sea. Appl Environ Microbiol. 1998;64(11):4128–33. https://doi.org/10.1128/AEM.64.11.4128-4133.1998 PMID: 9797256

8. Breitbart M, Rohwer F. Here a virus, there a virus, everywhere the same virus?. Trends Microbiol. 2005;13(6):278–84. https://doi.org/10.1016/j.tim.2005.04.003 PMID: 15936660

9. Kortright KE, Chan BK, Koff JL, Turner PE. Phage Therapy: A Renewed Approach to Combat Antibiotic-Resistant Bacteria. Cell Host Microbe. 2019;25(2):219–32. https://doi.org/10.1016/j.chom.2019.01.014 PMID: 30763536

10. Neu HC. The crisis in antibiotic resistance. Science. 1992;257(5073):1064–73. https://doi.org/10.1126/science.257.5073.1064 PMID: 1509257

11. Ten threats to global health in 2019. World Health Organization. 2019.

12. Antimicrobial Resistance Collaborators. Global burden of bacterial antimicrobial resistance in 2019: a systematic analysis. Lancet. 2022;399(10325):629–55. https://doi.org/10.1016/S0140-6736(21)02724-0 PMID: 35065702

13. Kutter EM, Kuhl SJ, Abedon ST. Re-establishing a place for phage therapy in western medicine. Future Microbiol. 2015;10(5):685–8. https://doi.org/10.2217/fmb.15.28 PMID: 26000644

14. Cohan FM, Zandi M, Turner PE. Broadscale phage therapy is unlikely to select for widespread evolution of bacterial resistance to virus infection. Virus Evol. 2020;6(2):veaa060. https://doi.org/10.1093/ve/veaa060 PMID: 33365149

15. Smith HW, Huggins MB. Effectiveness of phages in treating experimental Escherichia coli diarrhoea in calves, piglets and lambs. J Gen Microbiol. 1983;29(8):2659–75. https://doi.org/10.1099/00221287-8-129-8-2659 PMID: 6355391

16. Roach DR, Leung CY, Henry M, Morello E, Singh D, Di Santo JP, et al. Synergy between the host immune system and bacteriophage is essential for successful phage therapy against an acute respiratory pathogen. Cell Host Microbe. 2017;22(1):38-47.e4. https://doi.org/10.1016/j.chom.2017.06.018 PMID: 28704651
17. Jennes S, Merabishvili M, Soentjens P, Pang KW, Rose T, Keersebilck E, et al. Use of bacteriophages in the treatment of colistin-only-sensitive *Pseudomonas aeruginosa* septicaemia in a patient with acute kidney injury-a case report. Crit Care. 2017;21(1):129.
https://doi.org/10.1186/s13054-017-1709-y PMID: 28583189

18. Schooley RT, Biswas B, Gill JJ, Hernandez-Morales A, Lancaster J, Lessor L, et al. development and use of personalized bacteriophage-based therapeutic cocktails to treat a patient with a disseminated resistant Acinetobacter baumannii infection. Antimicrob Agents Chemother. 2017;61(10):e00954-17. https://doi.org/10.1128/AAC.00954-17 PMID: 28807909

19. Leitner L, Ujmajuridze A, Chanishvili N, Goderdzishvili M, Chkonia I, Rigvava S, et al. Intravesical bacteriophages for treating urinary tract infections in patients undergoing transurethral resection of the prostate: a randomised, placebo-controlled, double-blind clinical trial. Lancet Infect Dis. 2021;21(3):427-36. https://doi.org/10.1016/S1473-3099(20)30330-3 PMID: 32949500

20. Jault P, Leclerc T, Jennes S, Pirnay JP, Que Y-A, Resch G, et al. Efficacy and tolerability of a cocktail of bacteriophages to treat burn wounds infected by *Pseudomonas aeruginosa* (PhagoBurn): a randomised, controlled, double-blind phase 1/2 trial. Lancet Infect Dis. 2019;19(1):35-45. https://doi.org/10.1016/S1473-3099(18)30482-1 PMID: 30292481

21. Blasco L, López-Hernández I, Rodríguez-Fernández M, Pérez-Florido J, Casimiro-Soriguer CS, Djebara S, et al. Case report: Analysis of phage therapy failure in a patient with a *Pseudomonas aeruginosa* prosthetic vascular graft infection. Front Med (Lausanne). 2023;101199657. https://doi.org/10.3389/fmed.2023.1199657 PMID: 37275366

22. Pirnay J, Djebara S, Steurs G, Griselain J, Cochez C, De Soir S. Personalized bacteriophage therapy outcomes for 100 consecutive cases: a multicentre, multinational, retrospective observational study. Nature Microbiology. 2024-20.

23. Marchi J, Zborowsky S, Debarbieux L, Weitz JS. The dynamic interplay of bacteriophage, bacteria and the mammalian host during phage therapy. iScience. 2023;26(2):106004.
https://doi.org/10.1016/j.isci.2023.106004 PMID: 36818291

24. Dabrowska K. Phage therapy: What factors shape phage pharmacokinetics and bioavailability? Systematic and critical review. Med Res Rev. 2019;39(5):2000-25.
https://doi.org/10.1002/med.21572 PMID: 30887551

25. Zborowsky S, Seurat J, Balacheff Q, Minh CNN, Tit'ecat M, Evrard E. Macrophage-induced reduction of bacteriophage density limits the efficacy of in vivo pulmonary phage therapy. bioRxiv. 2024.

26. Leung CYJ, Weitz JS. Modeling the synergistic elimination of bacteria by phage and the innate immune system. J Theor Biol. 2017;429241-52. https://doi.org/10.1016/j.jtbi.2017.06.037 PMID: 28668337

27. Cairns BJ, Timms AR, Jansen VAA, Connerton IF, Payne RJH. Quantitative models of in vitro bacteriophage-host dynamics and their application to phage therapy. PLoS Pathog. 2009;5(1):e1000253. https://doi.org/10.1371/journal.ppat.1000253 PMID: 19119417

28. Tiwari BR, Kim S, Rahman M, Kim J. Antibacterial efficacy of lytic Pseudomonas bacteriophage in normal and neutropenic mice models. J Microbiol. 2011;49(6):994-9.
https://doi.org/10.1007/s12275-011-1512-4 PMID: 22203564

29. Duerkop BA, Huo W, Bhardwaj P, Palmer KL, Hooper LV. Molecular Basis for Lytic Bacteriophage Resistance in Enterococci. mBio. 2016;7(4):e01304-16. https://doi.org/10.1128/mBio.01304-16 PMID: 27578757

30. Oechslin F. Resistance Development to Bacteriophages Occurring during Bacteriophage Therapy. Viruses. 2018;10(7):351. https://doi.org/10.3390/v10070351 PMID: 29966329

31. Borin JM, Avrani S, Barrick JE, Petrie KL, Meyer JR. Coevolutionary phage training leads to greater bacterial suppression and delays the evolution of phage resistance. Proc Natl Acad Sci U S A. 2021;118(23):e2104592118. https://doi.org/10.1073/pnas.2104592118 PMID: 34083444

32. Chan BK, Sistrom M, Wertz JE, Kortright KE, Narayan D, Turner PE. Phage selection restores antibiotic sensitivity in MDR *Pseudomonas aeruginosa*. Sci Rep. 2016;6:626717.
https://doi.org/10.1038/srep26717 PMID: 27225966

33. Filippov AA, Sergueev KV, He Y, Huang X-Z, Gnade BT, Mueller AJ, et al. Bacteriophage-resistant mutants in Yersinia pestis: identification of phage receptors and attenuation for mice. PLoS One. 2011;6(9):e25486. https://doi.org/10.1371/journal.pone.0025486 PMID: 21980477

34. Laanto E, Bamford JKH, Laakso J, Sundberg L-R. Phage-driven loss of virulence in a fish pathogenic bacterium. PLoS One. 2012;7(12):e53157.
https://doi.org/10.1371/journal.pone.0053157 PMID: 23308090
35. Sumrall ET, Shen Y, Keller AP, Rismondo J, Pavlou M, Eugster MR, et al. Phage resistance at the cost of virulence: Listeria monocytogenes serovar 4b requires galactosylated teichoic acids for InIB-mediated invasion. PLoS Pathog. 2019;15(10):e1008032.
https://doi.org/10.1371/journal.ppat.1008032 PMID: 31589660

36. Gurney J, Pradier L, Griffin JS, Gougat-Barbera C, Chan BK, Turner PE, et al. Phage steering of antibiotic-resistance evolution in the bacterial pathogen, *Pseudomonas aeruginosa*. Evol Med Public Health. 2020;2020(1):148-57. https://doi.org/10.1093/emph/eoaa026 PMID: 34254028

37. Gordillo Altamirano F, Forsyth JH, Patwa R, Kostoulias X, Trim M, Subedi D, et al. Bacteriophage-resistant Acinetobacter baumannii are resensitized to antimicrobials. Nat Microbiol. 2021;6(2):157-61. https://doi.org/10.1038/s41564-020-00830-7 PMID: 33432151

38. Gaborieau B, Delattre R, Adiba S, Clermont O, Denamur E, Ricard J-D, et al. Variable fitness effects of bacteriophage resistance mutations in Escherichia coli: implications for phage therapy. J Virol. 2024;98(10):e0111324. https://doi.org/10.1128/jvi.01113-24 PMID: 39213164

39. Kortright KE, Chan BK, Turner PE. High-throughput discovery of phage receptors using transposon insertion sequencing of bacteria. Proc Natl Acad Sci U S A. 2020;117(31):18670-9.
https://doi.org/10.1073/pnas.2001888117 PMID: 32675236

40. Tanji Y, Shimada T, Fukudomi H, Miyanaga K, Nakai Y, Unno H. Therapeutic use of phage cocktail for controlling Escherichia coli O157:H7 in gastrointestinal tract of mice. Journal of Bioscience and Bioengineering. 2005;100(3):280-7.

41. Hall AR, De Vos D, Friman V-P, Pirnay J-P, Buckling A. Effects of sequential and simultaneous applications of bacteriophages on populations of *Pseudomonas aeruginosa* in vitro and in wax moth larvae. Appl Environ Microbiol. 2012;78(16):5646-52. https://doi.org/10.1128/AEM.00757-12 PMID: 22660719

42. Borin JM, Lee JJ, Gerbino KR, Meyer JR. Comparison of bacterial suppression by phage cocktails, dual-receptor generalists, and coevolutionarily trained phages. Evol Appl. 2022;16(1):152-62.
https://doi.org/10.1111/eva.13518 PMID: 36699129

43. Wright RCT, Friman V-P, Smith MCM, Brockhurst MA. Functional diversity increases the efficacy of phage combinations. Microbiology (Reading). 2021;167(12):001110.
https://doi.org/10.1099/mic.0.001110 PMID: 34850676

44. Wang M, Wei J, Jiang L, Jiang L, Zhang J, He X, et al. Coevolutionary phage training and Joint application delays the emergence of phage resistance in *Pseudomonas aeruginosa*. Virus Evolution. 2023;9(2):vead067. https://doi.org/10.1093/ve/vead067

45. Gillet-Markowska A, Louvel G, Fischer G. bz-rates: A Web Tool to Estimate Mutation Rates from Fluctuation Analysis. G3 (Bethesda). 2015;5(11):2323-7. https://doi.org/10.1534/g3.115.019836 PMID: 26338660

46. Hamon A, Ycart B. Statistics for the Luria-Delbrück distribution. Electron J Statist. 2012;6(none).
https://doi.org/10.1214/12-ejs711

47. Lenski RE, Levin BR. Constraints on the Coevolution of Bacteria and Virulent Phage: A Model, Some Experiments, and Predictions for Natural Communities. The American Naturalist. 1985;125(4):585-602. https://doi.org/10.1086/284364284364

48. Levin BR, Bull JJ. Population and evolutionary dynamics of phage therapy. Nat Rev Microbiol. 2004;2(2):166-73. https://doi.org/10.1038/nrmicro822 PMID: 15040264

49. de Kievit TR, Iglewski BH. Bacterial quorum sensing in pathogenic relationships. Infect Immun. 2000;68(9):4839-49. https://doi.org/10.1128/IAI.68.9.4839-4849.2000 PMID: 10948095

50. Sully EK, Malachowa N, Elmore BO, Alexander SM, Femling JK, Gray BM, et al. Selective chemical inhibition of agr quorum sensing in Staphylococcus aureus promotes host defense with minimal impact on resistance. PLoS Pathog. 2014;10(6):e1004174.
https://doi.org/10.1371/journal.ppat.1004174 PMID: 24945495

51. Costerton JW, Stewart PS, Greenberg EP. Bacterial biofilms: a common cause of persistent infections. Science. 1999;284(5418):1318-22. https://doi.org/10.1126/science.284.5418.1318
PMID: 10334980

52. Gellatly SL, Hancock REW. *Pseudomonas aeruginosa*: new insights into pathogenesis and host defenses. Pathog Dis. 2013;67(3):159-73. https://doi.org/10.1111/2049-632X.12033 PMID:
23620179

53. Drusano GL, Fregeau C, Liu W, Brown DL, Louie A. Impact of burden on granulocyte clearance of bacteria in a mouse thigh infection model. Antimicrob Agents Chemother. 2010;54(10):4368-72.
https://doi.org/10.1128/AAC.00133-10 PMID: 20516275

54. Hodyra-Stefaniak K, Miernikiewicz P, Drapała J, Drab M, Jończyk-Matysiak E, Lecion D, et al. Mammalian Host-Versus-Phage immune response determines phage fate in vivo. Sci Rep.
2015;514802. https://doi.org/10.1038/srep14802 PMID: 26440922
55. Naknaen A, Samernate T, Wannasrichan W, Surachat K, Nonejuie P, Chaikeeratisak V. Combination of genetically diverse *Pseudomonas phages* enhances the cocktail efficiency against bacteria. Sci Rep. 2023;13(1):8921. https://doi.org/10.1038/s41598-023-36034-2 PMID: 37264114

56. Summers C, Rankin SM, Condliffe AM, Singh N, Peters AM, Chilvers ER. Neutrophil kinetics in health and disease. Trends Immunol. 2010;31(8):318-24. https://doi.org/10.1016/j.it.2010.05.006 PMID: 20620114

57. Guo B, Abdelraouf K, Ledesma KR, Chang KT, Nikolaou M, Tam VH. quantitative impact of neutrophils on bacterial clearance in a murine pneumonia model. Antimicrobial Agents and Chemotherapy. 2011;55(10):4601-5. https://doi.org/10.1128/aac.00508-11

58. Hurtado PJ, Kirosingh AS. Generalizations of the "Linear Chain Trick": incorporating more flexible dwell time distributions into mean field ODE models. J Math Biol. 2019;79(5):1831-83. https://doi.org/10.1007/s00285-019-01412-w PMID: 31410551

59. Weitz J. Quantitative viral ecology: dynamics of viruses and their microbial hosts. 2015.

60. Mitarai N, Brown S, Sneppen K. Population dynamics of phage and bacteria in spatially structured habitats using phage λ and *Escherichia coli*. J Bacteriol. 2016;198(12):1783-93. https://doi.org/10.1128/JB.00965-15 PMID: 27068593

61. Eriksen RS, Mitarai N, Sneppen K. Sustainability of spatially distributed bacteria-phage systems. Sci Rep. 2020;10(1):3154. https://doi.org/10.1038/s41598-020-59635-7 PMID: 32081858

62. Dominguez-Mirazo M, Harris JD, Demory D, Weitz JS. Accounting for cellular-level variation in lysis: implications for virus-host dynamics. mBio. 2024;15(8):e0137624. https://doi.org/10.1128/mbio.01376-24 PMID: 39028198

63. Nabergoj D, Modic P, Podgornik A. Effect of bacterial growth rate on bacteriophage population growth rate. Microbiologyopen. 2018;7(2):e00558. https://doi.org/10.1002/mbo3.558 PMID: 29195013

64. Nguyen TVP, Wu Y, Yao T, Trinh JT, Zeng L, Chemla YR, et al. Coinfecting phages impede each other's entry into the cell. Curr Biol. 2024;34(13):2841-2853.e18. https://doi.org/10.1016/j.cub.2024.05.032 PMID: 38878771

65. Abedon ST. Kinetics of phage-mediated biocontrol of bacteria. Foodborne Pathog Dis. 2009;6(7):807-15. https://doi.org/10.1089/fpd.2008.0242 PMID: 19459758

66. Guo Y, Chen P, Lin Z, Wang T. Characterization of Two *Pseudomonas aeruginosa* Viruses vB_PaeM_SCUT-S1 and vB_PaeM_SCUT-S2. Viruses. 2019;11(4):318. https://doi.org/10.3390/v11040318 PMID: 30939832

67. Zborowsky S, Balacheff Q, Theodorou I, Kane R, Delattre R, Weitz JS, et al. A nanoluciferase-encoded bacteriophage illuminates viral infection dynamics of *Pseudomonas aeruginosa* cells. ISME Commun. 2024;4(1):ycae105. https://doi.org/10.1093/ismeco/ycae105 PMID: 39296778

68. Mägälie A, Marantos A, Schwartz DA, Marchi J, Lennon JT, Weitz JS. Phage infection fronts trigger early sporulation and collective defense in bacterial populations. 2024. https://doi.org/10.1101/2024.05.22.595388

69. Romeyer Dherbey J, Bertels F. The untapped potential of phage model systems as therapeutic agents. Virus Evol. 2024;10(1):veae007. https://doi.org/10.1093/ve/veae007 PMID: 38361821

70. Rossy T, Distler T, Meirelles LA, Pezoldt J, Kim J, Talà L, et al. *Pseudomonas aeruginosa* type IV pili actively induce mucus contraction to form biofilms in tissue-engineered human airways. PLoS Biol. 2023;21(8):e3002209. https://doi.org/10.1371/journal.pbio.3002209 PMID: 37527210

71. Kolaczkowska E, Kubes P. Neutrophil recruitment and function in health and inflammation. Nat Rev Immunol. 2013;13(3):159-75. https://doi.org/10.1038/nri3399 PMID: 23435331

72. Singer M, Deutschman CS, Seymour CW, Shankar-Hari M, Annane D, Bauer M, et al. The third international consensus definitions for sepsis and septic shock (Sepsis-3). JAMA. 2016;315(8):801-10. https://doi.org/10.1001/jama.2016.0287 PMID: 26903338

73. Lourenço M, Chaffringeon L, Lamy-Besnier Q, Pédron T, Campagne P, Eberl C, et al. The spatial heterogeneity of the gut limits predation and fosters coexistence of bacteria and bacteriophages. Cell Host Microbe. 2020;28(3):390-401.e5. https://doi.org/10.1016/j.chom.2020.06.002 PMID: 32615090

74. Barr JJ, Auro R, Furlan M, Whiteson KL, Erb ML, Pogliano J, et al. Bacteriophage adhering to mucus provide a non-host-derived immunity. Proc Natl Acad Sci U S A. 2013;110(26):10771-6. https://doi.org/10.1073/pnas.1305923110 PMID: 23690590

75. Joiner KL, Baljon A, Barr J, Rohwer F, Luque A. Impact of bacteria motility in the encounter rates with bacteriophage in mucus. Sci Rep. 2019;9(1):16427. https://doi.org/10.1038/s41598-019-52794-2 PMID: 31712565

76. Rodriguez-Gonzalez RA, Balacheff Q, Debarbieux L, Marchi J, Weitz JS. Metapopulation model of phage therapy of an acute *Pseudomonas aeruginosa* lung infection. bioRxiv. 2024.024-01.
77. Li G, Leung CY, Wardi Y, Debarbieux L, Weitz JS. Optimizing the timing and composition of therapeutic phage cocktails: a control-theoretic approach. Bull Math Biol. 2020;82(6):75. https://doi.org/10.1007/s11538-020-00751-w PMID: 32533350

78. Wright RCT, Friman V-P, Smith MCM, Brockhurst MA. Resistance evolution against phage combinations depends on the timing and order of exposure. mBio. 2019;10(5):e01652-19. https://doi.org/10.1128/mBio.01652-19 PMID: 31551330