Investigating the impact of combination phage and antibiotic therapy: a modeling study

Antimicrobial resistance (AMR) is a serious threat to global health today. The spread of AMR, along with the lack of new drug classes in the antibiotic pipeline, has resulted in a renewed interest in phage therapy, which is the use of bacteriophages to treat pathogenic bacterial infections. This therapy, which was successfully used to treat a variety of infections in the early twentieth century, had been largely dismissed due to the discovery of easy to use antibiotics. However, the continuing emergence of antibiotic resistance has motivated new interest in the use of phage therapy to treat bacterial infections. Though various models have been developed to address the AMR-related issues, there are very few studies that consider the effect of phage-antibiotic combination therapy. Moreover, some of biological details such as the effect of the immune system on phage have been neglected. To address these limitations, we utilized a mathematical model to examine the role of the immune response in concert with phage-antibiotic combination therapy compounded with the effects of the immune system on the phages being used for treatment. We explore the effect of phage-antibiotic combination therapy by adjusting the phage and antibiotics dose or altering the timing. The model results show that it is important to consider the host immune system in the model and that frequency and dose of treatment are important considerations for the effectiveness of treatment. Our study can lead to development of optimal antibiotic use and further reduce the health risks of the human-animal-plant-ecosystem interface caused by AMR.

