Towards modeling phage therapy.

1 Introduction

Due to the evolution of bacterial resistance to multiple types of antibiotics, viruses that can infect and kill bacteria (called bacteriophages), are receiving renewed interest as a treatment option for bacterial infections. Recent clinical successes with patients recovering from rampant bacterial infections due to personalized phage therapy, have sparked interest into the design, efficacy, and safety of cocktails of phages that can be given to patients [1,2]. Phage therapy has been used as ‘living antibiotics’ in the Soviet Union over extended periods of time [3], but these recent successes were ‘compassionate use’ cases of life threatening bacterial infections that were treated with personalized cocktails of phages hastily selected based on their capacity to lyse the bacteria causing the infection [4–7]. Unfortunately, the very limited host range of most bacteriophages calls for such a personalized approach, and is a major challenge for generalized treatments and clinical trials [1].

Phage therapy may work for several reasons: (1) phages may simply infect and eliminate all bacteria (which is unlikely due to the rapid evolution of resistance [8,9]), (2) phages may reduce bacterial concentrations to levels where the host immune system regains control [10], (3) by evolving resistance to phages bacteria may regain sensitivity to particular antibiotics [2,4,5,11], (4) by modifying their capsule to prevent infection [7] bacteria may become more sensitive to killing by neutrophils, and (5) bacteria that evolved phage resistance may have such a crippled fitness that their replication rate drops below the rate at which the host immune system eliminates them. It seems virtually impossible to distinguish between these mechanisms on the basis of observational clinical data alone [1,6]. Because we expect many more people with multiple drug resistant (MDR) infections to be treated with phage therapy in the near future [2], we here develop mathematical models implementing these mechanisms, with the aim to employ these models for better interpretation of future data. Such models may help one optimize the design of therapy, and by fitting variations of these models to the data, one may become able to identify the mechanism(s) underlying the control of rampant bacterial infections by phage therapy in individual cases.

The infection of bacteria by phages, and the rapid evolution of resistance by bacteria, has been studied extensivelyin vitro, and it is known that these systems tend to a approach an attractor where susceptible and resistant bacteria co-exist with the phages [8,9]. It was therefore surprising that phage therapy can be successfulin vivo, and by elegantly combining modeling with experiments in mice, Roachet al.showed that the innate immune response, in this case mostly neutrophils, is largely responsible for thein vivosuccess of phage therapy in their experiments [10,9]. Apparently, neutrophils are usually incapable of controlling rampant bacterial infections, and a major effect of phage therapy is to reduce bacterial populations to a level where the killing of bacteria by neutrophils suffices to control bacterial growth. This mechanism requires the existence of a threshold in the bacterial density above which neutrophils can no longer control bacterial growth. Such a threshold is readily brought about in a model where at high bacterial densities (1) the neutrophils approach a maximum density, and (2) the number of bacteria killed per neutrophil per unit of time approaches a maximum, e.g., due to the fact that killing takes time [10,12,13].

Neutrophil densities rapidly increase during inflammation [14]. Because they are released more rapidly from the bone marrow, they can approach a new (higher) steady state level on a time scale of less than a day [15,16]. Thus, there should be two thresholds: (1) the maximum bacterial density that can be controlled by the normal density of neutrophils, e.g., small bacterial invasions that are controlled on a daily basis without triggering additional inflammation, and (2) the much higher threshold above which infections can grow uncontrolled despite the massive killing by neutrophils at their maximal density. In between an innate immune response is required to control the infection before it breaches the upper threshold.

After developing a novel mathematical model having both thresholds, and allowing for several bacterial strains differing in their sensitivity to various phages and antibiotics, we illustrate its usefulness by describing a case of successful phage therapy of a patient suffering from a rampant infection with multi-drug resistant (MDR)Acinetobacter baumanniibacteria [3,4,17]. We present a few scenarios suggesting why this particular phage therapy might have been so successful. This patient was treated with several phage cocktails that were given sequentially, leading to rapid sequential development of bacterial resistance [7]. Most importantly, we find that the level of resistance increases rapidly when sensitive bacteria decline due to the (unsuccessful) phages in the first cocktail, and that this may markedly reduce the efficacy of subsequent cocktails. Thus, similar to lessons learned from treating rapidly evolving viruses like HIV, we find that it is important to ‘hit hard and early’ [18] by providing all phages at the same time, and to increase the ‘genetic barrier’ [19] by providing as many phages as possible.

2 Results

The killing of bacteria by neutrophils has been studied widelyin vitro[20–23] andin vivo[24,25]. Importantly, these data reveal that there are two types of thresholds, one for the neutrophils and one for the bacteria. The control of bacterial growth requires a minimum number of neutrophils, which has been called the critical neutrophil concentration (CNC) [20–22]. Thein vitrodata in these papers have been modeled successfully with exponential growth of bacteria,B, and mass-action killing by neutrophils,N, i.e.,, revealing that the neutrophil density at which, i.e.,, defines a CNC that is independent of the bacterial density [20–22]. Conversely, Drusanoet al.[24,25] inoculated the thighs or lungs of mice with increasing doses of bacteria, and found a threshold in the initial bacterial density above which the infection grows uncontrolled (one could call this a critical bacterial concentration, CBC). A CBC is a natural outcome of a model in which the number of bacteria that can be killed by the maximum number of neutrophils approaches a maximum [12,13]. For instance, a model with a saturated killing rate,, in which a neutrophil can maximally killkhbacteria per unit of time, has a CBC of. Since this saturation model approaches the mass-action model at low bacterial densities, i.e., when, this model also has a CNC at low bacterial densities. Malkaet al.[23] vary both neutrophil and bacterial densities inex vivoexperiments, and find evidence for a double saturation model,, which also obeys mass-action kinetics at low densities, but decreases the killing rate per neutrophil as a function of the neutrophil density (for which there is additional evidence [26]). However, for reasons of simplicity, we here adhere to the conventional model in which the killing rate only saturates with the bacterial density.

2.1 Results: Mouse models

Leunget al.[26] developed a model allowing for a critical bacterial concentration, CBC, by also using the conventional saturated killing term, and in a second paper Roachet al.[10] showed that this model can explain how phage therapy can cure mice suffering from a rampant bacterial infection in the lungs by bringing the infection below the CBC. Because this model accounts for most of thein vivoobservations during successful phage therapy in mice [10], we take their model as a starting point. The model defines phage-sensitive bacteria,S, phase-resistant bacteria,R, and bacteriophages,P,

whereNis the density of neutrophils, andB=S+Ris the total density of bacteria. Bacteria grow logistically, with a net replication rate,r, have a carrying capacity,K, and a mutation rate,. Resistant bacteria may suffer from a fitness cost,. Mutations are not reversible, i.e., sensitive bacteria mutate into resistant bacteria, but mutations of resistant bacteria back to sensitive strains are not considered in the model (other mutations are included in the net replication rate). Time units are per hour and densities are defined in terms of numbers per gram of tissue [10]. At low phage densities the infection rate obeys mass action kinetics, where individual susceptible bacteria are infected at a rate. At high phage densities, i.e., when, sensitive bacteria are infected at a maximal rate, and whenthe infection rate is half-maximal. The neutrophil killing rate was also defined by a conventional saturated killing term [26], saturating as a function of the bacterial density. At low bacterial densities the killing rate obeys mass action kinetics, with individual bacteria being killed at a ratekN, i.e., proportional to the density of neutrophils. At high bacterial densities individual bacteria are killed at a rate(which is proportional to the effector:target ratio, i.e.,N/B). At very high bacterial densities, i.e., when, the killing rate by neutrophils will therefore have a negligible effect, and bacterial densities will approach their carrying capacity,K. InEq. (3)phages are cleared at aper capitarate(which is independent of their absorption to bacteria), andbdefines a burst size, i.e., the number of phages produced due to the infection of a susceptible bacterium. InEq. (4)the recruitment of neutrophils into the lungs was modeled by logistic growth at a maximum rate,, which saturated as a function of the bacterial density,B. The carrying capacity,, of this logistic equation defined the maximum density of neutrophils per gram of tissue, and their initial density (in the absence of bacteria) was defined as an independent parameter,N(0) (note that in the absence of bacteria).

We continue by mildly modifying the model and re-estimating some of the parameters. First, as we later require an explicit equation for the infected cells, we define

whereis the rate at which infected cells burst. This is a minor change as making the quasi steady state assumption (QSSA),, gives

and substituting this intoinEq. (5)we again retrieveEq. (3), with a burst sizeb.

Lytic phages typically take about half an hour to kill the bacteria they infected and produce a new generation of viral particles [27,28]. Having about two generations per hour, and a burst size of approximatelyb= 100 particles [10], an upper bound for the replication rate of phages would be that allb= 100 phages produced by an infected cell are infectious, and together rapidly infect a total of a hundred bacteria, which after half an hour each produce another generation ofb= 100 phages. If this were a synchronous process during the first few generations, we can index the number of infected cells for each half hour, i.e.,, by the map

and so on. Thus, starting with one infected cell,I0= 1, and ignoring viral clearance, there would bephages after one hour. To translate this level of expansion into the differential equations we are using, we can solve for the initial replication rate,, from, i.e., fromgiving an upper boundh–1. Note that this estimate ofis quite robust because we are taking a logarithm by solving, whereb= 100 andL= 0.5 is the latent period in hours. For instance, if the latent period were to be 1 hour, we would estimateh–1. One could also allow for infected cell clearance, by consideringI(1)=103, and obtainh–1.

For a bacterial population at carrying capacity,S=K, the initial replication rate in the Roachet al.[10] model follows fromEq. (3),

