
==== Front
bioRxiv
BIORXIV
bioRxiv
2692-8205
Cold Spring Harbor Laboratory

39282405
10.1101/2024.09.07.611814
preprint
1
Article
Multi-strain phage induced clearance of bacterial infections
Marchi Jacopo *Department of Biology, University of Maryland, College Park, MD, USA

Ngoc Minh Chau Nguyen Institut Pasteur, Université Paris Cité, CNRS UMR6047, Bacteriophage Bacterium Host, Paris, France and Sorbonne Université, Collège Doctoral, Paris, France

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

Weitz Joshua S. †Department of Biology, University of Maryland, College Park, MD USA
Department of Physics, University of Maryland, College Park, MD USA
University of Maryland Institute for Health Computing, North Bethesda, MD and Institut de Biologie, École Normale Supérieure, Paris, France

* jmarchi@umd.edu
jsweitz@umd.edu - Former affiliation: School of Biological Sciences, Georgia Institute of Technology, Atlanta, GA, USA
07 9 2024
2024.09.07.611814https://creativecommons.org/licenses/by-nc-nd/4.0/ This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which allows reusers to copy and distribute the material in any medium or format in unadapted form only, for noncommercial purposes only, and only so long as attribution is given to the creator.
nihpp-2024.09.07.611814.pdf
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 phage (i.e., ‘cocktails’) have 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.
==== Body
pmcI. 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–29]. 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 [30]. 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 Pa; hence joint use of phage OMKO1 and antibiotics can lead to therapeutic success [31]. 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, 32–37]. 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 [38–41]. Cocktails with phages targeting different cell receptors have higher efficacy against P. aeruginosa [42]. 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 [43].

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.

II. METHODS

A. 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 SI Sec. II. 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 [44], which learns the model parameters via the Generating Function estimator from [45] (see SI Sec. III for more details on the inference).

B. Mathematical models of phage therapy

1. Summary

We study two different models of single- and double-phage therapy of bacterial infections in immunomodulated hosts leading to four different scenarios. Figure 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] (Figure 1, A1). We first analyze a single-phage treatment model (Figure 1, B1), then we add a second phage to the cocktail to address its impact on treatment outcome (Figure 1, B2 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 (Figure 1, A2). 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 (Figure 1, B2 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.

2. 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: (1) dS(t)dt=rS-d˜SS+R+μrR-S-κIS1+S+RKD-ΦSP,dR(t)dt=rR-d˜RS+R+μrS-R-κIR1+S+RKD,dP(t)dt=βΦSP-ωP,

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 KC=rd˜. They mutate between the susceptible and resistant types, and the susceptible type can be infected by phage with adsorption rate ϕ. At the same time phage lyse susceptible bacteria producing β virions per lysis event, and decay at a constant rate ω. 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 KD, 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 ϕ). Table I reports the parameters meaning and values.

3. 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 P1 and P2 (Figure 1, B2 top). The pathogen types include the wild-type infecting bacteria S, susceptible to both phages, which can in turn mutate to a type R1 resistant to P1 or to a type R2 resistant to P2: (2) dS(t)dt=rS-d˜SBtot+μrR1+R2-S-κIS1+BtotKD-ΦSP1-ΦSP2,dR1(t)dt=rR1-d˜R1Btot+μrS-R1-κIR11+BtotKD-ΦR1P2,dR2(t)dt=rR2-d˜R2Btot+μrS-R2-κIR21+BtotKD-ΦR2P1,dP1(t)dt=βΦS+R2P1-ωP1,dP2(t)dt=βΦS+R1P2-ωP2.

We initially assume that the two phages have the same adsorption rate ϕ, 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 [41, 42, 54]. 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 viceversa), 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.

4. 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] (Figure 1, A2): (3) dS(t)dt=rS-d˜SBtot+μrR-S-κiS1+BtotKD-SΦFP1,dR(t)dt=rR-d˜RBtot+μrS-R-κiR1+BtotKD,dE11(t)dt=ΦSFP1-LηE1,1dE1i∈{2…L}(t)dt=LηE1i-1-LηE1i,dP1(t)dt=βLηE1L-ΦSFP1-ωP1,di(t)dt=αi1-iIBtotBtot+KN.