1. Introduction. Antimicrobial resistance (AMR) is a serious threat to global health. The Centers for Disease Control and Prevention (CDC) estimates that at least 2 million people become infected by antibiotic-resistant bacteria and at least 23,000 people die each year as a direct result of these infections, costing the United States $55 billion annually [3]. Infections caused by bacteria are usually treated with antibiotics, however, due to over-prescribing and mis-prescribing, many strains of bacteria have become resistant to currently available antibiotics. A list of antibiotic-resistant pathogens, a catalog of 12 families of bacteria, for which new antibiotics are urgently needed, has been provided by the World Health Organization (WHO) [21]. However, since bacteria evolve resistance to antibiotics at a relatively rapid rate, there has been less commercial interest in developing new antibiotics. Only 6 new antibiotics were approved by the Food and Drug Administration (FDA) for use in the United States from 2010 to 2016, an obvious downward trend compared to the 16 new antibiotics approved by FDA between 1983 and 1987 [41]. In 2015, a global action plan on antimicrobial resistance (GAP-AMR) was endorsed at the World Health Assembly, and one of the five strategic objectives of the GAP-AMR is to optimize the use of antimicrobial agents [65]. In 2018, the U.S. government launched the Antimicrobial Resistance Challenge to call for leaders from around the world to work together to improve antibiotic use, accelerate research on new antibiotics and antibiotic alternatives [3]. The spread of antimicrobial resistance combined with the lack of new drug classes in the antibiotic pipeline has resulted in a resurgence of interest in phage therapy.
Phage therapy is the use of bacteriophages to treat pathogenic bacterial infections. Before the widespread use of antibiotics, phage therapy was successfully applied in treating a variety of infections in the 1920s and 1930s [45]. Due to a poor understanding of the biological nature of phages, medical limitations of the day, and introduction of broader spectrum antibiotics, phage therapy was largely dismissed by most of western medicine in the 1940s [39]. However, the rise of antibiotic resistance has resulted in renewed interest in using phage to treat bacterial infections [53]. One of the first international, single-blind clinical trials of phage therapy, which aimed to target 220 burn patients with wounds infected by Escherichia coli or Pseudomonas aeruginosa, was launched in 2015 [17,33,47]. Furthermore, clinical trials are currently underway to explore phage treatment for infections caused by Staphylococcus aureus, particularly for respiratory tract infection (e.g., pneumonia), and to reduce the population of pathogens in ready-to-eat foods and meat [1,26,27,25,34,51]. In contrast to antibiotics, bacteria sensitivity to phages is largely specific for both species and strain, which can be considered as a major advantage, since the effects of antibiotics on commensal gut microbes are notorious for secondary outcomes such as antibiotic-associated diarrhea and C. difficile infection [46]. See Figure 1 for a timeline of important events in the development and use of phage therapy.
Because the problem of antibiotic resistant bacteria is complex and growing, with no known solution, various mathematical models have been proposed to explore the dynamics of the variety of systems involved. Most models focus on the transmission dynamics of antibiotic-resistant bacteria at the host population level [4, 5, 7, 8, 9, 10, 11, 14, 18, 19, 24, 28, 31, Timeline | Highlights in the development of phage as a potential therapeutic agent for bacterial infections 1915 1926 1930s 1932 1933 1937 1940s 1980s 1990s 1996 2015 [52,38,6,36,2]. Now, with increasing interest in phage therapy as an alternative or supplement to antibiotic treatment [39], mathematical models incorporating phage therapy have been developed [22,13,43,44,54,57,64,16,35,37,50,49]. In particular, Rodriguez-Gonzalez et al. [50] developed a mathematical model of phage-antibiotic combination therapy, representing the interactions among bacteria, phage, antibiotics, and the innate immune system, but ignoring the effect of immune system on phages. Some evidence shows that even though phages are not able to boost innate immunity, bacteria-boosted innate immunity acts against the phages [29], which is an important finding to further explain instances of phage ineffectiveness and to suggest better protocols for using phage therapy. To include this important component, we extended earlier models, in particular, the model developed by Rodriguez-Gonzalez, et al. [50]. The goal is to understand the role of the immune response in concert with phage-antibiotic combination therapy by introducing immune activity related to phages to the model. We aim to explore the effect of phage-antibiotic combination therapy by adjusting the phage and antibiotic doses and/or altering the timing of the dose(s). Details of the system of nonlinear, ordinary differential equations which take into account the interactions among bacteria, phage, antibiotics, and the immune system are given in Section 2. In Sections 3 and 4, equilibria and sensitivity analysis of some reduced cases are pro- vided. In Section 5, the simulation results exploring various infection and treatment scenarios are presented, followed by a discussion in Section 6.
2. Mathematical Model. We present a deterministic antibioticphage combination therapy model that describes density-dependent interactions between two strains of bacteria, phage, antibiotics, and the host immune response. The model development builds on the work by Leung and Weitz [34], and Rodriguez-Gonzalez, et al. [50]. The model in [34] includes phage-sensitive bacteria, phage, and a saturating innate immune response. The phage therapy model in [50] extended [34] to include two bacteria strains, phage therapy, antibiotic treatment and some immune response components. The model presented here adds biological functions not included in these previous models: interactions of the immune response with phage, and the decay of the immune response. See Figure 2 for a schematic diagram of our model. The phage-sensitive (antibiotic-resistant) bacteria, denoted as B P , respond to treatment by phages, P whereas the antibiotic-sensitive (phageresistant) bacteria, B A , responds to treatment by antibiotics, A. It can be assumed that bacteria is either sensitive to phages or to antibiotics due to conservation of evolutionary resources in the bacteria [15]. The total immune response (I) is activated by the presence of bacteria and phages, and attacks both bacterial strains.
The bacteria grow logistically with growth rate r i , carrying capacity K c , and density dependence B tot = i B i , where i ∈ {A, P }. We assume that phage-sensitive bacteria mutate to become antibiotic sensitive bacteria with probability µ P . Similarly, µ A represents the probability of emergence of phage-sensitive mutants from antibiotic-sensitive bacteria. Therefore the growth of the bacteria population is modeled as: where i, j ∈ {A, P } and i = j. As in [50] both populations of bacteria are killed by an activated innate immune response which includes a densitydependent immune evasion by bacteria. That is, the mass action killing

EC50
Concentration of antibiotic at which the killing rate is half its maximum 0.3697ug/ml [50] term, IB i , is scaled by the parameter 1 + Btot . See [34] for more details. The decrease in density of B A by the antibiotic treatment is approximated by a Hill function as in [50]. The phage-sensitive bacteria are infected and lysed by phage at a rate of F (P ). Following the work in [50] two phage infection modalities, F (P ), are considered -homogeneous mixing and heterogeneous mixing. The homogeneous mixing modality is given by F (P ) = φP so that the infection rate is proportional to the phage density. The second modality is given by F (P ) = φP γ where γ is the power-law exponent. The homogeneous mixing modality is assumed for our analytical results in Section 3 and sensitivity analysis in Section 4, whereas the heterogeneous mixing model is used for the numerical analysis in Section 5. Mathematical analysis is not valid with fractional exponents, but it is likely that the phage distribution is heterogeneous.
We assume that once the antibiotic treatment is administered it is injected at a constant rate where A * = A I /θ.
The growth in the phage density is due to the release of phage through lytic infection of B P at a rate of β. Free phage particles decay at a rate ω. One of the novel biological features included in this model is the effect the immune response on the phage virus. The differences in the effectiveness of phage therapy between in vitro and in vivo suggest that the infected mammalian host's immune response may be responsible for bacterial phage resistance [29]. The per capita kill rate of phage by the immune response is denoted by κ.
As in [34] we assume there is a saturated innate immune response that is activated by bacteria. We have included a saturated innate immune response that is activated by the presence of phage where β is the maximum growth rage, K I is the maximum capacity, and K M is the phage concentration at which the immune response growth rate is half its maximum. In addition, we assume that d is the rate of decay in the immune response. Parameter values are given in Table 1. 3. Analytical results. Here, we analytically explore possible treatment outcomes via equilibria analysis.
In Table 2, we summarize possible equilibria of the system, suggesting that infection dynamics can be resulted in any of the following cases under combination treatment: I. Combination treatment fails. II. Partial success is gained, since antibiotic sensitive bacteria die out as a result of combination (or only drug treatment). Yet, the equilibria analysis, detailed below, suggests that there might be up to three outcomes, indicating that the system might have bistable dynamics; i.e., the treatment outcomes might depend on the initial bacteria density, treatment doses and timing (see numerical results section). III. Successful treatment. Both phage and antibiotic-sensitive bacteria get cleared. IV. Phage treatment completely fails. It decays before clearing the phage-sensitive bacteria. V. Phage treatment fails, yet drug treatment successfully eradicates the antibiotic sensitive bacteria. VI. Drug treatment fails, yet phage therapy eradicates phage sensitive bacteria. A rigorous mathematical analysis and feasibility of these outcomes require stability analysis. We provide the detailed analysis of some of the cases. Below, we provide the details from the analysis of Case (II). Cases (IV) and (V) are detailed in Appendix. Due to complexity of the system, we explore the possible outcomes derived here using numerical experiments in Section 5. Case II. Antibiotic Sensitive Bacteria (ASB)-free equilibrium (E ASBf ). In the absence of antibiotic-sensitive (phage-resistant) bacteria, we obtain the following subsystem:

2)
where F (P ) = φP and B A = 0. Equilibria of the system are the timeindependent solutions. Here we are interested in phage treatment only, i.e., coexistence equilibrium By setting the left hand of the system equal to zero, from the first equation, we obtain, By the second equation, we also have In addition, by the third equation, we get Rearranging the equality (3.7), we have By the equality (3.8), we also have Substituting (3.9) into (3.5), we get 10) where (B † P := f 4 (I † )). Finally, by substituting the the right hand side of the equation (3.6) into (3.10), we get the following equality as a function of immune equilib-rium component, I † : .
and two asymptotes I 1,2 : Under distinct cases with respect to sign and order of the critical points I 0,1,2 , the subsystem (3.1) might have zero or up to three possible positive equilibria. Note that whenever I † > 0, we have B † P > 0 by the equation (3.6). Therefore we are looking for immune equilibrium component, I † : I † > 0 ⇒ P † > 0. This result indicates that the system might have bistable dynamics; i.e., the treatment outcomes might depend on the initial bacteria/phage density, treatment doses and timing (see Section 5).

