Rapid resistance evolution against phage cocktails

When bacteria are treated with multiple antibiotics simultaneously, resistance is exceedingly unlikely to evolve. In stark contrast, resistance against multiple phages frequently arises during therapy. Why does resistance against multi-phage cocktails evolve so easily? Using a mathematical model, we show how the bacterial evolutionary dynamics and phage replicative dynamics uniquely intertwine, facilitating the rapid evolution of multi-phage resistance. As different phages replicate and become inhibitory at varying time points, bacteria can sequentially acquire resistance rather than simultaneously – increasing the chance of multi-resistance by orders of magnitude. Additionally, we identify a regime where multi-phage resistance is robustly prevented. Our findings provide a framework for the rational design of phage cocktails to minimize resistance development.

In stark contrast to this prediction, resistance against all components of a multi-phage cocktail is frequently observed during therapy 2,5 , suggesting that multi-phage resistance evolves much more rapidly than predicted by current evolutionary theory.
Here, we use a deliberately minimalistic mathematical model to investigate how phage replication dynamics affect bacterial resistance evolution.To set the stage, we first model a bacterial population predated by a single type of phage.In line with previous phage replication models 8 , we consider a generic process where phages bind and burst bacteria, after which new phage particles are released (Methods).Using realistic parameters (Suppl.Table 1), we find that the bacterial population initially grows seemingly unperturbed by the presence of the phage until, after several hours, the bacterial population suddenly collapses due to phage predation (Fig. 1a).The time point of collapse shortens as the burst size (number of phages produced per bacterium) increases (blue to red).In all cases, the collapse is followed by a relapse of growth by resistant bacteria, which naturally emerge during the growth of the bacteria.These simulations closely match experimentally observed dynamics of bacterial populations undergoing phage predation 9 .
We next investigate the dynamics of a bacterial population undergoing predation by two different phages simultaneously and find strikingly different behavior depending on the ratio of burst sizes between both phages.We focus first on the dynamics in the absence of the emergence of multi-phage resistant bacteria (Methods).When phages are combined with a similar burst size, we only observe a single and definitive population collapse (Fig. 1b, dark curves).By contrast, when phages with sufficiently different burst sizes are combined, we observe an initial collapse and regrowth of a resistant mutant, followed by a second and final population collapse (Fig. 1b, blue and red curves).
The message of this Brief Report consists of two parts: 1) that and why there are these two qualitatively different regimes for phage combinations, and 2) how these regimes massively impact the chance of escape mutants against the phage cocktail.To make this concrete, we introduce the phage selection pressure (Methods Eq. 5), which is proportional to the phage concentration.In essence, the phage concentration at which the selection pressure equals 1 is comparable to the minimal inhibitory concentration (MIC) of an antibiotic, with the important distinction that the phage concentration is strongly time-dependent due to the replicative nature of phages.
The time-resolved selection pressure (Fig. 1c) explains the difference between the single-peaked and double-peaked regimes.Initially, the selection pressure grows super-exponentially as phages are self-replicative, with a rate that depends on the bacterial density -which in turn grows exponentially at first (Fig. 1b).When either of the two phages crosses the MIC (Fig. 1c, dashed line), the bacterial population collapses due to predation by that phage.If both phages reach the MIC before the population has collapsed, the selection pressure of both phages is exerted simultaneously and only a single and definitive collapse is observed.In contrast, when the phage with the lower burst size has not yet reached the MIC, its further amplification is halted until bacteria resistant to the high-burst size phage emerge and reach a large population size.A second collapse is observed when the low-burst size phage now also crosses the MIC threshold -in this case, the selective pressure of both phages is exerted sequentially.
While simultaneous double-resistance is hard to evolve, sequential resistance evolution is much more likely.We show that despite the simultaneous application of both phages at t=0h, the selection pressure can still be exerted sequentially depending on the burst size (Fig. 1c).In the single-collapse regime of simultaneous selection pressure, the chance that a double-escape mutant arises (Methods Eq. 7) is extremely small (<1 in a million, Fig 1d).In contrast, in the double-collapse regime, the chance of escape mutants rapidly increases to near-certainty as double resistance can be acquired step-wise.Therefore, we conclude that the evolution of multi-phage resistance is only prevented within a window of opportunity where both phages are selective simultaneously (Fig. 2).
We next investigate what determines the size of this window of opportunity, where the chance of escape mutants is small, by systematically varying the model's parameters and initial conditions (Suppl.Fig. 1).We find that the window size is hardly influenced by the phage kinetic parameters, nor by the initial conditions.Instead, the size of the window of opportunity is almost entirely determined by the bacterial growth rate (Suppl.Fig. 1): the generation of a double-resistant mutant is extremely unlikely only when the collapse times of the two separate phages coincide within ~1 bacterial doubling time of each other.
Lastly, we use this convenient criterion to investigate the chance that two arbitrary phages have similar enough infection parameters to operate in the window of opportunity where the selection pressure is exerted simultaneously.We investigate a recently published set of experiments, where the collapse time has been measured for many different phages on the same bacterial host 9 and find that 9.1% of the phage pairs have a collapse time that happens within 0.5 h (roughly the cell doubling time) of each other (Methods).We then calculate the chance that any two phages are within this window as a function of the size of a cocktail, where phages are randomly picked from the above-mentioned pool (Methods Eq. 8).We find that this chance increases roughly linearly with the number of phages and that, for this specific library in vitro, ~10 phages are required to robustly prevent resistance evolution (Fig. 1e).This linear scaling is qualitatively different from the exponential scaling of chance of resistance for antibiotics 4 , and explains why multi-phage resistance is so frequently observed during therapy.
What can we learn from the results of this mathematical model for designing phage therapy cocktails?At first glance, it seems like phage resistance can be prevented by in vitro screening different phages and/or their initial titers such that they cause a collapse at similar time points (Suppl.Fig. 1).However, the pharmacokinetics of phages inside the body 10 are substantially different from the dynamics in a test tube.Thus, we are skeptical that the right in vivo regime can be reliably predicted from in vitro data.Instead, the chance of being in the right regime can be substantially increased by further expanding the cocktail composition (Fig. 1e).Off-the-shelf cocktails typically already have a large size, some including more than 20 different phages 12 , but these are not tailored to the bacterium of interest so the number of effective phages is likely much lower for most patients.The frequent use of these large cocktails does suggest that many different phages can be safely administered simultaneously 12 .Personalized cocktails are currently much smaller (2-5 different phages 2 ), but larger phage banks in combination with more high throughput phage screening methods 13 would allow for substantial expansion, thereby providing a robust, resistance-proof therapy.Alternatively, resistance against some phages causes tradeoffs 7 , such as decreased virulence 14 , or increased sensitivity to antibiotics 2 or other phages 3 .In these cases, resistance evolution does not equal therapy failure.Indeed, some researchers propose including phage resistance tradeoffs as a criterion for picking therapeutic phages 15 .Our work shows the inherent difficulty in preventing phage resistance evolution and therefore emphasizes the clinical importance of investigating and leveraging evolutionary tradeoffs.1b-c as a function of the burst size ratio between phages 1 and 2. Only when the phages exert a simultaneous selection pressure, there is a strong suppression of double-resistance evolution.Panels a-d use the same color coding for the burst size (right).e) As more phages are combined in a cocktail, the chance increases that any combination of two phages applies a simultaneous selection pressure.Here we calculate that chance (Methods Eq. 8) for a public library of phages 16 using previously published data 9 .Fig. 2 Schematic of the resistance evolution pathways against phage cocktails.When a multi-phage cocktail is applied to a bacterial population, the different phages in the cocktail each replicate at their own rate.Consequently, each phage reaches its critical, inhibitory concentration at distinct time points.when multiple phages reach inhibitory concentrations simultaneously (bottom), is the evolution of cocktail resistance prevented.
Suppl.Fig. 1 The size of the window of opportunity mainly depends on the bacterium, not on the phage or the initial conditions.Dual-phage predation assays were performed, using two phages with identical kinetic properties ( , and ) and the concentration of double-resistant mutants evolved at the end was recorded (y-axes, similar to Fig. 1d) for various initial phage concentrations of one phage, keeping the other phage concentration constant.For each dual-phage predation simulation, single-phage predation simulations were performed to measure the difference in collapse time between both phages (x-axis).These simulations were performed from low (light) to high (dark) values of various parameters (text in the graph) in 10 logarithmically spaced steps (Suppl.Table 1).
[cfu] Suppl.Table 1 Parameters used for modeling In Suppl.Fig. 1, we simulate a dual-phage predation assay.Both phages have the same parameters, which we vary one by one over a large range in 10 logarithmically spaced steps, keeping the other parameters fixed.In all of our simulations, we use a resistance mutation rate for both phages of 10 -8 .To tune the difference in collapse time, we vary the initial concentration of one phage in 20 logarithmically spaced steps from 10-10 7 , keeping the initial concentration of the other phage fixed at 10 -4 except for the center bottom panel where we vary the initial concentration of both phages.