which at low phage densities would approach the unrealistically fasth–1(parameter values are provided in the caption ofFig 1). However, because the mice were treated by administering a large number of phages before the bacteria approached their carrying capacity (seeFig 1), the initialper capitareplication rate of the phages remained reasonably low, i.e., inEq. (3)and. Since we consider treatments of rampant human bacterial infections in the second half of this paper, where the bacteria have approached their carrying capacity, we limit the initial replication rate by modifying the model in two ways. First, we slow the initial replication of the phages by not making the QSSA performed inEq. (6), which limits the initial replication rate by the time scale at which bacteria are lysed (h–1[27,28]). Second, because at high bacterial densities the infection rate should no longer be limited by bacteria, we modify the infection term into a double saturation function

which naturally approacheswhen,when, and a mass-action termwhen bothSandPare small [29,30]. The model therefore changes into

whereEq. (11)has remained the same asEq. (2). The initial replication rate of this model can be computed by considering the initial phase of a small population of phages () during a rampant bacterial infection (). In the Methods (Section4.1) we derive that in mice the replication rateh–1is approached whenCFU/g (which is similar to the other saturation constant,PFU/g, inEq. (9)). Thus, if we were to use this saturation constant value forin the double saturation function,Eq. (9), we are guaranteed to have a feasible initial phage growth rate in mice.

The characteristic behavior of the Roachet al.[10] model.Panels A–C illustrate their previous results with a large bacterial infection,S(0)=7.4×107colony forming units (CFU/g) andR(0) = 0. Panel A depicts an infection without phages. Panel B depicts the same infection treated with a large dose of phages,P(2)=7.4×108plaque forming units per gram (PFU/g), given after two hours. Panel C depicts the failure of the same treatment when neutrophils cannot expand beyond their initial level. Red lines depict sensitive bacteria,S, blue lines resistant bacteria,R, purple lines total bacteria,B=S+R, green lines phages,P, and orange lines neutrophils,N.The horizontal dashed black lines depict the CBC thresholdB^Nmax=6.7×107(this CBC is absent from Panel C because neutrophils cannot expand). Panel D illustrates that in the absence of phages bacterial infections starting below the CBC grow initially, and can only be controlled after the neutrophils have expanded significantly, and that bacterial infections starting above the CBC grow uncontrolled to their logistic limit. Infections are started withS(0)=106,107toS(0) = 108CFU/g, and the line width of the curves identifies the three initial conditions. Roachet al.[10] parameters:r= 0.75 h–1,K= 1010CFU/g,ϕ=0.9,μ=2.85×10−8,β=5.4×10−8g/(PFUh),hP=1.5×107PFU/g,b= 100,dP=0.07h–1,k=8.2×10−8g/(cellh),hk=4.1×107CFU/g,rN=0.97h–1,hN=107CFU/g,N(0)=2.7×106cells g–1andNmax=2.4×107cells g–1.

2.1.1 Neutrophil recruitment.

In the original model [26], the recruitment of neutrophils into the lungs was modeled by logistic growth at a maximum rate,, which saturated as a function of the bacterial density,B; seeEq. (4). The neutrophil equation was fitted to neutrophil dynamics data obtained from the lungs of mice stimulated with LPS [14], which suggested an initial density ofneutrophils/g, and maximal density ofcells/g (i.e., an almost 9-fold increase). We modify the neutrophil model for several reasons. First, recruitment of neutrophils into tissues actually consists of two phases, i.e., a very fast phase where circulating neutrophils leave the blood locally to enter an inflamed tissue, and (2) a slower phase where the number of circulating neutrophils increases due to an increased release of immature and mature neutrophils from the bone marrow [15,16]. The 24 h time scale of the Reutershanet al.[14] experiments, and the fact that the number of neutrophils in the lungs exceeded the initial number present in the circulation, implies that the data largely reflect an increased release from the bone marrow. Second, since in Eq. (4)whenB= 0, neutrophil densities fail to return to their initial density after a bacterial infection has been cleared. We therefore propose

whereis the normal production rate of circulating neutrophils by the bone marrow, andwould allow this to be increased approximately 9-fold at high bacterial densities. The neutrophil death rateallows the neutrophil density to return to the initial state, after clearing the bacterial infection. The maximum neutrophil density that is approached for large bacterial densities is. Thus, we replace their logistic growth withEq. (14)and adopt the previous parameters by defining,and, such that we have the same initial and maximal density of neutrophils. The loss rate of the neutrophils,, defines the time scale at which the newis approached, and we have tunedandto have a similar increase of neutrophils during the first phase of the infection. Below we will adapt this model for fitting human data. The more mechanisticEq. (14)has the additional advantage that we have reasonable estimates forin humans.

2.1.2 Neutrophil killing rates.

The killing rates of the original model were based uponin vivoexperiments in the lungs of mice in which bacterial expansion was studied, as a function of the number of bacteria inoculated [24]. The data revealed that initial bacterial densities smaller thanB(0) = 107colony forming units per gram (CFU/g) can be controlled by the neutrophil response, and that larger densities grow uncontrolled [24]. By describing the data with a model with a saturated killing term similar to the one inEq. (1), i.e.,

where we combine the killing rate and the neutrophil density into, we have a model allowing for such a threshold [13]. This threshold can be found fromEq. (15)by computing the bacterial density at which. Solving the full quadratic equation defined by setting,Leunget al.[26] and Marchiet al.[31] derive the CBC threshold(the same threshold can also be found by considering bacterial densities way below the carrying capacity, i.e., by solvingforB). Since, this CBC depends on the neutrophil density, and because the model approaches two typical neutrophil densities,N(0) and, we also define two thresholds. The lower CBC threshold,, defines bacterial infections that are immediately controlled, and need no inflammatory response, and the upper CBC threshold,, defines bacterial densities above which the maximum neutrophil response fails to control the infection. The two thresholds will be depicted as dashed or dotted lines in all figures.

Drusanoet al.[24] fitted data they collected on the rate ofPseudomonas aeruginosakilling by neutrophils in a murine model of pneumonia with a mathematical model, and some of the parameters they estimated were adopted by Roachet al.[10], i.e.,r= 0.75 h–1(a doubling time of 0.9 h) andh–1. Assuming that this estimated total killing rate,, was due to a neutrophil density at its maximum, Leunget al.[26] subsequently estimated the mass action killing rate[10]. By setting the saturation constant toCFU/g, their model was given the property that bacterial infections starting at a density belowCFU/g are controlled, whereas larger densities approach the carrying capacity (seeFig 1D). For these parameters individual neutrophils can maximally killCFU/g bacteria per hour, and the upper threshold is located atCFU/g. However, the combination of their bacterial replication rate,r, their neutrophil killing rate,k, and their normal neutrophil level,, does not allow small bacterial infections to be controlled without an inflammatory response, i.e.,kN(0)/r< 1, implying that their lower threshold,, is negative, and thatN(0) is below the CNC.

To improve upon this, and to resolve issues with the identifiability of the parameters, we refit the Drusanoet al.[24] data with somewhat simpler models in the Methods (Section4.2). Since we obtain new parameter estimates from the same data, with a faster mass-action killing rate,k, and a lower saturation constant,, we obtain that small bacterial infections can be controlled without increasing the neutrophil density, and maintain the feature that large bacterial infections cannot be controlled. For these new parameter estimates the upper CBC becomesCFU/g (which seems reasonable because this is just 2-fold lower than the previous estimate [10]; see Table 2 in the Methods).

2.2 Results: Mouse experiments

The most important contribution of the Roachet al.[10] paper was to demonstrate that the neutrophil immune response is essential to make phage therapy forPseudomonas aeruginosawork in mice. Mice were infected intranasally with large doses ofP. aeruginosa, and two hours later a large dose of a Pseudomonas phage was given intranasally. Time courses of bacteria and phages were followed over 72 hours. Importantly, the phage therapy was beneficial only when the mice had normal neutrophil numbers, and when these were able to respond to the infection [10]. We illustrate their results inFig 1A–1C. A large bacterial infection that starts just above the CBC approaches carrying capacity in 10 h (Fig 1A). The same infection is controlled by a high dose of phages administered after two hours, by bringing it below the CBC within ten hours (Fig 1B). The same phage therapy fails when neutrophils cannot expand (Fig 1C), resulting in a chronic infection with phage-resistant bacteria. These model predictions were elegantly confirmed by experiments, where the failure of the therapy inFig 1Cwas tested with MyD88-deficient hosts that lack immune cell activation and recruitment to the site of infection [10]. The fact that for the original parameters the model required expansion of the neutrophils to control even very small bacterial infections is illustrated inFig 1D, where we depict a series of bacterial challenges, fromS(0) = 104to 108CFU/g, that all grow initially, and all evoke an increase in the neutrophil densities.Fig 1Dalso confirms that challenges aboveS(0) = 107grow uncontrolled (despite the fact that this is below the CBC located atCFU/g; see the horizontal black line).

We have modified the Roachet al.[10] model by replacing the logistic growth term for the increase of the neutrophils by a more natural model having a source from the bone marrow (seeEq. (14)), by making the infection rate a saturation function of both the phages and the bacteria, by not making the QSSA for the infected bacteria, by increasing the killing rate,k, and by decreasing the saturation constant,, such that bacterial infections smaller thanB(0) = 106CFU/g can be cleared without requiring increased production of neutrophils in the bone marrow. Because we have tuned the parameters of the new neutrophil model to that of the previous model, we observe very similar establishment of a chronic infection with a large challenge of bacteria (Fig 2A). The new model has two CBCs,CFU/g andCFU/g, for the normal and the maximum number of neutrophils, respectively (see Table 2). The effect of the two thresholds is illustrated inFig 2B, where we depict a similar series of bacterial challenges, fromB(0) = 105to 108CFU/g, in the absence of phages. Infections starting belowCFU/g decline immediately, and infections starting withB(0) = 107CFU/g or more grow uncontrolled (similar to the original model).