The innate immune cells i are recruited in the lungs at a maximum rate per capita α 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 Btot with half saturation constant KN. These saturation profiles in immune recruitment within infected tissues are supported by empirical evidence that neutrophils can only be produced up to a limit [59], and that innate immune responses saturate as a function of bacterial loads [53], as already mentioned in the previous section.

Additionally, we structure the bacteria population introducing L infected bacteria stages, E1i∈{1…L}, between phage adsorption and lysis to model a finite phage infection time distributed as an Erlang distribution with average 1η, an approach applied in several models of epidemiological dynamics and phage-bacteria interactions [57, 60, 61]. Finally, as proposed in [16], we add a saturation term in phage adsorption FP1=P11+P1Pc, as if phage adsorbed into bacteria through Michaelis-Menten reaction kinetics [62, 63], producing a phage-adsorption profile similar to Monod growth [64]. Table II provides a list of parameter definitions and their values.

5. 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 P1 and P2, 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 R1 resistant to P1 or to a type R2 that can be infected by P1. Compared to section IIB3, now we introduce a parameter p continuously modulating the evolutionary relations between the two phages and the target bacteria, in an asymmetric fashion. P2 can either infect R1 with probability p, or R2 with probability 1-p. So when p=1 the two phages are complementary in infecting each other’s resistant type with rate ϕ, and Ri 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 P2 capable of infecting either bacterium, but with lower specificity to either type since the sum of the effective adsorption rates to all bacteria is always ϕ. The equations for this model are: (4) dS(t)dt=rS-d˜SBtot+μrR1+R2-S-κiS1+BtotKD-SΦFP1-SΦFP2,dR1(t)dt=rR1-d˜R1Btot+μrS-R1-κiR11+BtotKD-R1pΦFP2,dR2(t)dt=rR2-d˜R2Btot+μrS-R2-κiR21+BtotKD-R2ΦFP1-R2(1-p)ΦFP2,dE11(t)dt=ΦSFP1+ΦR2FP1-LηE11,dE1iϵ{2…L}(t)dt=LηE1i-1-LηE1i,dE21(t)dt=ΦSFP2+(1-p)ΦR2FP2+pϕR1FP2-LηE21,dE2iϵ{2…L}(t)dt=LηE2i-1-LηE2i,dP1(t)dt=βLηE1L-ΦSFP1-ΦR2FP1-ωP1,dP2(t)dt=βLηE2L-ΦSFP2-(1-p)ΦR2FP2-pϕR1FP2-ωP2,di(t)dt=αi1-iIBtotBtot+KN,

where E1i denotes a bacterium infected by P1 (at infection stage i), and E2i denotes a bacterium infected by P2. We start by considering phages with the same life history traits, and then relax this assumption studying a case where P1 adsorbs at a rate ϕ1 while P2 adsorbs at a rate ϕ2. In this case we study a scenario where FP1=P11+P1+P2Pc, 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.

III. RESULTS

A. 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 IIB2. In the absence of phage, the model admits a stable fix point BIS and an unstable one BIU for bacteria densities. After the addition of phage, whenever bacteria are driven below the infectious dose BIU 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 BP=ωβΦ< KDκIKCrKD-12 (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 I0: (5) I>I0=d˜KDκ1+ωβΦKD2,

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 SI Sec. III B). The mathematical condition BIU>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 BIU before the resistant bacteria proliferate from rare to densities above this threshold. When bacteria densities are depleted below BIU, 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 BIU>0, crossing the threshold I>I0 drives a transition from therapeutic failure to success as phage drive the susceptible bacteria below BIU 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 KD is lower than bacteria carrying capacity KC 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 immune strength attaining BIU>0, necessary for immunophage synergy (see SI Sec. IIIB for the derivation): (6) I>Ib=rκ.