Model
To maximize the generality and tractability of our model, we use a deliberately minimalistic model of phage predation.In this model, strictly lytic phages bind and burst, releasing several new phage particles, which can start subsequent infection cycles.Bacteria grow in a steady state and can spontaneously evolve resistance.The number of bacteria that is sensitive to all phages, , grows at a rate and spontaneously evolves resistance against phage at a rate .We include predation by phage , which binds bacteria at a rate : (1) Next, we model the population dynamics of bacteria resistant to phage ( ) and sensitive to all other phages.For simplicity, we do not consider resistance tradeoffs or phage (co-)evolution: (2) Bacteria infected by phage ( ) burst at a rate : ( For each burst, phages are released: From the phage concentration, we calculate a selection pressure as a competition between phage binding and bacterial cell division: (5) Lastly, we count the number of bacterial mutants resistant to both phage and phage that arise in the population: As we are only interested in the chance of a double-resistant mutant arising in the population, we do not model the subsequent growth of these double-resistant mutants.Instead, we calculate the average number of bacterial mutants that arise in the population from Methods Eq. 6. Next, we approximate the chance that a double-mutant evolves, , assuming that the number of mutants follows a Poisson distribution: The set of differential equations was solved numerically using solve_ivp (scipy 1.11.1,Python 3.11.5)and the collapse time of the bacterial populations was detected using find_peaks (scipy 1.11.1,Python 3.11.5).

Chance of simultaneous selection pressure in a phage cocktail
The model shows that whether resistance emerges largely depends on whether two phages exert a sequential or simultaneous selection pressure -which depends on the phage kinetic properties and initial phage concentrations.When more than two phages are combined in a cocktail, the chance increases that any two phages exert a simultaneous selection pressure.We now investigate the chance of being in the simultaneous selection pressure regime, as a function of the cocktail size, when we randomly pick phages from a biobank of phages.Assuming the number of combinations leading to a simultaneous selection pressure is small compared to the total number of combinations ( , then one can compute the chance that at least one pair of the cocktail exerts a simultaneous selection pressure: (8)   The value depends on the composition of the phage biobank.We calculate this parameter for a recently measured dataset 9 , where the collapse time was measured for each phage of the publicly available "Basel collection" 16 , which consists of 68 phages all infecting the same strain of E. Coli (MG1655).We note that this collection is not clinically used, that the bacterial strain is non-pathogenic, and that all data was acquired from in vitro cultures.However, we are not aware of any other data set where the collapse time of many different phages has been measured for the same bacterial strain.In this data set, we find that 9.1% of the phage pairs have a collapse time that happens within 0.5 h (roughly the cell doubling time) of each other -which is roughly the criterion for simultaneous selection pressure (Suppl.Fig. 1).Therefore, in Fig. 1e we plot Methods Eq. 8 with .

Fig. 1
Fig. 1 Phage cocktails exert either a simultaneous or sequential selection pressure, which dictates the chance of multi-phage resistance.a) Bacterial population dynamics under predation by a single type of phage.b) Bacterial population dynamics under predation by two types of phage.Phage 1 has a varied burst size, whereas phage 2 has a fixed burst size of 100.c) The selection pressure of both phages was calculated for the simulations in Fig. 1b.The dashed line is where phage predation and bacterial growth are equal (minimal inhibitory concentration).d) The chance of the emergence of a double-resistant mutant is calculated for the simulations of Fig.1b-cas a function of the burst size ratio between phages 1 and 2. Only when the phages exert a simultaneous selection pressure, there is a strong suppression of double-resistance evolution.Panels a-d use the same color coding for the burst size (right).e) As more phages are combined in a cocktail, the chance increases that any combination of two phages applies a simultaneous selection pressure.Here we calculate that chance (Methods Eq. 8) for a public library of phages16 using previously published data9 .