Similar behavior of the new model given byEqs. (10)–(13)and(14).Panel A is likeFig 1A, an infection without phages. Panel B is likeFig 1Dwith infections starting withS(0)=105,106,107toS(0) = 108CFU/g. Infections starting belowB(0)=B^N0decline monotonically, and those aboveS(0)=9×106CFU/g grow uncontrolled (in the absence of phages). Panel C is likeFig 1B, a controlled infection due to a large dose of phages,P(2)=7.4×108PFU/g, given after two hours. Panel D depicts a controlled infection due to the same large dose of phages, in the absence of an increase in neutrophils (α=0). Panel E depicts the evolution of resistance (likeFig 1C) obtained by decreasing the infection rate 2-fold, and starting withR(0)=5000 CFU/g resistant bacteria. Panel F is like Panel C, a controlled infection, but due to a small dose of phages,P(2) = 1 PFU/g, given after two hours. The horizontal black dashed lines represent the low and the high CBC,B^N0andB^Nmax, respectively (the high CBC is absent whenα=0). Colors as inFig 1,S: red,R: blue,B: purple,P: green, andN: orange. New parameters:σ=2.25×105cells h–1,dI=2h–1,dN=0.0833h–1(i.e., 0.5 d–1),r= 0.3 h–1,k=1.5×10−7g/(cellh),α=7.9,dN=1/12h–1,hB=9.63×106CFU/g,hN=105CFU/g, andhk=2.86×106CFU/g.

We obtain a very similar successful phage therapy with a large dose of phages after two hours (Fig 2C). A novel feature of the new model is that phage therapy can also work when the bacterial densities drop below the lower CBC; seeFig 2D, where we have setto prevent expansion of the neutrophils. Increased production of neutrophils is therefore not strictly required to make phage therapy work. The new model can also account for a similar failure of phage therapy, and the evolution of resistance (Fig 2E), when we increase the initial number of resistant bacteria and decrease the infection rate. Finally, we observe in the new model that the initial dose at which the phages are administered hardly matters because they expand vigorously in the presence of a large density of susceptible bacteria (compareFig 2CwithFig 2F), despite the fact that we have bounded the maximal replication rate by not making the QSSA, and by using the new infection term ofEq. (9). Likewise, replacing the single dose of phages with a small influx starting after two hours, representing the intravenous administration used to treat patients [2,4,6], yields a time course that is very similar to that shown inFig 2F(see below).

What have we learned from modeling phage therapy in mice? First, we confirm that the dose at which phages need to be administered to successfully control a bacterial infection hardly matters [26]. If the phages are able to infect the bacteria they will expand rapidly in a few generations, even in a model where we limit the replication rate of phages toh–1. Clinically this is interesting because it is important that the phage preparations are sterile, and have endotoxin levels below FDA-approved levels [1,2], which is easier to achieve with small than with large doses. Second, by re-estimating the rate at which neutrophils kill bacteria, we now have two thresholds in the model, and successful phage therapy no longer requires an increased neutrophil production since the bacterial infection will also become controlled when bacterial densities drop below the level that can be controlled by normal (non-inflamed) neutrophil densities. This could also be important clinically because it is not known how long neutrophil responses remain elevated during chronic bacterial infections. As a next step we develop this mouse model further by applying it to the human situation where patients are treated with cocktails of phages that the bacteria can independently evolve resistance to.

2.3 Results: A general human model

Because patients are treated with cocktails of phages, and bacteria independently evolve resistance to all of them, we have to extend the model by allowing for several strains of bacteria, and define the mutation rates at which these strains evolve from one another. Because patients are infused with phages we also add a time dependent source,V(t), to the ODE for the phages. We also add a strain-specific loss of bacteria to model treatment with antibiotics, because patients are often treated with various antibiotics, as sensitive and resistant strains frequently differ in their sensitivity to antibiotics [2]. For instance, Chanet al.[5] used phages binding the membrane proteins of a multi-drug efflux system, such that bacteria developing phage-resistance by down regulating these receptors became more sensitive to antibiotics. Because the concentration of phages in human serum can decline quite rapidly [4], we estimate the loss rate of phages,, by modeling the pharmacokinetics of the phages in serum after intravenous administration, as measured by Schooleyet al.[4] in the Methods (Section4.5). Fitting simple models to that data we estimate that in humansh–1(see Fig 9). An alternative would be to allow for a loss of phages by absorption to bacteria, and/or a loss of phages due to the (slow) development of phage-specific antibody responses by the host [2]. Finally, since one of the potential mechanisms underlying the success of phage therapy is a loss of bacterial fitness due to their evolution of phage resistance, we extend the logistic bacterial growth model by including a density independent death rate,, of the bacteria, i.e.,. This allows the steady state of a bacterial strain,, to depend on its fitness. In the previous model the steady state of bacteria in the absence of phages and neutrophils was defined as, and it did not depend on their fitness.

The Roachet al.[10] model was recently extended to study treatment with two phages [31], and since we aim to model bacterial strains resistant to several ‘functionally different’ groups of phages [32], we generalize upon this by using a conventional bit string notation, where the strings ‘000’, ‘001’, ‘010’, ‘011’, ‘100’ have 0, 1, 2, 3, and 4 as their decimal interpretation, respectively. Also ordering the phage groups from right to left in the string, we define bacteriaB1with bit strings1= 00 to be susceptible to phage groupsP1andP2, bacteriaB2with the bit strings2= 01 to be resistant to group one only, bacteriaB3with the bit strings3= 10 to be resistant to group two only, and bacteriaB4with bit strings3= 11 to be resistant to both groupsP1andP2(i.e., the bacterial strain numberiis one plus the decimal interpretation of its string). A model withdifferent phage groups, each requiring a unique resistance mechanism, therefore needs bit strings ofbits definingdifferent types of bacteria, andtypes of infected cells, and asusceptibility matrixdefining all possible routes of infection.

Assuming that each phage group requires it own unique resistance mechanism, evolving resistance to them requires different ‘mutations’, i.e., deletions, insertions and point mutations, with different mutation rates. We therefore generalize the single mutation rate,, into a vector,, where. Since each bacterial strainican evolve from any other strainjwe combine these mutations into a matrix,. To keep the model simple we also only consider ‘gain of resistance’ mutations [10,26], i.e., the upper diagonal of the mutation matrix is zero. For example, a fully susceptibleB0strain will mutate toB1andB10strains, at ratesand, respectively, and into strainB11at a rate. For these four bacterial strains the mutation matrix would look like

where each column defines the rate at whichmutates into thes, which is a loss term forand a gain for a(see below). Following the original model developed for mice [10,26], mutations occur when bacteria divide. In chronically infected patients resistant bacteria are expected to be present before phages are introduced, as bacteria have evolved rapid resistance mechanisms over millions of years of co-evolution with phages [33]. The fitness,, of bacterial strainican be estimated from data (see Methods; Section4.4), or can be defined as a weighted sum (or a product) of the costs associated with each resistance mechanism, e.g.,, whereis the cost associated with mutation,, andis the bit string of bacterial straini, andjis the phage number (and the position in the string when counting from right to left). Finally, bacteria may become more sensitive to killing by neutrophils when they become resistant by modifying their outer membrane, [7]. We allow for this by giving each bacterial strain its own rate,, of being killed by neutrophils.

Since we now consider several phage groups we need to reconsider the saturation in the infection term. As different phages,, may use the same receptor molecule on the bacteria to enter and infect, we mechanistically derive a high-dimensional infection term by ordering the phages into groups using the same receptor, and deriving a saturated infection term for all phages within one group (see the cartoon (Fig 7) and the full derivation in the Methods; Section4.3). Due to its complexity we propose to lump all phages within one functional group into one “super phage” with an “average” infection rate,, and saturation constant,. This has the additional advantage that the model can be written with one phage-resistance mechanism per group, i.e., if the bacteria down regulate the receptor used by groupi, all phages in that group would be affected. Thanks to this simplification, the form of the infection term remains similar to the one used previously [10].

Summarizing, we consider a population of sensitive and resistant bacteria,, growing logistically. Bacteria,, are infected by different groups of phages,, to form infected cells,, which burst and produce more phages,. Bacteria are killed by a population of neutrophils,N, and the rate of killing is half-maximal when. Bacteria can only be infected by the phage groups they are susceptible to, which is summarized in the susceptibility matrixdefining all possible routes of infection. All of this can be formalized into the following model:

forand, and wheredefines the number of phage groups. The vectordefines the infectivity of phage strainj,defines which bacterial strain is resistant () or sensitive () to what phage, andis the number of potentially effective antibiotics a patient is treated with, whereis a constant defining the concentration of theantibiotic, anddefines the sensitivity of bacterial strainifor antibioticj.Eq. (19)for the phages is extended with a time-dependent and phage-specific intravenous infusion rate,, and a mass-action absorption to all bacteria at a ratea(as phages can bind bacteria that they cannot infect). The latter is a simplification when bacteria down-regulate the receptor molecule that specific phages use for binding and entering as a resistance mechanism.

2.3.1 Human parameters.

The value and the biological interpretation of all parameters is specified in the Methods (Table 4). The original model was parametrized on mouse data per gram of tissue, and since bacterial and phage densities per gram of tissue are probably similar in mice and men, we use most parameters of the mouse model for the human situation (except when human estimates are available). Since patients with rampant bacterial infections typically have several tissues infected with bacteria, the neutrophil equation of the model can best be interpreted as the number of neutrophils in the circulation and peripheral tissues. In the absence of inflammation we define a normal value,N(0), and we allow this to be increased 10-fold at high bacterial densities, by setting. Because we use the faster loss rate of phages that we estimated for humans,h–1(see Methods; Section4.5), the phage densities now decline to much lower levels when the bacteria become resistant or disappear due to successful therapy. To allow a population of phages to go extinct (and not grow back from unrealistic ‘atto-fox’ levels when sensitive bacteria recover [34]), we multiply the infection rate with a steep Hill function, i.e., we use a ‘truncated’ infection rate,(note that Marchiet al.[35] use a similar function to limit the diffusion of phages in a spatial model). Thus, the infection rate is halved when, it rapidly approaches zero when, and it rapidly approacheswhen. Since we are modeling one gram of tissue, and one PFU should approximately correspond to a single phage, we setPFU/g to allow phages to go extinct when there is less than one particle in a total body weight of about 10 kg of relevant tissue (settingeither 10-fold higher or lower does not affect our main results).