If phage drive the susceptible bacteria below BIU fast enough, the simultaneous satisfaction of the conditions in Eqs. (5) and (6), or more compactly (7) I>IS=maxI0,Ib,

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 ϕ. 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 SI Sec. III A). 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).

B. 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.

Figure 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 ϕ). This is expected as the added phage will lyse the bacteria resistant to the other one, as shown in the population dynamics in Figure 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 Figure 4B).

We propose a heuristic derivation of the phase space of immunophage synergy, as detailed in SI Sec. III C. 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 ϕ>ϕc, where ϕc satisfies the transcendental equation: (8) 2ϕcP0r=2+lnϕcP0r+lnωϕcΩβTωr.

Here P0 is the initial phage density. Ω is a numerical factor ensuring bacteria exponential decay before extinction (see SI Sec. III C), which needs to satisfy Ω≫1 and Ω<ωϕcβT.Ω affects the equation only logarithmically, and ϕc decreases with Ω. The condition in Eq. (8) gives the lowest ϕ above which phage can drive bacteria below the extinction threshold T before bacteria resume growth (due to phage decay at rate ω). The leading factors determining the phage cocktail success are the initial rate of phage killing ϕP0 and bacteria growth r, which drive the transition to phage-bacteria coexistence. The dashed vertical green line in Figure 4A shows ϕc for Ω=10. We observe that the numerical simulations agree with the heuristic prediction Eq. (8) when I≪Ib. Note that when phage are not very efficient ϕ≪ϕc the presence of resistance is not the major treatment failure driver, therefore Eq. (5) is an upper bound 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 bacterial mutants are present, insofar as sufficiently effective phage – according to the quantitative condition ϕ>ϕc – are utilized.

C. 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 assuming 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 [62, 63]), incorporating 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 ϕ. We find a qualitatively similar pattern as in Section III A, 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] (blue diamonds). 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 (upper diamonds), 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 IIIA 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 BIU. To progress further in the analysis we assume that phage infection dynamics are fast compared to other timescales, disregarding the dynamics of infected bacteria E1i∈{1…L}, which corresponds exactly to the model in [16]. Assuming also a scenario in which phage are so abundant that P≫Pc, and that the immune cells dynamics relax fast to I, we find (details of the derivation in SI Sec. III D) a new immune threshold IS that constrains the host immune capacity for which the single-phage therapy succeeds: (9) I>IS=maxI1,Ib,

with (10) I1=r-ϕPcκ1+B0KD.

As in Section IIIA 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 κ regardless of the new model ingredients. Eq. (9) depends on the infecting bacteria density at the time of therapy administration B0, highlighting the importance of prompt intervention when treating the infection. The dash blue line in Fig. 5 shows Eq. 9 as a function of ϕ. The analytical result agrees with the numerical results of the general model in Eq. 3 and with the findings of [16] where therapy worked 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 B0. We note that the highest B0 scenario brings the success to failure transition very close to the parameters that saw successful treatment in immunocompetent mice in [16].

D. 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 IIB5, we predict that treatment would succeed if (see SI Sec. IV) (11) I>maxr-ϕPcκ1+B0KD,r-pϕPc/2κ.

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 and in Fig. S1 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 IIB3. 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 ϕc=2rpPc 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 [42], 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 P2 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 SI Sec. IV for the derivation). Therefore, it would be better to deliver P2 alone than together with another phage. Fig. S2 shows that, in the special case p=0.5, a single P2 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 SI Sec. V). 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.

IV. DISCUSSION

In this work, analyzed simple tripartite population dynamic 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 phage 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 [65], 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, 59]. 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 [66] 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 [67], while neglecting the other components of the innate immune response and the adaptive immune response altogether. Notably, inhibition of phage via immune responses [46] and by the reduction of circulating infectious phage by macrophages [25]. In the future it will be important to increase our quantitative understanding of the impact of the complex immune dynamics arising during infections.