Sensitivity Analysis.
Building on the work in [50], we adopt many parameter values from literature estimates and behavior fitting. However, several of our parameter values are not experimentally measurable. To determine the relative effect of fluctuations in parameter values on the model output, we use Matlab and Simbiology to implement the model and run a sensitivity analysis (similar to the process described in [55]). Sensitivity analysis of parameters for our model will inform us about changes to which parameters would have the most affect on the model transients. The following general steps were performed to produce a global sensitivity analysis for all the parameters over the simulation time period.
First, we established a set of reasonable parameters. The model needs to start at an admissible point in parameter space. We then used this fitted model to generate the discretized sensitivity matrix S. We then used S to rank parameters by sensitivity and set a threshold such that parameters with sensitivity below the threshold (insensitive) are fixed and parameters with sensitivity above the threshold (sensitive) are explored.
To apply this process to our model, we used the referenced values as a starting point as listed in Table 2. Most of these parameter values were used in [50] and we estimated the new parameter values for the full model to achieved biologically reasonable transient output for the model. All four observable model outputs (B A , B P , P , and I) were sampled at 10 time points (16 hours, 20 hours, 24 hours, 40 hours, 48 hours, and days 3-7). Given that there are 18 model parameters explored, a 40 × 18 discretized sensitivity matrix S is produced.
Next, we ranked the impact of each parameter on all four observable model outputs (B A , B P , P , and I) by calculating a root mean square sensitivity measure, as defined in Brun et al. [12]. For each column j of the normalized sensitivity matrix, we get Parameter j is deemed insensitive if RM S j is less than 5% of the value of the maximum RM S value calculated over all parameters. By this measure, 12 parameters were deemed insensitive, as shown in Figure 3, and fixed at their nominal values in later investigations. The strongly sensitive parameters areβ and φ. These modulate the rate of phage replication in the phage equation and also the burst rate of phage infected bacteria. Sinceβ only appears in the P equation and it actually multiplies φ, we explored the effect of φ on the model outcome.
Although not deemed as sensitive, we also chose to explore the effect of κ, the rate of removal of phage by immune cells, on the model outcomes since it is a new parameter in our extended system.  Fig. 3: Relative sensitivities. Values below 5% of the maximum sensitivity value (indicated with the dashed line) are considered insensitive.
In the numerical results below, we explore the changes to model transients that result from different choices of these sensitive parameters.