We estimate the infection rate of the various phages by the successfulness of the treatment. Phage therapy will be successful when it reduces the bacterial density several orders of magnitude, from its carrying capacity to breach a CBC below which neutrophils can eliminate the infection (this CBC may depend on the bacterial strain when strains experience different killing rates,). The bacteria will only decrease when the force of infection,F, exceeds the maximal net replication rate,. At high phage densities, i.e.,, the maximum infection rate per phage group approaches(summing up to a force of infectionfor all phages). Since the phages rapidly approach high densities, we define a critical infection rate above which one phage group can control the bacteria, by solving, i.e.,g per PFU per hour. To maximize the initial replication rate of phages toh–1, we now substitute the human parameters from Table 4 intoEq. (24)and estimate thatCPU/g (for).

2.3.2 Single phage group.

To explore the behavior of this new ‘human’ model we first simulate a simplified phage therapy with a single phage group, while varying the infection rate,, and the mutation rate,. We consider a bacterial infection at steady state, with a corresponding increased number of neutrophils, and start with an arbitrary small number of phages,P1(0) = 1 PFU/g, on day zero (seeFig 3). Because bacterial densities are high the phages approach their maximum replication rate ofh–1, and although their replication rate slows down when the phage densities approach, they very rapidly expand to densities exceeding 1010PFU/g (Fig 3). Because phages replicate so rapidly when bacteria are abundant, this confirms that the initial dose hardly matters.

Exploring the human model with a single phage.Sensitive bacteria are now denoted asB0and resistant bacteria asB1(instead ofSandRin the previous model). In Panels A & B, C & D, and E & F, we setβ=10−8,β=2×10−8, andβ=3×10−8g per PFU per hour, respectively, i.e., below, at, and aboveβcrit. In Panels A and C the mutation rate is set to zero, in Panels B, D and E we have kept theμ=2.85×10−8h–1that was used by Roachet al.[10], and in Panel F the mutation rate is increased 10-fold. In the absence of resistance, the phage clears the bacterial infection whenβ≥βcrit; compare Panels A andC.In the presence of resistance, an infection rate larger thanβcritneed not be sufficient, because the resistantB1bacteria may breach the upper CBC threshold (B^Nmax); see Panel D and PanelF.In Panel B the phages do not go extinct: the system slowly approaches a steady state dominated by a large population of resistant bacteria, with neutrophils at their maximum density, and the phages and the sensitive bacteria approach this steady state in a dampened oscillation (see a longer simulation inS2A Fig). In Panels D and F the phages are declining while the sensitive bacteria,B0, are slowly recovering. Here the phages do go extinct, and the sensitive bacteria slowly take over from the resistant bacteria due to their 10% higher fitness (seeS2BandS2C Fig). For better visibility the neutrophil concentration is divided byN(0) (i.e.,Nis scaled between one and ten). Novel parameter values:ϕ1=1,ϕ2=0.9,r= 0.31 h–1,dB=0.01h–1,dP=3.5h–1,hB=3.56×107CFU/g,k1=k2=0.405h–1,a= 0, andμ^=(00μ0); see Table 4 in the Methods.

Phage therapies with phages having a low infection rate (i.e.,) fail because the sensitive bacteria just approach a somewhat lower level (seeFig 3A; where there are no resistant bacteria because we set), or only have a similar small transient effect because the resistant bacteria rapidly take over, and the bacterial infection returns to a similar high steady state (Fig 3B). Since the phages are declining after most of the bacteria have become resistant, the sensitive bacteria start to recover slowly (from day 7 onward inFig 3B). Ultimately, this also allows the phages to recover, and the system slowly approaches a steady state dominated by resistant bacteria and high neutrophil levels, and about 107CFU/g sensitive bacteria maintaining a population of about 105PFU/g phages (S2A Fig). In the absence of resistance, phages with a higher infection rate ofsuffice to control the bacterial infection (Fig 3C), but in the presence of resistance this infection rate need not be sufficient (Fig 3D), because the bacteria rapidly become resistant and the phages go extinct (S2B Fig). Disturbingly, the ‘resistant’ steady state is approached faster when the phages are more infectious, because phages having a higher infection rate lead to a steeper decline of the bacteria, which reduces the logistic competition, and allows the resistant bacteria to replicate faster (compare the total bacterial load,B, inFig 3Bwith that inFig 3D).

Increasing the infection rate further allows the phages to clear the bacterial infection in the presence of resistance (Fig 3E), because more infectious phages more rapidly deplete the sensitive bacteria, which here brings the total bacterial density below the upper threshold CBC,, before the resistant bacteria have taken over. Note that resistant bacteria start declining when the total,B, breaches the upper threshold, that the phages disappear afterwards quite rapidly due to their relatively high turnover, and that the neutrophils return to their normal steady state (Fig 3Cand3E). However, a therapy with the same more infectious phages fails when the bacteria more easily become resistant (Fig 3F), because it then takes less time for the resistant bacteria to breach the upper threshold CBC.

Summarizing, whether or not phages can control a rampant bacterial infection depends both on their infection rate and on the rate by which the bacteria develop resistance, as phage therapy will fail as soon as the resistant bacteria breach the upper CBC. Although it should always be better to treat with highly infectious phages, an important new lesson learned is that resistant bacteria expand more quickly when the phages are more infectious. Having a basic understanding of the new model we next use it to simulate phage therapy with two phages.

2.3.3 Two phage groups.

Consider a phage therapy consisting of two phages (or two cocktails of different phage groups) that are given sequentially. For simplicity, we consider phage groups with identical infection rates, and bacterial strains with identical mutation rates, killing rates, and fitness costs (seeFig 4). Since a therapy with a single phage failed because of rapid development of resistance (Fig 3C), we add on a second phage group on day five (Fig 4A), or day two of the treatment (Fig 4Band4C). The former fails because the frequency of bacteria resistant to both phage groups (B11denoted by orange lines), markedly increased when the bacterial strains resistant to the first group (B1denoted by blue lines) took over. This is a natural result because fully resistant bacteria are just one mutational step away from bacteria resistant to the first cocktail, and two steps away from susceptible bacteria (seeEq. (16)). One can prevent double-resistant bacteria to breach the upper CBC by giving the second treatment earlier (e.g., at day two;Fig 4B), which makes the combined treatment successful. Obviously, this timing depends on the initial frequencies of resistant bacteria (seeFig 4Cwhere we increased the mutation rates 10-fold), and the optimal treatment would be to give both phage groups at the same time.

Simulated phage therapy with two phages administered sequentially.The second phage is either given fairly late (at day five (Panels A and D)) or early (at day two: Panels B, C, E & F). Sensitive bacteria are now denoted asB0, and the 3 types of resistant bacteria asB1,B10orB11, respectively. Both phages have the critical infection rateβ=2×10−8g/(PFU h). In the upper row the phages are given as a small single shot at day 0, 2 or 5, i.e.,P1(0) = 1 andP2(2) = 1 orP2(5) = 1 PFU/g (see the color-matched arrows at the bottom of the graph) while in the lower row the treatment is given in the form of a continuous infusion, i.e.,V1(t)=V¯at all time points, andV2(t)=0 at the start until we setV2(t)=V¯whent> 5 days (Panel D) ort> 2 days (Panels E & F); see the color-matched horizontal lines atbottom of the graph. The first five days in Panel (A), and the first two days in Panel (B) are therefore identical tothe unsuccessful therapy shown inFig 3D. Adding a secondtherapy at day five fails to prevent the evolution of resistance (Panels A & D), starting earlier makes the same therapy successful (Panels B & D), but a 10-fold higher mutation rate again leads to failure due to the evolution of resistance (Panels C & F). Importantly, there is hardly any difference in the total bacterial load,B, between single shots of the phages (Panels A–C), and a therapy comprised of continuous infusion (Panels D–F). Phages do not go extinct in Panel (C), and not when they are administered by continuous infusion (seeS3C–S3F Fig). Novel parameter values:ϕ1=1,ϕ2=ϕ3=0.9,ϕ4=0.8,V¯=103andμ1=μ2=2.85×10−8h–1in the matrix shown inEq. (16).

This resembles the ‘hit hard and early’ concept from antiretroviral therapy for HIV [18], where combinations of drugs requiring different escape mutations increases the genetic barrier to full resistance. Similarly, experiments with phage cocktails containing various combinations of functionally different phages, demonstrate that increasing the functional diversity markedly increases the efficacy of the phage therapy [32]. A major difference between antiretroviral therapy and phage therapy is that phages are replicating, and that phage therapy therefore relies much less on the dosing schedule. This is illustrated in the bottom panels ofFig 4where we simulate a treatment consisting of a continuous infusion of the phages (rather than a small single ‘shot’ at a particular day). The high similarity in the total bacterial load,B, between the panels in the upper and lower row ofFig 4confirms that dosing makes little difference. Continuous infusion only prevents the phages from going extinct when the availability of sensitive bacteria declines (seeS3 Fig).

The most important new lesson learned is that phage therapies should be given as early as possible to slow down the evolution of full resistance. Continuous infusion may reduce the reversal of resistance in a bacterial population, but even this would not require high doses as phages start replicating as soon as susceptible bacteria reappear.

2.4 Results: The Patterson case