The overall modeling structure used here adopts an implicit view of complex spatial processes taking place during phage treatment of respiratory lung infections [68]. 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 [69, 70], so that future models can incorporate spatial components [71]. 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 [72]. 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 [40, 73], even though the generality of these results is unclear. In the future it will be important to integrate empirically constrained population dynamics 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 combinations of phages that target different bacteria receptors in order to improve treatment efficacy [41, 42, 54]. This hypothesis is supported by fluctuation test findings in which P. aeruginosa did not randomly develop simultaneous resistance to LUZ19v and PAK_P1 (see SI Section III). 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 preclude the potential for a population already selected to persist given exposure to one of the two phages may evolve to become double resistant. The current study shows that phage cocktails can perform robustly even when modulating rules regarding phage complementarity. In the future 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 provide 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 immunocompromised hosts. Importantly, the parameters inferred in [16] fall right at the boundary between treatment failure and success in immunodeficient hosts, which could result in drastic 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 specific 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 [42, 43]. Whether we seek to exploit in vitro phage training [30, 43] or evolutionary trade-offs [31–36], in the future it will be essential to integrate further experimental evidence into quantitative models, 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.

Supplementary Material

Supplement 1

Acknowledgements.

The work was supported by the National Institutes of Health (R01 AI146592 to JSW, LD). JSW was supported, in part, by the Chaires Blaise Pascal program of the Ile-de-France region. We thank Jeremy Seurat and Rogelio Rodriguez-Gonzalez for insights in the development of the in vivo phage therapy model.

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 A1 shows the interactions in the simple phage therapy model in a signaling deficient host introduced in Section IIB2, whereas Subpanel A2 sketches the addition of more complex dynamics in the model introduced in Section IIB4, 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 B1 sketches the single-phage treatment structure used in Sections IIB2 and IIB4, whereas B2 summarizes the cross-infection structure in the phage-cocktail treatment models introduced in Sections IIB3 and IIB5 (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(⋅)∝˙⋅+C.

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 I0 (left panel) or above (center, right panels), with phage (left, center) or without (right). Increasing the immune system strength above I0, with BIU>0, phage can drive susceptible bacteria below BIU 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 I with ϕ=1011g/(hPFU),I=1.5⋅107 and 2.9⋅107 cells/g respectively in the left and the other two panels.

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 I.

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 ϕc obtained solving Eq. (8) with Ω=10 The black triangles correspond to the parameters yielding the population dynamics in panels B),C),D) that show the evolution of phage and bacteria strains during the simulated infection. Simulation parameters are reported in Table I.

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. 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 blue diamonds correspond to the I and ϕ inferred in [16] from in vivo experiments on a model without infected bacteria classes. 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 II.

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. 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 blue diamonds correspond to the I and ϕ inferred in [16] from in vivo experiments. Simulation parameters are reported in Table II.

TABLE I: Model parameters.

Most values are taken from [26]. The mutation rate is inferred from the fluctuation test (Section II A). For the viral decay ω we elect a mid-value between the parameter reported in [26, 46] 1h-1 and the one inferred in [16] 0.07h-1.

Variable	Meaning	Value	Unit	
r	Maximum bacterial growth rate	1	h-1	
d˜	Bacteria density-driven death rate	10−10	g/(h CFU)	
μ	Mutation probability per replication	4.3 · 10−8		
ϕ	Phage adsorption rate	[10−12, 10−6]	g/(h PFU)	
β	Burst size	100	(PFU/cell)	
ω	Viral decay	0.4	h-1	
κ	Immune killing rate	8.2 · 10−8	g/(h CFU)	
I	Immune cells density	[104, 109]	(immune cell)/g	
KD	Cell shielding density	107	(CFU)/g	

TABLE II: Model parameters.

Most values are taken from [16]. The mutation rate is inferred via a fluctuation test (Section II). The lysis rate 2h-1 is compatible with reported lysis times of phages infecting P. aeruginosa [55, 56]. The number of infected stages L is the same as several previous models of phage-bacteria dynamics [57, 58].