Numerical Results.
In this section, we explore numerically computed transients for some biological relevant cases of the system and apply our proposed model to investigate the interactions between bacteria, phages, antibiotics and the immune system.
Exploring the immune response. Without any treatments (either phage or antibiotic), Fig. 4 shows that the activated immune response can clear bacteria when bacterial densities (cell densities) are low enough; however, when bacterial densities are sufficiently high, the immune system cannot mount a sufficient response to clear the infection. The complicated role of the immune response in therapeutic application of phage and antibiotics are still overgeneralize here and will be further expanded upon in later version of the model.
In this model, we have included terms to allow immune cells act against the phages during the treatment [29]. This is a relevant inclusion to the model because it helps explain instances of phage ineffectiveness. One component of the new model terms, is the parameter κ, which describes the clearance of phages by the host immune system. Because this is a new addition to the model, we investigate the effect the value of κ has on the effectiveness of the combination phage-antibiotic therapy.  Table 1.  Table 1.
In Fig. 5, we use three different values of κ, with other parameter values fixed in Table 1 and initial bacterial levels B P (0) = B A (0) = 1/2 * 7.4 * 10 8 CFU/g. Also during our experiments, both phage and antibiotic therapy are received after 2 hours of the infection (P = 7.4 * 10 8 PFU/g, A = 0.035 ug/ml). As can be seen in Fig. 5(A)-(C), our results show that higher values of κ, the killing rate of phages by immune response, results in lower availability of phages at the equilibrium state, but that the final patient outcome is not different.  Table 1.
Effect of nonlinear phage absorption rate φ. In our simulations, we assume that phage infects and lyses B P bacteria at a rate F (P ), where the function F (P ) = φP 0.6 is used to account for heterogeneous mixing. In the above sensitivity analysis shown in Fig. 3, the system's transients are sensitive to the nonlinear phage absorption rate φ. We therefore have explored the predicted effectiveness of phage therapy, as it changes with altering φ. In Fig. 6, we have shown transients for three different choices of φ, with other parameter values fixed in Table 1 and initial bacterial levels B P (0) = B A (0) = 1/2 * 7.4 * 10 8 CFU/g. In all panels, both phage and antibiotic therapy are received 2 hours after the start of the simulation. In Fig. 6(A), where phage absorption rate is 2φ, the B A goes to near zero but B P stays high; while in Fig. 6(B), where phage absorption rate is 2.5φ, the same initial dose of phages is able to bring down the level of B P to zero, and we attain the trivial equilibrium; while in Fig. 6(C), where phage absorption rate is 5φ, the B P goes to zero and the process occurs faster than compared to (B).
Effect of time of administration of phage dose. Next, we investigate the effects of timing of phage therapy on the outcome of the infection. In all the simulations, antibiotics are given at the start of the simulation, the initial bacterial levels are B P (0) = B A (0) = 1/20 * 7.4 * 10 8 CFU/g (a relatively low level), the nonlinear phage absorption rate is 2φ, and the other parameter values are fixed as shown in Table 1. In Fig. 7(A), the phage dose is given 2 hours after infections. We can see that the antibioticsensitive bacteria, B A , decays quickly and goes to equilibrium near zero. Even though phage therapy lowers the B P bacteria, B P does not get completely removed from the system and a non-zero equilibrium is achieved  Table 1.
for both B A and B P . However, in Fig. 7(B), the phage is administered after 10 hrs after the start of infection: the density of B P is already high, which helps the phages to grow and in turn phages are able to reduce the density of B P . Then in the absence of B P , the phage level also goes to zero. These experiments indicates that the timing of phage therapy can be an important factor because phage effectiveness depends on the density of bacteria present in the system.
Varying time and quantity of phage dose in multi-dose regimen. We continue our experiments by varying the frequency and quantity of the phage therapy dose to explore possible outcomes of phage therapy. For both simulations in Fig. 8 the same initial infection level B P (0) = B A (0) = 1/2 * 7.4 * 10 8 CFU/g are used, the same antibiotic therapy is administered after 2 hours, and the parameter values are same as in Table 1. In the first experiment (Fig. 8(A)), we use only one dose of phage. The dose of phage (P = 7.4 * 10 8 PFU/g) is administered after 2 hours of infection. It is shown that the B A goes to nearly zero, but B P goes to a positive equilibrium (B P 0). That is, we do not have a successful treatment. In the second experiment, we explore two doses of phage. As in (Fig. 8(A)), the first dose (P = 7.4 * 10 8 PFU/g) is administered after 2 hours of infection. Then 10 hours after the first dose, the second dose (P = 2.4 * 10 12 PFU/g) is given. We found that if the amount of second dose of phage is high enough, then the density of B P goes to zero rapidly, and we obtain a successful treatment at the end. Otherwise, you need to  Table 1.
do more doses of phage treatment (See Fig. 9).
In Fig. 9, three experiments are shown. In all simulations, the same initial infection level B P (0) = B A (0) = 1/2 * 7.4 * 10 8 CFU/g is used, antibiotic therapy is administered after 2 hours, and parameter values are fixed as shown in Table 1. To conduct the comparison study, the first experiment ( Fig. 9(A)), is same as Fig. 8, i.e., only one dose of phages (P = 7.4 * 10 8 PFU/g) is administered after 2 hours of infection and it did not lead to a successful treatment. In the second experiment ( Fig. 9(B)), we administer two doses of phages. The first dose (P = 7.4 * 10 8 PFU/g) is administered after 2 hours of infection. Then 10 hours after the first dose, the second dose (P = 1.8 * 10 12 PFU/g, a relatively low value compared to 8(B)) is given. We found that even though the B P density decreases quickly after the second dose of phage, it eventually rebound. This indicates we still fail the treatment. In the third experiment ( Fig. 9(C)), we have three doses of phage therapy. The first two doses are administered as in the second experiment, i.e., the first dose (P = 7.4 * 10 8 PFU/g) is administered after 2 hours of infection. Then 10 hours after the first dose, the second dose (P = 1.8 * 10 12 PFU/g, a relatively low value compared to 8(B)) is given. Now, 10 hours after the second dose, we try the third dose (P = 4.5 * 10 11 < 1.8 * 10 12 PFU/g), and see that we can obtain a The first dose 7.4 * 10 8 PFU/g was administered after 2 hours of the beginning of infection, the second dose P = 1.8 * 10 12 PFU/g was administered after 10hrs of first dose, and the third dose P = 4.5 * 10 11 PFU/g was administered after 10hrs of second dose. Here, the initial bacterial level B P (0) = B A (0) = 1/2 * 7.4 * 10 8 CFU/g is used, antibiotic therapy is administered after 2 hours, and parameter values are fixed as shown in Table 1.
successful treatment. Hence, we believe that the number of doses and the size of the dose of phages have significant impacts on the clearance of the bacterial infection.
6. Discussion. In this work, we have analyzed a prior model, developed in [50], for the use of combination antibiotic and phage therapy for the treatment of a systemic bacterial infection. We extended it by including immune response to circulating phages. While phages are not "infectious" to humans, they are a foreign substance in the body and will elicit an inflammatory response. Additionally, the innate immune response of the patient will clear some of the phages either through filtration or through phagocytosis. Therefore, we would like to see if this dynamic is important to consider for predicting the effectiveness of the combination therapy.
By utilizing equilibria analysis, we find that the model proposed by Rodriguez-Gonzalez, et al. [50] can have six possible steady-state cases for model outcomes. In addition, our analysis suggests that in some cases, the system might display bistable dynamics; i.e., the treatment outcomes can depend on initial conditions, determined by dose of drug or phage cocktail, or timing of any of these treatments, or the frequency of these treatments in combinations. Therefore, we numerically explores outcomes of treatment options using phage therapy in combination with antibiotic treatment in order to gain insights of how to optimize the treatment outcomes.
We performed a sensitivity analysis to determine which parameters are likely to affect the transient behavior and the overall outcome of the system. To that end, we found that two parameters were of the most interest. The one with the most biological meaning (φ) was investigated and found to have an effect on the outcome of the system. It will be important in future modeling work to better estimate the number of phage released during lysis while in a human host.
The timing of the phage treatment was also important for determining patient outcome. Because phages replicate inside the bacteria, the level of bacterial infection at the time of treatment initiation influences the effectiveness of the phage therapy. Repeated dosing with phages is also helpful in clearing the bacterial infection. Determining dosing protocols and quantifying the related risks will be important for future studies.
This initial investigation has been fruitful for understanding some of the competing dynamics observed in antibiotic/phage combination therapy, and has opened the work up to further questions and lines of research: • Are there further interactions with the host immune system that need to be explored? (Innate/adaptive/filtering) • Can we determine an optimal treatment strategy?
• Do different bacterial infections require different parameter values or are there other considerations that need to be made? Some bacteria have "broad spectrum" response to phages and some require treatment with very specific phages. • How fast do bacteria develop or lose immunity to phages?
• What additional complications occur in immunocompromized individuals? There is hope that phage therapy will usher in a new line of treatment for difficult bacterial infections but there are much work needed to understand the complex dynamics and to devise effective, broadly implementable treatment protocols.
Acknowledgements. The work described herein was initiated during the Collaborative Workshop for Women in Mathematical Biology hosted by the Institute for Pure and Applied Mathematics at the University of California, Los Angeles in June 2019. Funding for the workshop was provided by IPAM, the Association for Women in Mathematics' NSF ADVANCE "Career Advancement for Women Through Research-Focused Networks" (NSF-HRD 1500481) and the Society for Industrial and Applied Mathematics. The authors thank the organizers of the IPAM-WBIO workshop (Rebecca Segal, Blerta Shtylla, and Suzanne Sindi) for facilitating this research.

APPENDIX
Case IV. Phage-free equilibrium (E + P f ). Setting P = 0, at the steadystate we obtain the following equation system: By the last equation in (.1), we have Then rearranging (.2), we obtain .
Also note that B + P = B + tot − B + A . Then by (.9), we obtain where the expressions of a i for i = 0, 1, 2 are given in (.10). Therefore the system has at most one positive phage-free equilibrium E + P f .
Case V. Phage & Antibiotic Sensitive Bacteria (ASB)-free equilibrium E P &ASBf = (B + P , 0, A + , 0, I + ). Setting P = B A = 0, we obtain the following equation system: At the steady state, by the second equation, we have Rearranging it, we obtain By the first equation, we also have Note that the equalities (.13) and (.14) are functions of I + , and intersection of both equations give the equilibrium I + component of the equilibria of the system, and the other component of the equilibria can be found by substituting the component I + into the equation (.13). It is clear that the system can have more than one phage & antibiotic sensitive bacteria-free equilibrium E P &ASBf = (B + P , 0, A + , 0, I + ).