The first spectacular case of a patient recovering from a rampant bacterial infection with multidrug resistant (MDR)Acinetobacter baumanniiafter the administration of various cocktails of phages was described as the ‘Patterson case’ by Schooleyet al.[4] and by Strathdeeet al.[17]. Phage therapy was initiated about a 100 days after the initial infection while the patient was in coma, and although the bacteria developed resistance to all phages in the first two cocktails on a time scale of days, the patient woke up and started recovering [4,7,17]. Our aim here is to model this recovery during the first ten days of treatment.

The first two days of treatment used a cocktail of phages calledPC, that was given via the catheters draining the intra-abdominal abscess cavities. It remained unclear whether or not the phages were able to replicate in this fairly acidic environment, and the patient did not improve. A second cocktail consisting of four other phages (IV) became available on day three, and because the two days ofPC treatment seemed to have little effect, the treatment was switched to dosing intravenously. After these four days of phage therapy the patient awoke. The therapy was stopped the next day because of clinical problems, and as soon as this was diagnosed to be due to an unrelated bacterial infection, the treatment was reinstalled two days later, and continued for more than a 100 days [4]. After a week, minocycline was added to the treatment, but because the bacteria later turned out to be resistant to minocycline [7] we do not include any antibiotics in our model (i.e., we setc= 0 inEq. (17)). A third cocktail was started on day 113, and because we only model the first 10 days during which the patient recovered, this is also not considered. Both the intra-abdominal and the intravenous infusion is accounted for with theVparameter inEq. (19), and this particular time schedule of the treatments is implemented by settingon day 1 and 2 for the phages in thePC cocktail, and on days 5 and 6 for both cocktails.

The genomes of the phages were sequenced and this revealed that they were all closely related to T4-like myophages, and that three out of the four phages in thePC cocktail were actually identical [7]. Thus, the first cocktail consisted of two related phages called Maestro and AC4, while the secondIV cocktail consisted of four related phages called Navy1, Navy4, Navy71 and Navy97 (seeTable 1). T4-like myophages use long tail fibers to bind to the phage’s receptor on the cell surface, and tail fiber analysis showed the myophages used in the cocktails had two different types of tail fibers, with Maestro, AC4, and Navy71 belonging to one group, and Navy1, Navy4, and Navy97 belonging to the other cluster. This agreed well with the observation that at day two, the bacteria became resistant to all three phages in the first cluster, and remained sensitive to the three in the second cluster (Table 2in Liuet al.[7]). Importantly, this suggests that although the phages in the firstPC cocktail had little clinical effect, they were able to replicate and select for resistance. On day six, i.e., after four days of treatment, bacteria were resistant to all phages [7].

Table: Composition of the two cocktails used by Schooleyet al.[4]. The 6 phages can be classified in two resistance groups by their tail fibers [7], which is indicated in the rows.

Table: The killing parameters of Roachet al.[10], and new estimates obtained by refitting the Drusanoet al.[24] experiments in mice usingEq. (15).

The fact that resistance comes in two clusters suggests that this can be modeled with two phage groups. Since Navy71, belonging to the first resistance cluster, was only present in the second cocktail (seeTable 1), which was administered after the bacteria were already resistant to it, we did not model it separately. Because resistance to Navy1, Navy4 and Navy97 evolved simultaneously [7], we also treat the second cocktail as one uniform group of phages. Since bacteria became resistant to the phages in the second cluster by a common 6 base pair deletion in genes predicted to be involved in the synthesis of their capsule [7], a ‘mutation’ which led to the loss of the bacterial capsule should also make them resistant to the phages in first cluster. Thus, we ignore theB10bacteria and consider only two ‘mutation’ rates inEq. (16). The sequencing data confirmed that the resistant strains evolved from a common ancestor [7], and in the model we let sensitive bacteria develop resistance to cluster one at a rate, and letB0andB1strains shed their capsule at a rate. All other elements in the mutation matrix ofEq. (16)are zero (i.e., theB10strain remains absent), and we will varyandto match the data.

In the first paper on the Patterson case [4], the development of phage-resistance was measured by culturing bacterial isolates from the patient in the presence and absence of the phages in the cocktails. We fit thesein vitrodata in the Methods (Section4.4) to estimate the relative growth rates of the bacterial strains resistant to the first and second cocktail. The estimated fitness of pretreatment bacteria resistant to the first cocktail was estimated to be, and their initial abundance was estimated to be high (21%) (seeTable 3). Culturing bacterial isolates taken at day six revealed that the estimated fitness of bacteria resistant to both cocktails is about(seeTable 3), with an estimated abundance of about 30% (at day six). The mutation rate leading to theB1strain resistant to the first cocktail was set toto match the high initial abundance in the respiration experiments. The mutation rate leading to resistance to both cocktails was set to a 20-fold lower value to account for the later appearance of theB11strains (that are also resistant to the first cocktail). Note that these mutation rates are orders of magnitude higher than those used above.

Table: Fitting the respiration data of Schooleyet al.[4]. The values between the brackets denote standard errors computed from the Hessian matrix by the FME package [41]. Common parameters were similar for the pretreatment and the day-six experiments, i.e.,dB=0.042ande= 0.0036, anddB=0.051ande= 0.0029, respectively.

At the estimated abundance of about 20% resistantB1strains in the pretreatment bacterial population, the resistant strains are way above the upper CBC (see day zero inFig 5A), and hence they rapidly take over during the first two days of treatment (whatever the infection rate of the phages). Because the first cocktail was administered into the fairly acidic environment of the abdominal abscess cavities, we assigned it an infection rate below the critical infection rate, i.e.,, and nevertheless observe that this unsuccessful cocktail does select for resistance in the bacteria. After two days of treatment the bacterial population changed from 81% susceptibleB0strains, to 98%B1and 1%B11strains (see day two inFig 5A). Because this 1% abundance of fully resistantB11strains also exceed the upper CBC, treatment with the second cocktail would also be unsuccessful (whatever the infection rate of the phages); seeFig 5A. Upon the introduction of the secondIV cocktail, the dominantB1strains immediately decline, which accelerates the replication rate of the fully resistantB11strains, that rapidly become dominant. Thus, the first cocktail markedly increases the resistance to the second cocktail even before thisIV cocktail is administered. Note that increasing the infectivity of the first cocktail would only speed up the development of full resistance, without making it successful.

The Schooleyet al.[4] patient modeled with two phage groups.The first cocktail is given on day 1 (i.e., att= 0) and a second on day 3 (i.e., att= 2); see the color-matched horizontal lines at the bottom of the graph. In Panel A the infection rate of cocktails is set toβ1=βcrit/2=10−8andβ2=2βcrit=4×10−8g/(PFUh), which leads to failure of the treatment because at the estimated pretreatment level of resistance both types of resistant bacteria are above the upper CBC (withμ1=0.05andμ2=μ1/20=0.0025). In Panel B we allow for a 10-fold faster killing rate of the fully resistantB11bacteria, which clears the infection. In Panel C we decrease the second mutation rate 5-fold (toμ2=μ1/100=0.0005), which keeps the fully resistant bacteria below the upper CBC, which also makes the second cocktail successful. In Panel D the combination therapy of Panel A becomes successful when both cocktails are given at the same time, because the fully resistant bacteria have less time to expand. We use the same parameters as inFig 3(seeTable 4), but use the estimated abundances and growth rates to set the mutation ratesμ1andμ2, and the fitness valuesϕ1=1,ϕ2=0.74andϕ4=0.33(seeTable 3). For better visibility the neutrophil concentration is divided byN(0) (i.e.,Nis scaled between one and ten).

SinceFig 5Ais not in agreement with the data as the secondIV cocktail was successful, we consider two possibilities. First, the fully resistantB11strains are much more sensitive to being killed by neutrophils because they have a deficient capsule [7]. InFig 5Bwe therefore increase the killing rate ofB11bacteria,k4, by a factor 10, which is sufficient to explain the success of the treatment. Increasing the killing rate increases the upper CBC for fully resistant bacteria, which allows them to be controlled by neutrophils as soon as the total density of partially resistant bacteria has declined sufficiently by infection with theIV phages (seeFig 5B). Note that the increased killing rate also decreases the initial abundance of fully resistant bacteria (see day zero inFig 5B), which would be helpful if their abundance were to remain below the upper CBC.

A second possibility is that the initial abundances estimated from the respiration data are overestimated because (1) the relationship between the true bacterial density and the observed respiration need not be linear [36], (2) the background respiration that we subtracted from the data was large compared to the initial respiration of the bacteria, and (3) we do not have an estimate for the pretreatment abundance of the fully resistant bacteria because they hardly replicated in the respiration experiments [4]. Explaining the success of the second cocktail just requires a much lower abundance of the fully resistant bacteria, i.e., by setting, the fully resistant bacteria remain below the upper CBC, even after their rapid expansion during the first two days of treatment. This also suffices to explain the clinical success of the second cocktail (Fig 5C).

Because the resistance to the second cocktail markedly increases during the treatment with the first cocktail, we improve the unsuccessful treatment ofFig 5Aby giving both cocktails at the same time (seeFig 5D). Simultaneous administration makes the phage therapy successful because the fully resistant strain is not given the opportunity to expand during the treatment with the first cocktail. Summarizing, whether or not phage therapy can clear a rampant bacterial infection requires sufficiently low pretreatment levels of resistance, and administration of all phages at the same time may prevent the rapid expansion of fully resistant bacteria. Experiments comparing the efficacy of simultaneous and sequential treatments confirm that adding phases simultaneously tends to be the most successful [36]. Our modeling adds the disturbing insight that increasing the infection rate in an earlier cocktail leads to a more rapid development of resistance to the next one.

2.5 Results: Functional diversity of the phages