Variable	Meaning	Value	Unit	
r	Maximum bacterial growth rate	0.75	h-1	
d˜	Bacteria density-driven death rate	7.5 · 10−9	g/(h CFU)	
μ	Mutation probability per replication	4.3 · 10−8		
ϕ	Phage adsorption rate	[10−9, 10−6]	g/(h PFU)	
Pc	Phage adsorption saturation density	1.5 · 107	(PFU)/g	
η	Lysis rate	2	h-1	
L	Number of infected stages	10		
β	Burst size	100	(PFU/cell)	
ω	Viral decay	0.07	h-1	
κ	Immune killing rate	8.2 · 10−8	g/(h CFU)	
i	Immune cells capacity	[104, 108]	(immune cell)/g	
KD	Bacteria shielding density from immune action	4.1 · 107	(CFU)/g	
KN	Bacteria density when immune recruitment is half its maximum	107	(CFU)/g	
α	ximum immune cells recruitment rate	0.97	h-1
==== Refs
[1] Suttle CA (2007) Marine viruses — major players in the global ecosystem. Nature Reviews Microbiology 5 :801–812.17853907
[2] Koskella B , Brockhurst MA (2014) Bacteria–phage coevolution as a driver of ecological and evolutionary processes in microbial communities. FEMS Microbiology Reviews 38 :916–931.24617569
[3] Luria SE , Delbrück M (1943) Mutations of bacteria from virus sensitivity to virus resistance. Genetics 28 :491–511.17247100
[4] Labrie SJ , Samson JE , Moineau S (2010) Bacteriophage resistance mechanisms. Nature Reviews Microbiology 8 :317–327.20348932
[5] Hampton HG , Watson BNJ , Fineran PC (2020) The arms race between bacteria and their phage foes. Nature 577 :327–336.31942051
[6] Bratbak G , Heldal M , Norland S , Thingstad TF (1990) Viruses as partners in spring bloom microbial trophodynamics. Applied and Environmental Microbiology 56 :1400–1405.16348190
[7] Wichels A , (1998) Bacteriophage diversity in the north sea. Applied and environmental microbiology 64 :4128–4133.9797256
[8] Breitbart M , Rohwer F (2005) Here a virus, there a virus, everywhere the same virus? Trends in microbiology 13 :278–284.15936660
[9] Kortright KE , Chan BK , Koff JL , Turner PE (2019) Phage therapy: a renewed approach to combat antibiotic-resistant bacteria. Cell host & microbe 25 :219–232.30763536
[10] Neu HC (1992) The crisis in antibiotic resistance. Science 257 :1064–1073.1509257
[11] WHO (2019) Ten threats to global health in 2019. (https://www.who.int/news-room/spotlight/ten-threats-to-global-health-in-2019).
[12] Murray C , Ikuta K , Sharara F , Moore C (2022) Global burden of bacterial antimicrobial resistance in 2019: a systematic analysis. The Lancet 399 :P629–655 Publisher: Elsevier.
[13] Kutter EM , Kuhl SJ , Abedon ST (2015) Re-establishing a place for phage therapy in western medicine. Future Microbiology 10 :685–688 26000644
[14] Cohan FM , Zandi M , Turner PE (2020) Broadscale phage therapy is unlikely to select for widespread evolution of bacterial resistance to virus infection. Virus Evolution 6 :veaa060.33365149
[15] Smith HW , Huggins MB (1983) Effectiveness of phages in treating experimental escherichia coli diarrhoea in calves, piglets and lambs. Microbiology 129 :2659–2675.
[16] Roach DR , (2017) Synergy between the host immune system and bacteriophage is essential for successful phage therapy against an acute respiratory pathogen. Cell host & microbe 22 :38–47.28704651
[17] Jennes S , (2017) Use of bacteriophages in the treatment of colistin-only-sensitive pseudomonas aeruginosa septicaemia in a patient with acute kidney injury-a case report. Critical Care 21 .
[18] Schooley RT , (2017) Development and use of personalized bacteriophage-based therapeutic cocktails to treat a patient with a disseminated resistant acinetobacter baumannii infection. Antimicrobial Agents and Chemotherapy 61 :10.1128/aac.00954-17
[19] Leitner L , (2020) Intravesical bacteriophages for treating urinary tract infections in patients undergoing transurethral resection of the prostate: a randomised, placebo-controlled, double-blind clinical trial. The Lancet. Infectious diseases 21 .
[20] Jault P , (2019) 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. The Lancet. Infectious diseases 19 1 :35–45.30292481
[21] Blasco L , (2023) Case report: Analysis of phage therapy failure in a patient with a pseudomonas aeruginosa prosthetic vascular graft infection. Frontiers in Medicine.
[22] Pirnay JP , (2024) Personalized bacteriophage therapy outcomes for 100 consecutive cases: a multicentre, multinational, retrospective observational study. Nature Microbiology pp 1–20.
[23] Marchi J , Zborowsky S , Debarbieux L , Weitz JS (2023) The dynamic interplay of bacteriophage, bacteria and the mammalian host during phage therapy. iScience 26 :106004.36818291
[24] Dabrowska K (2019) Phage therapy: What factors shape phage pharmacokinetics and bioavailability? systematic and critical review. Medicinal research reviews 39 :2000–2025.30887551
[25] Zborowsky S , (2024) Macrophage-induced reduction of bacteriophage density limits the efficacy of in vivo pulmonary phage therapy. (https://www.biorxiv.org/content/early/2024/01/16/2024.01.16.575879) Under review.
[26] Leung J , Weitz J (2017) Synergistic elimination of bacteria by phage and the innate immune system. Journal of Theoretical Biology 429 .
[27] Tiwari B , Kim S , Rahman M , Kim J (2011) Antibacterial efficacy of lytic pseudomonas bacteriophage in normal and neutropenic mice models. Journal of microbiology (Seoul, Korea) 49 :994–9.22203564
[28] Duerkop BA , Huo W , Bhardwaj P , Palmer KL , Hooper LV (2016) Molecular basis for lytic bacteriophage resistance in enterococci. mBio 7 :10.1128/mbio.01304-16.
[29] Oechslin F (2018) Resistance development to bacteriophages occurring during bacteriophage therapy. Viruses 10 .
[30] Borin JM , Avrani S , Barrick JE , Petrie KL , Meyer JR (2021) Coevolutionary phage training leads to greater bacterial suppression and delays the evolution of phage resistance. Proceedings of the National Academy of Sciences 118 :e2104592118.
[31] Chan BK , (2016) Phage selection restores antibiotic sensitivity in mdr pseudomonas aeruginosa. Scientific reports 6 :26717.27225966
[32] Filippov AA , (2011) Bacteriophage-resistant mutants in yersinia pestis: identification of phage receptors and attenuation for mice. PloS one 6 :e25486.21980477
[33] Laanto E , Bamford JKH , Laakso J , Sundberg LR (2013) Phage-driven loss of virulence in a fish pathogenic bacterium. PLOS ONE 7 :1–7.
[34] Sumrall ET , (2019) Phage resistance at the cost of virulence: Listeria monocytogenes serovar 4b requires galactosylated teichoic acids for inlb-mediated invasion. PLOS Pathogens 15 :1–29.
[35] Gurney J , (2020) Phage steering of antibiotic-resistance evolution in the bacterial pathogen, pseudomonas aeruginosa. Evolution, medicine, and public health 2020 :148–157.34254028
[36] Gordillo Altamirano F , (2021) Bacteriophage-resistant acinetobacter baumannii are resensitized to antimicrobials. Nature microbiology 6 :157–161.
[37] Gaborieau B , (2024) Variable fitness effects of bacteriophage resistance mutations in escherichia coli: implications for phage therapy. Journal of Virology.
[38] Kortright KE , Chan BK , Turner PE (2020) High-throughput discovery of phage receptors using transposon insertion sequencing of bacteria. Proceedings of the National Academy of Sciences 117 :18670–18679.
[39] Tanji Y , (2005) Therapeutic use of phage cocktail for controlling escherichia coli o157:h7 in gastrointestinal tract of mice. Journal of bioscience and bioengineering 100 3 :280–7.16243277
[40] Hall AR , Vos DD , Friman V , Pirnay JP , Buckling A (2012) Effects of sequential and simultaneous applications of bacteriophages on populations of pseudomonas aeruginosa in vitro and in wax moth larvae. Applied and Environmental Microbiology 78 :5646 – 5652.22660719
[41] Borin JM , Lee JJ , Gerbino KR , Meyer JR (2023) Comparison of bacterial suppression by phage cocktails, dual-receptor generalists, and coevolutionarily trained phages. Evolutionary Applications 16 :152–162.36699129
[42] Wright RC , Friman VP , Smith MC , Brockhurst MA (2021) Functional diversity increases the efficacy of phage combinations. Microbiology 167 .
[43] Wang M , (2023) Coevolutionary phage training and Joint application delays the emergence of phage resistance in Pseudomonas aeruginosa. Virus Evolution 9 :vead067.38089014
[44] Gillet-Markowska A , Louvel G , Fischer G (2015) bz-rates: A Web Tool to Estimate Mutation Rates from Fluctuation Analysis. G3 Genes—Genomes—Genetics 5 :2323–2327.26338660
[45] Hamon A , Ycart B (2012) Statistics for the luria-delbrück distribution. Electronic Journal of Statistics 6 .
[46] Hodyra-Stefaniak K , (2015) Mammalian host-versus-phage immune response determines phage fate in vivo. Scientific reports 5 :14802.26440922
[47] Lenski RE , Levin BR (1985) Constraints on the coevolution of bacteria and virulent phage: a model, some experiments, and predictions for natural communities. The American Naturalist 125 :585–602.
[48] Levin BR , Bull JJ (2004) Population and evolutionary dynamics of phage therapy. Nature Reviews Microbiology 2 :166–173.15040264
[49] De Kievit TR , Iglewski BH (2000) Bacterial quorum sensing in pathogenic relationships. Infection and immunity 68 :4839–4849.10948095
[50] Sully EK , (2014) Selective chemical inhibition of agr quorum sensing in staphylococcus aureus promotes host defense with minimal impact on resistance. PLoS pathogens 10 :e1004174.24945495
[51] Costerton JW , Stewart PS , Greenberg EP (1999) Bacterial biofilms: a common cause of persistent infections. science 284 :1318–1322.10334980
[52] Gellatly SL , Hancock RE (2013) Pseudomonas aeruginosa: new insights into pathogenesis and host defenses. Pathogens and disease 67 :159–173.23620179
[53] Drusano G , Fregeau C , Liu W , Brown D , Louie A (2010) Impact of burden on granulocyte clearance of bacteria in a mouse thigh infection model. Antimicrobial agents and chemotherapy 54 :4368–4372.20516275
[54] Naknaen A , (2023) Combination of genetically diverse pseudomonas phages enhances the cocktail efficiency against bacteria. Scientific Reports 13 .
[55] Guo Y , Chen P , Lin Z , Wang T (2019) Characterization of two pseudomonas aeruginosa viruses vb_paem_scut-s1 and vb_paem_scut-s2. Viruses 11 .
[56] Zborowsky S , (2024) A nanoluciferase-encoded bacteriophage illuminates viral infection dynamics of Pseudomonas aeruginosa cells. ISME Communications p ycae 105 .
[57] Mitarai N , Brown S , Sneppen K (2016) Population dynamics of phage and bacteria in spatially structured habitats using phage λ and escherichia coli. Journal of bacteriology 198 :1783–1793.27068593
[58] Măgălie A , (2024) Phage infection fronts trigger early sporulation and collective defense in bacterial populations. bioRxiv.
[59] Summers C , (2010) Neutrophil kinetics in health and disease. Trends in immunology 31 :318–324.20620114
[60] Hurtado P , Kirosingh A (2019) Generalizations of the ‘linear chain trick’: incorporating more flexible dwell time distributions into mean field ode models. Journal of Mathematical Biology 79 .
[61] Weitz J (2015) Quantitative Viral Ecology: Dynamics of Viruses and Their Microbial Hosts.
[62] Nabergoj D , Modic P , Podgornik A (2017) Effect of bacterial growth rate on bacteriophage population growth rate. MicrobiologyOpen 7 :e00558.29195013
[63] Nguyen TVP , (2024) Coinfecting phages impede each othrer’s entry into the cell. Current Biology 34 :2841–2853.38878771
[64] Abedon ST (2009) Kinetics of phage-mediated biocontrol of bacteria. Foodborne pathogens and disease 6 :807–815.19459758
[65] Romeyer Dherbey J , Bertels F (2024) The untapped potential of phage model systems as therapeutic agents. Virus Evolution 10 :veae007.38361821
[66] Rossy T , (2023) Pseudomonas aeruginosa type iv pili actively induce mucus contraction to form biofilms in tissueengineered human airways. PLOS Biology 21 :1–32.
[67] Kolaczkowska E , Kubes P (2013) Neutrophil recruitment and function in health and inflammation. Nature reviews immunology 13 :159–175.
[68] Lourenço M , (2020) The spatial heterogeneity of the gut limits predation and fosters coexistence of bacteria and bacteriophages. Cell Host & Microbe 28 :390–401.e5.
[69] Barr JJ , (2013) Bacteriophage adhering to mucus provide a non-host-derived immunity. Proceedings of the National Academy of Sciences 110 :10771–10776.
[70] Joiner KL , Baljon A , Barr J , Rohwer F , Luque A (2019) Impact of bacteria motility in the encounter rates with bacteriophage in mucus. Scientific reports 9 :16427.31712565
[71] Rodriguez-Gonzalez RA , Balacheff Q , Debarbieux L , Marchi J , Weitz JS (2024) Metapopulation model of phage therapy of an acute pseudomonas aeruginosa lung infection. bioRxiv pp 2024–01.
[72] Li G , Leung CY , Wardi Y , Debarbieux L , Weitz JS (2020) Optimizing the timing and composition of therapeutic phage cocktails: a control-theoretic approach. Bulletin of mathematical biology 82 :1–29.
[73] Wright RC , Friman VP , Smith MC , Brockhurst MA (2019) Resistance evolution against phage combinations depends on the timing and order of exposure. MBio 10 :10–1128.
[74] Moir DT , Di M , Opperman T , Schweizer HP , Bowlin TL (2007) A high-throughput, homogeneous, bioluminescent assay for pseudomonas aeruginosa gyrase inhibitors and other dna-damaging agents. Journal of Biomolecular Screening 12 :855–864 17644773
[75] Debarbieux L , (2010) Bacteriophages Can Treat and Prevent Pseudomonas aeruginosa Lung Infections. The Journal of Infectious Diseases 201 :1096–1104.20196657
[76] Ferran A , (2022) The selection of antibiotic- and bacteriophage-resistant pseudomonas aeruginosa is prevented by their combination. Microbiology Spectrum 10 .
[77] Boulanger P (2009) Purification of Bacteriophages and SDS-PAGE Analysis of Phage Structural Proteins from Ghost Particles, eds Clokie MR, Kropinski AM (Humana Press, Totowa, NJ), pp 227–238.
[78] Lea D , Coulson C (1949) The distribution of the numbers of mutants in bacterial populations. Journal of genetics 49 :264–285.24536673
[79] Virtanen P , (2020) SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods 17 :261–272.32015543
[80] Dormand JR , Prince PJ (1980) A family of embedded runge-kutta formulae. Journal of Computational and Applied Mathematics 6 :19–26.
[81] Gil A , Segura J , Temme NM (2007) Numerical Methods for Special Functions (Society for Industrial and Applied Mathematics).