In experiments, increasing the number of functionally different phages in cocktails added to bacterial cultures increased the efficacy with which the cocktails suppressed bacterial growth in a saturating manner [32]. In our model, the major advantage of an immediate parallel treatment with several functionally different phage groups is that the initial density of the fully resistant bacteria will still be at the pretreatment level, and that for each additional phage group requiring another resistance mutation, the expected pretreatment concentration will drop by its mutation rate. The lower the pretreatment concentration the longer it will take a fully resistant strain to breach the upper CBC, and hence the more likely that the treatment works by bringing the total bacterial load below the upper CBC. To address the question whether it is better to improve the infection rate of the phages, or to increase the diversity of a cocktail of phage groups (that are all given immediately), we performed a parameter sweep increasingand(seeFig 6). Assuming identical mutation ratesfor all phage groups, we useEq. (16)to define the full mutation matrix. We vary the mutation rate in the range suggested by the cultured isolates, i.e.,and 0.01 inFig 6A–6C, respectively, and we vary the common infection rate128-fold by increasing it in 2-fold steps. This reveals that increasing the infection rate abovehardly improves the successfulness of the therapy. Note that we do not adjustwhenchanges, and that phages with an infection rate exceeding, will therefore initially replicate faster thanh–1. This does not affect our result that increasing the infection rate hardly matters, because these initial phage replication rates allow for even better control. Even extremely infectious phages require a minimum diversity of phage groups to make phage therapy successful at these pretreatment resistance levels. Thus selecting for functionally different phages, requiring different resistance mechanisms, is much more important than optimizing individual phages. This is in good agreement with the experimental results of Wrightet al.[32].

Parameter sweep, varyingβ¯and the number of phage groups.In Panel A, B and C we setμ¯=0.0025,0.005and 0.01, such that about 5%, 10% and 20% of the pretreatment bacteria are resistant to one of the phage groups, respectively. Each heat map depicts a parameter sweep, varyingβ¯128-fold, and the number of phage groups from 1 to 6 (all start att= 0). The fitness of the bacteria declines proportional to the number of resistance mutations of the bacterial strain, i.e.,ϕi=1−0.05∑j=1nPsij, wheresiis the bit string of strainiandjdenotes the position in the string. Green boxes denote successful therapy and red boxes depict failure. Wheneverβ¯>βcrit, increasing the infection rate hardly affects the successfulness of the therapy. Even at very high infection rates,β¯>100βcrit(whereβcrit=2×10−8g/(PFU h), the pretreatment resistance levels, as set by the mutation rates, require a minimum number of phage groups, i.e., 1, 2 and 3, forμ¯0=0.0025,0.005and 0.01, respectively.

3 Discussion

Having developed a novel mathematical model for phage therapy, and after confirming that this model can successfully describe clinical data, we now have anin silicotool enabling us to compare potential scenarios for treatment. Most importantly we have observed in our models that during the first days of successful treatment the replication rate of (partially) resistant bacterial strains markedly increases because they are released from the severe competition with the previously dominant pool of sensitive bacteria. Whether or not a therapy will be successful depends on time it takes resistant strains to breach the critical bacterial concentration (CBC) above which the immune system can no longer control the resistant bacteria. Unfortunately, the more infectious the phages in a first cocktail, the faster the release from competition, and the faster phage-resistant bacteria will breach this CBC. Indeed, once the infection rate of the phages suffices for controlling sensitive bacteria, a further increase of the infection rate is less important than extending the therapy with additional phages (Fig 6). It is therefore better to provide all phages at the same time [7,36], and use diverse cocktails of functionally different phages [32], to avoid this early expansion of phage-resistant bacterial strains.

These findings mimic the lessons learned from antiretroviral therapy to HIV, where one combines several drugs, each requiring different resistance mutations [19], to hit the virus hard and early [18]. The main difference with antiretroviral and antibiotic therapy is that the dosing is less important in phage therapy because phages are rapidly ‘replicating drugs’ (the maximum replication rate of phages,h–1, by far exceeds that of their bacterial hosts). This is clinically important because using lower doses would reduce the problem of residual endotoxin in phage preparations. Current recommendations of the Antibacterial Resistance Leadership Group (ARLG) in the United States are to use the highest safe and tolerated dose of bacteriophage products, although it is not clear that clinical outcome always improve with higher doses [2]. For instance, phage therapies in the United States typically use much higher doses that those in Europe, with similar clinical success rates between 55 and 70% [2]. Trials investigating the effect of dosage are ongoing [37], but the lowest doses are typically kept above 107PFU/ml. The main concern with using lower doses is that the phages may not reach the site of infection (which in our well-mixed model is never an issue). Nevertheless, because lower doses should reduce current problems with residual endotoxin in bacteriophage products, it would be interesting to further investigate this experimentally and/or clinically.

To explain the very rapid development of phage-resistance in the patient suffering from theAcinetobacterinfection, bothin vivoandin vitro[4,7], we required that a fairly large fraction (>1%) of the pretreatment bacteria were resistant to the first cocktail (seeFig 5). This is high compared to other estimates for the frequency of phage-resistant bacteria in bacterial populations [1,10]. One potential explanation for this difference would be thatAcinetobacter baumanniihas been involved in a long arms race with T4-like myophages, and has evolved multiple distinct defense systems [33]. Mechanisms allowing someAcinetobacter baumanniito rapidly shed their envelope could be an example of such an adaptation, as this confers resistance to many different phages. Not all defense systems will be present in every individual bacterium, but since several are definitely expected to be present in a large bacterial population, it seems best to chose phages that are equipped to overcome those [33].

Although we have tried to use realistic parameter values by (1) adopting parameter values estimated from phage experiments in mice [10], (2) improving the estimated killing rate of neutrophils by refitting the experiments of Drusanoet al.[24] (seeS1 Fig), (3) estimating the fitness costs and pretreatment levels of resistance by fitting time courses fromin vitroexperiments [4] (see Fig 8), (4) estimating the clearance rate of phages from human data [4] (see Fig 9), and (5) bounding the maximum replication rate of phages toh–1, several parameter estimates remain uncertain. We also tried to develop a novel model of phage therapy in humans by (1) adopting an appropriate model developed for describing phage experiments in mice [10], (2) allowing for several phages groups and many bacterial strains being sensitive, partially resistant, or fully resistant, (3) improving the neutrophil model by making it dependent on a neutrophil source from the bone marrow, but nonetheless the model remains a caricature of the complex reality of a bacterial infection distributed over various organs. The model lacks other anti-bacterial immune responses, like macrophages that can play an important role during phage therapy [38], the complement system, antibodies, and B and T cells, and lacks antibodies elicited by the phages [39]. It therefore remains unclear whether or not our description of the phage therapy of the patient in Schooleyet al.[4] is realistic. We prefer to see these simulations as examples of mechanisms that are most likely to have played a role, and we have used these results as an inspiration for further improvement of future phage therapies.

To make this model useful for interpreting future phage therapies it would probably be good to also allow for partial resistance (as the current resistance matrix is binary). Importantly, we also need better estimates of several of its crucial parameters, which requires detailed measurements during a phage therapy. For instance, knowing the bacterial resistance mechanisms for each of the cocktails [7], was essential for establishing the actual number of phage groups in our model. It would therefore be very useful to measure the time courses of the densities of bacterial strains having different resistance mechanisms, and the density of the various phages and neutrophils over time. Such longitudinal data before, during and after treatment, would allow the model to better quantify the various processes determining the successfulness of future phage therapies, and help to determine exactly how such a phage therapy manages to clear a rampant bacterial infection.

Despite these uncertainties on the parameter values and mechanistic details of our models, our major observation that the outgrowth of phage-resistant bacterial strains is fast because they are released from competition with the dominant pool of sensitive bacteria, is so general that it does not depend on the details of our model. The hitting hard and early [18,36] with cocktails setting a high genetic barrier [19,32], is a recommendation applying to any species that rapidly evolves resistance to a treatment.

4 Methods

4.1 The initial replication rate

The initial replication rate of this model can be computed by considering the initial phase of a small population of phages () during a rampant bacterial infection (), i.e., by simplifyingEq. (12)into. The initial phase of the infection is then described by the linear system,

with eigenvalues

which has only one positive root defining the initial replication rate

Solving for the new saturation constant, we find

(See Gadhamsettyet al.[40] who used a similar approach for estimating the infection rate,, during the acute phase of an HIV-1 infection). Substituting the parameters given in the captions ofFigs 1and2, and considering the upper boundh–1, we find that in mice this replication rate is approached whenCFU/g (which is similar to the other saturation constant,PFU/g, inEq. (9)).

4.2 Neutrophil killing rates

Unfortunately, the estimates for the killing rates suffer from several uncertainties. First, Drusanoet al.[24] actually fitted their data with a different model

having a saturatedper capitakilling rate, rather than a saturated total killing rate. Second, the data they used was relatively sparse and only consisted of 5 data points, i.e., 5 pairs of an initial bacterial density,B(0), and a corresponding density,B(24), after 24 h (see the circles inS1A Fig). Apparently their fitting procedure suffered from a parameter identifiability problem, because the standard deviations of most parameter estimates are larger than their mean [24]. We have refitted the data withEq. (15), while reducing the identifiability problem by (1) fixing the carrying capacity to the one estimated previously, and (2) by estimating the growth rate,r, from an experiment where Drusanoet al.[24] depleted the neutrophils and observed bacterial densities to increase fromtoCFU/g (see the bullets inS1A Fig). Since the latter is well below the carrying capacity,K= 1010CFU/g, we assume exponential growth, and solve forrfromgivingr= 0.306 h–1(i.e., a doubling time of 2.3 h). FixingrandKwe obtain the estimatesh–1andCFU/g (seeTable 2and the legend ofS1 Figfor additional information). Even the current estimates ofandare not completely reliable, because even after fixingrthe estimates ofkandare corrected negatively.

Finally, having somewhat better estimates for the total killing rate, we still need to estimate the killing rate per neutrophil,, as the neutrophil density should also be increasing in the Drusanoet al.[24] experiments (i.e.,). To estimate a new killing rate from our new estimateh–1forr= 0.3 h–1, we first apply the constraint that small bacterial infections should be controlled at normal neutrophil levels, i.e.,kN(0)>r, implying thatg/(h cell). Settingg/(h cell), which is almost 2-fold larger than the previous estimate, we next constrain the estimate offrom the fact that the lower CBC,, should be much smaller than 107CFU/g (as higher densities grow uncontrolled [24]). SolvingfromyieldsCFU/g, which translates into a maximum killing rate ofCFU/g per neutrophil per hour. For these new parameter estimates the upper CBC becomesCFU/g.

4.3 Derivation of the model

To derive the infection term mechanistically we separate the time scale of phages binding bacteria from the time scale governing the cellular dynamics. First, consider a group of phages,with, that are using the same type of receptor molecule to enter and infect the bacteria. It is important to sub-divide the phages into such groups, because there will be competition amongst phages binding the same type of receptor, and if the bacteria down-regulate this receptor as a phage-resistance mechanism, it will equally affect all phages in this particular group. Next consider one bacterium that is expressing a total ofTcopies of this particular receptor molecule. The rate at which the cell will become infected by one of the phages is determined by the number of phages bound to a receptor, and becoming internalized. UsingRfor the number of unbound receptors,for the binding (or absorption) rate of phagei, andfor the entry rate of phagei, we consider the scheme

wherewith, defines the number of receptors bound by phages of typei, anddefines the corresponding number of internalized phages (see the cartoon inFig 7A). The infection rate,F, of the bacterium is assumed to be proportional to the number of internalized phages, i.e.,.

A cartoon explaining the model ofEq. (33)(Panel A) andEq. (34)(Panel B).The heavy horizontal line depicts the bacterial cell wall. The black boxes with different indentations are receptors having different binding sites (triangular, rectangular or oval). Panel A depicts one group of phages,Pi,i=1,…,4(where the type,i, is reflected by the color), all interacting with the same triangular ligand, that can bind receptors with a triangular indentation. We depict a total ofT= 4 receptors carrying a triangular binding site. Bound receptors are indicated as a complexes,Ci, free receptors for this group asR, and internalized phases are calledPi′. The force of infection,F, is proportional to the total number of internalized phages (seeEq. (33)). Panel B extends this into three phage groups, with a total ofTi=3receptors per group. Free phagesPikonly bind receptors of groupk(see the complementary shapes), forming complexesCik, that are internalized to becomePik′phages that all contribute to the rate of infection (seeEq. (34)).

Defining the number of unbound receptors by the conservation equationwe write for the receptors bound by phages in groupi

Making the QSSAfor all thenphages binding this receptor, we arrive at

where. Summing these expressions over the phage groups yields

or

Hence

where. For the internalized phages we write, whereis their intracellular decay rate. Also making a QSSAwe finally arrive at

(For a single type of phage this can be written as the saturated infection rate ofEq. (1), i.e., as, whereand).

Sinceprovides the number of internalized phages of typei, we still need to sum over allnphage types to arrive at the infection rate imposed by all phages binding this receptor, i.e.,

whereis an infection rate of the bacterium,Fis its force of infection, andcombines the infection rate, with the binding rate, the entry rate, and the decay rate, andremains the ratio of the binding rate over the entry rate of phagej(seeFig 7A). Since different phage groups may bind different receptors (seeFig 7B), each expressed withcopies, the infection rate ofEq. (33)still needs to be summed over all different groups of phages, to obtain the overall infection rate, i.e.,

wheremis the number of phage groups,is the number of different phages in groupk, andis density of phageiin groupk, with corresponding ratesand.

Because we typically have no information about the receptor molecules the phages used in therapies bind, and sinceEq. (34)is complicated and requires two different (typically unknown) parameters for each phage, we decided to simplify the infection term by considering all phages within one group as a single phage, having an infection rate,, and saturation constant,. Summing over the different groups, each containing a single “super phage” the infection term becomes

wherennow stands for the number of groups, and we have assumed that all groups have the same saturation constant, and all receptors have the same copy number,. Thus, although the infection rate of each phage group approaches a maximum, adding additional phage groups allows the total infection rate to go up to an expected value, whereis the average infection rate.

Finally, to prevent the infection terms to become unrealistically fast when large numbers of bacteria are each expressingTreceptors, we also allow for saturation in the bacterial densities by rewritingEq. (35)into

This double-saturation function can be derived mechanically for systems where a single type of phages interact with a single type of bacteria [29], as inEq. (9), but for the high-dimensional system many types of phages and bacteria that we consider here this function remains somewhat phenomenological [30]. The main advantage of the double saturation is that the infection term approacheswhen, andwhen.

4.4 Fitting the sensitivity assays

The activity of the bacteriophage cocktails was tested by growing isolates ofA. baumanniisampled from the patient, in the presence and absence of the phagesin vitro(seeFig 3in Schooleyet al.[4]. Assuming that infected bacteria hardly replicate and that resistant bacteria are not infected, the relative growth rate observed in these experiments should reflect the relative fitness of the resistant strains. The density of the bacteria was measured over time for 20 h in “relative respirational units,” and assuming that this measure is proportional to the bacterial density, we fit a simple model for bacterial growth to these data. Finally, making the simplifying assumption that the bacterial growth rate depends linearly on the nutrient concentration, we write

whereBis the bacterial density, 0 <N< 1 the scaled nutrient concentration,eis the amount of resource required per bacterium,rthe maximum replication rate, andthe death rate. Because the observed respiration initially varied around 50 units even in the absence of bacterial growth, we subtracted the minimum of the respiration observed in each experiment, and fittedB(t) to the background-corrected data (seeFig 8). We fitthe data in each experiment simultaneously allowing only the replication rate,r, and the initial condition,B(0), to vary between the bacteria with and without all the phages (see the red and black lines inFig 8, respectively). Differences in the replication rate should reflect fitness costs, and differences in the initial condition should reflect initial densities of phage-resistant bacterial strains present within the isolates.

Fitting the sensitivity assays.Activity of the bacteriophage cocktails on bacterial isolates taken before treatment (Panel A) or after six days of treatment (PanelB). Symbols depict the data copied fromFig 3in Schooleyet al.[4] and lines depict the best fits ofEq. (37)to this data. The parameter estimates are provided inTable 3. Fitting was performed by using the Levenberg-Marquardt algorithm implemented in FME [41].

When we fit the model ofEq. (37)to the data of bacteria growing in the absence of phages we obtain an average replication rate of the various strains weighted by their frequencies (the match between the black symbols and the black lines inFig 8show that the model describes the data well). In the presence of phages (red symbols inFig 8), there is initially hardly any expansion of the bacteria, which provides support for our assumption that in these cultures sensitive bacteria are infected so rapidly that they fail to expand, and hence that the estimated replication rate should reflect the growth rate of strains that are resistant to the phages. The model also describes these data well; seeTable 3for the parameter estimates.

The first data set of Schooleyet al.[4] involved growing bacteria isolated from the patient before the onset of phage therapy with or without the phages comprising the first cocktailPC (seeFig 3Ain their paper, andFig 8Aof this paper). The model fits the data with and withoutPC phages quite well (seeFig 8A), and provides reasonable parameter estimates (Table 3), i.e., a doubling time of about an hour, 0.43 <r< 0.58 h–1, and a negligible death rate,h–1. The bacteria resistant to the full cocktail have an estimated replication rate,r= 0.43 h–1, which is about 74% of that of the average bacterial population havingr= 0.58 h–1(i.e., inTable 3we write). Since this is the relative growth with respect to theB0‘wild type’ having a fitness, we setfor theB1strain. The initial density of the resistant phages isB(0)=1.01, which is about 21% of that of the totalB(0)=4.89.

The second data set of Schooleyet al.[4] is similar except that it uses bacteria isolated from the patient at day six and the second “Navy” cocktail,IV (seeFig 3Fin their paper). These data are reproduced inFig 8Bfor the bacteria only (black bullets), and for bacteria in the presence of the fullIV cocktail (red). The model also fits these two data sets reasonably well (seeFig 8B), and provides similar parameter estimates (Table 3), i.e., a doubling time of about an hour in the absence of phagesh–1, and a negligible death rate,h–1. In the presence of phages we estimate a replication rate ofr= 0.22 h–1, or a relative growth rate ofh–1. Since this is the relative growth with respect to theB1strain having a fitness, we setfor theB11strain. At day six the density of the resistant phages isB(0)=1.95, which is about 30% of that of the totalB(0)=6.56.

4.5 Fitting the bacteriophage pharmacokinetics

The decay rate of phages in human serum was measured after the administration ofPFU of bacteriophage via intravenous injection (seeFig 4in Schooleyet al.[4]). We depict this data inFig 9by green circles, and fit them in two ways with the following model

where we track the concentration of the phages in the inoculum,, and phages in the serum,P. Phages are lost from the inoculum at rated0, a fractionfappears in the serum, where they are lost at rate. SettingandP(0)=0 we fittedP(t) as predicted byEq. (38)to the observed plasma levels, and estimatedd0= 16 h–1,h–1, and. This best fit is depicted as the red line inFig 9. We also fittedthe same data with simple exponential decay by settingI(0)=0, and estimatedandh–1, which is shown as the blue line inFig 9. Since these two estimates ofare similar, and the fit of the exponential decay model is somewhat better, we use the estimate based upon the simplest model,h–1, throughout (seeTable 4).

The phage titers measured by Schooleyet al.[4].The symbols depict the serum phage concentration data copied fromFig 4in Schooleyet al.[4], the red line the results of fitting of this data withEq. (38), and the blue line the result of fitting with an exponential decay model.

4.6 Parameter values

Table: Values, interpretations and dimensions of the parameters of the human model.

References

1. HatfullGF, DedrickRM, SchooleyRT. Phage Therapy for Antibiotic-Resistant Bacterial Infections. Annu Rev Med. 2022; 73: 197–211. doi.org/10.1146/annurev-med-080219-122208
2. PirnayJ-P, DjebaraS, SteursG, GriselainJ, CochezC, De SoirS, et al. Personalized bacteriophage therapy outcomes for 100 consecutive cases: a multicentre, multinational, retrospective observational study. Nat Microbiol. 2024; 9(6): 1434–53. doi.org/10.1038/s41564-024-01705-x
3. SchooleyRT, StrathdeeS. Treat phage like living antibiotics. Nat Microbiol. 2020; 5(3): 391–2. doi.org/10.1038/s41564-019-0666-4
4. SchooleyRT, BiswasB, GillJJ, Hernandez-MoralesA, LancasterJ, LessorL, 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. doi.org/10.1128/AAC.00954-17
5. ChanBK, TurnerPE, KimS, MojibianHR, ElefteriadesJA, NarayanD. Phage treatment of an aortic graft infected with Pseudomonas aeruginosa. Evol Med Public Health. 2018; 2018(1): 60–6. doi.org/10.1093/emph/eoy005
6. DedrickRM, Guerrero-BustamanteCA, GarlenaRA, RussellDA, FordK, HarrisK, et al. Engineered bacteriophages for treatment of a patient with a disseminated drug-resistant Mycobacterium abscessus. Nat Med. 2019; 25(5): 730–3. doi.org/10.1038/s41591-019-0437-z
7. LiuM, Hernandez-MoralesA, ClarkJ, LeT, BiswasB, Bishop-LillyKA, et al. Comparative genomics of Acinetobacter baumannii and therapeutic bacteriophages from a patient undergoing phage therapy. Nat Commun. 2022; 13(1): 3776. doi.org/10.1038/s41467-022-31455-5
8. CampbellA. Conditions for the Existence of Bacteriophage. Evolution. 1961; 15(2): 153. doi.org/10.2307/2406076
9. LevinBR, StewartFM, ChaoL. Resource-Limited Growth, Competition, and Predation: A Model and Experimental Studies with Bacteria and Bacteriophage. The American Naturalist. 1977; 111(977): 3–24. doi.org/10.1086/283134
10. RoachDR, LeungCY, HenryM, MorelloE, SinghD, Di SantoJP, 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. doi.org/10.1016/j.chom.2017.06.018
11. SchmidtC. Phage therapy’s latest makeover. Nat Biotechnol. 2019; 37(6): 581–6. doi.org/10.1038/s41587-019-0133-z
12. HollingCS. The components of predation as revealed by a study of small mammal predation of the European pine sawfly. Canadian Entomologist. 1959; 91: 293–320.
13. PilyuginSS, AntiaR. Modeling immune responses with handling time. Bull Math Biol. 2000; 62(5): 869–90. doi.org/10.1006/bulm.2000.0181
14. ReutershanJ, BasitA, GalkinaEV, LeyK. Sequential recruitment of neutrophils into lung and bronchoalveolar lavage fluid in LPS-induced acute lung injury. Am J Physiol Lung Cell Mol Physiol. 2005; 289(5): L807-15. doi.org/10.1152/ajplung.00477.2004
15. AthensJW, HaabOP, RaabSO, MauerAM, AshenbruckerH, CartwrightGE, et al. Leukokinetic studies. IV. The total blood, circulating and marginal granulocyte pools and the granulocyte turnover rate in normal subjects. J Clin Invest. 1961; 40(6): 989–95. doi.org/10.1172/JCI104338
16. TakT, WijtenP, HeeresM, PickkersP, ScholtenA, HeckAJR, et al. Human CD62Ldim neutrophils identified as a separate subset by proteome profiling and in vivo pulse-chase labeling. Blood. 2017; 129(26): 3476–85. doi.org/10.1182/blood-2016-07-727669
17. StrathdeeS, PattersonT. The Perfect Predator: A Scientist’s Race to Save Her Husband from a Deadly Superbug. New York, USA: Hachette Books. 2019.
18. HoDD. Time to hit HIV, early and hard. N Engl J Med. 1995; 333(7): 450–1. doi.org/10.1056/NEJM199508173330710
19. BeerenwinkelN, DäumerM, SingT, RahnenfuhrerJ, LengauerT, SelbigJ, et al. Estimating HIV evolutionary pathways and the genetic barrier to drug resistance. J Infect Dis. 2005; 191(11): 1953–60. doi.org/10.1086/430005
20. LiY, KarlinA, LoikeJD, SilversteinSC. A critical concentration of neutrophils is required for effective bacterial killing in suspension. Proc Natl Acad Sci U S A. 2002; 99(12): 8289–94. doi.org/10.1073/pnas.122244799
21. LiY, KarlinA, LoikeJD, SilversteinSC. Determination of the critical concentration of neutrophils required to block bacterial growth in tissues. J Exp Med. 2004; 200(5): 613–22. doi.org/10.1084/jem.20040725
22. SilversteinSC, RabadanR. How many neutrophils are enough (redux, redux)?. J Clin Invest. 2012; 122(8): 2776–9. doi.org/10.1172/JCI63939
23. MalkaR, WolachB, GavrieliR, ShochatE, Rom-KedarV. Evidence for bistable bacteria-neutrophil interaction and its clinical implications. J Clin Invest. 2012; 122(8): 3002–11. doi.org/10.1172/JCI59832
24. DrusanoGL, VanscoyB, LiuW, FikesS, BrownD, LouieA. Saturability of granulocyte kill of Pseudomonas aeruginosa in a murine model of pneumonia. Antimicrob Agents Chemother. 2011; 55(6): 2693–5. doi.org/10.1128/AAC.01687-10
25. DrusanoGL, FregeauC, LiuW, BrownDL, LouieA. Impact of burden on granulocyte clearance of bacteria in a mouse thigh infection model. Antimicrob Agents Chemother. 2010; 54(10): 4368–72. doi.org/10.1128/AAC.00133-10
26. LeungCYJ, WeitzJS. Modeling the synergistic elimination of bacteria by phage and the innate immune system. J Theor Biol. 2017; 429: 241–52. doi.org/10.1016/j.jtbi.2017.06.037
27. YoungR. Bacteriophage lysis: mechanism and regulation. Microbiol Rev. 1992; 56(3): 430–81. doi.org/10.1128/mr.56.3.430-481.1992
28. MillerES, KutterE, MosigG, ArisakaF, KunisawaT, RügerW. Bacteriophage T4 genome. Microbiol Mol Biol Rev. 2003; 67(1): 86–156. doi.org/10.1128/MMBR.67.1.86-156.2003
29. BorghansJA, de BoerRJ, SegelLA. Extending the quasi-steady state approximation by changing variables. Bull Math Biol. 1996; 58(1): 43–63. doi.org/10.1007/BF02458281
30. De BoerRJ, PerelsonAS. Towards a general function describing T cell proliferation. J Theor Biol. 1995; 175(4): 567–76. doi.org/10.1006/jtbi.1995.0165
31. MarchiJ, MinhCNN, DebarbieuxL, WeitzJS. Multi-strain phage induced clearance of bacterial infections. PLoS Comput Biol. 2025; 21(2): e1012793. doi.org/10.1371/journal.pcbi.1012793
32. WrightRCT, FrimanV-P, SmithMCM, BrockhurstMA. Functional diversity increases the efficacy of phage combinations. Microbiology (Reading). 2021; 167(12): 1110. doi.org/10.1099/mic.0.001110
33. BernheimA, SorekR. The pan-immune system of bacteria: antiviral defence as a community resource. Nat Rev Microbiol. 2020; 18(2): 113–9. doi.org/10.1038/s41579-019-0278-2
34. FowlerAC. Atto-Foxes and Other Minutiae. Bull Math Biol. 2021; 83(10): 104. doi.org/10.1007/s11538-021-00936-x
35. MarchiJ, KhalekC, GeorgeA, WeitzJS, ChaitR. Stable coexistence and transport of lytic phage infections with migrating bacterial hosts. bioRxiv. 2025. doi.org/10.1101/2025.04.18.649568
36. HallAR, De VosD, FrimanV-P, PirnayJ-P, BucklingA. 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. doi.org/10.1128/AEM.00757-12
37. TammaPD, SouliM, BillardM, CampbellJ, ConradD, EllisonDW, et al. Safety and microbiological activity of phage therapy in persons with cystic fibrosis colonized with Pseudomonas aeruginosa: study protocol for a phase 1b/2, multicenter, randomized, double-blind, placebo-controlled trial. Trials. 2022; 23(1): 1057. doi.org/10.1186/s13063-022-07047-5
38. ZborowskyS, SeuratJ, BalacheffQ, EcomardS, MuletC, MinhCNN, et al. Macrophage-induced reduction of bacteriophage density limits the efficacy of in vivo pulmonary phage therapy. Nat Commun. 2025; 16(1): 5725. doi.org/10.1038/s41467-025-61268-1
39. MiędzybrodzkiR, KasprzakH, LetkiewiczS, RogóżP, ŻaczekM, ThomasJ, et al. Pharmacokinetic and Pharmacodynamic Obstacles for Phage Therapy From the Perspective of Clinical Practice. Clin Infect Dis. 2023; 77(Suppl 5): S395–400. doi.org/10.1093/cid/ciad516
40. GadhamsettyS, CoorensT, De BoerRJ. Notwithstanding circumstantial alibis, cytotoxic T cells can be major killers of HIV-1-infected cells. J Virol. 2016; 90(16): 7066–83. doi.org/10.1128/JVI.00306-16
41. SoetaertK, PetzoldtT. Inverse Modelling, Sensitivity and Monte Carlo Analysis in R Using Package FME. Journal of Statistical Software. 2010; 33(1): 1–28. doi.org/10.18637/jss.v033.i03

---

Republished from the open web under CC0. Authors: de Boer RJ, Schooley R, Perelson AS. Read the original.