Suppressing evolution through environmental switching
SSuppressing evolution through environmental switching
Bryce Morsky ∗ and Dervis Can Vural † Department of Physics, University of Notre Dame, Nieuwland Science Hall,Notre Dame, IN, USA Department of Biology, University of Pennsylvania, Carolyn Lynch Laboratories,Philadelphia, PA, USAFebruary 12, 2021
Abstract
Ecology and evolution under changing environments are important in many subfields ofbiology with implications for medicine. Here, we explore an example: the consequences offluctuating environments on the emergence of antibiotic resistance, which is an immense andgrowing problem. Typically, high doses of antibiotics are employed to eliminate the infectionquickly and minimize the time under which resistance may emerge. However, this strategymay not be optimal. Since resistance typically has a reproductive cost to the microbes, therecan be trade-offs between resistance and competitiveness depending on the environment. Herewe show conditions under which such trade-offs can be exploited to prevent the emergence ofresistance. We develop a stochastic Lotka-Volterra model of a microbial system with competingphenotypes: a wild strain susceptible to the antibiotic, and a mutant strain that is resistant. Weinvestigate the impact of various pulsed applications of antibiotics on population suppression.Leveraging the trade-offs between susceptibility and resistance when the antibiotic is appliedor is not, we show how a strategy of environmental switching can suppress the infection whileavoiding resistant mutants. We discuss limitations of the procedure depending on the microbeand pharmacodynamics and methods to ameliorate them.
Populations will face a variety of environmental fluctuations of both biotic and abiotic nature. Sincephenotypes typically have different reproductive success in differing environments, the dynamicsof these fluctuations can be crucial in determining phenotypic composition. Here, we consider theeffects of varying environments on the emergence and maintenance of antibiotic resistance.The rise of microbial resistance is a looming catastrophe, and prudential use of antimicrobialsis a fundamental means to prevent it [30]. Such strategies to limit the chance of resistance canbe made at all levels of disease dynamics, from population level protocols to individual patienttherapies. Studies of antibiotic resistance in vivo, in hospitals, and in the community at large usingsimple mathematical models can help address the pharmacodynamics, pharmacokinetics, andepidemiology of resistance [4, 9, 14, 19, 31, 36]. Such models have found use in effectively modellingreal-world experimental data [8, 35, 44, 50]. In particular, modelling has been used in identifyingdosing regimens that suppress the emergence of resistance [48, 49]. ∗ [email protected] † [email protected] a r X i v : . [ q - b i o . P E ] F e b here are several mechanisms by which bacteria can be resistant to antibiotics [42], an exampleof which is overexpression of the efflux pump [10, 47, 54], which bacteria use to expel antibiotics.Typically, such resistant mechanisms have a fitness cost, which can result in trade-offs betweenresistance and growth [7, 16, 33, 53]. Example costs of resistance include less energy for othercellular processes, and impaired motility. When the antibiotic is present, the resistant mechanismpays for its cost by providing a fitness advantage relative to the susceptible strain. However, whenthe antibiotic is not present, the resistant mechanism incurs a fitness disadvantage. Resistance,therefore, can be reversed under evolutionary forces by altering the environment [3]. As such, apulsed protocol, where the antibiotic is periodically applied so that the environment switches fromantibiotic to antibiotic-free regimes, may be able to eliminate the bacteria. However, there is a riskthat resistant mutants evolve to reduce the fitness cost of resistance rather than lose the resistancemechanism, whereby they could be competitive whether the antibiotic is present or not [39]. Further,mutations in regulatory genes can produce phenotypes of irreversible resistance [51]. These risks cancripple pulsed protocols aimed at controlling the infection while preventing resistance. Preventinga sustained presence of resistance is therefore a high priority.Here, we develop a mathematical model of pulsed protocols of antibiotic and antibiotic-freeregimes, switching rapidly from one environment to the other, to control a bacterial populationof Escherichia coli . We consider concentration-independent (i.e. time-dependent) bactericides suchas 𝛽 -lactams (e.g. penicillins and cephalosporins), which require high maintained concentrationsto be effective [2]. We assume that there is a maximum benefit to the concentration amount (dueeither to the pharmacodynamics of the bactericidal mechanism or tolerance of the patient to theantimicrobial). Therefore, we fix the dose concentration when applied, and find the proper periodsfor each regime that prevent the emergence of resistance while eliminating the infection.Previous theoretical studies have shown that pulsed protocols of antibiotics can eliminate bac-teria [1,12,13, 29]. However, these studies feature only a “persistent phenotype”, that neither growsnor dies under application of the antimicrobial agent [6, 56]. Bacteria may transition between thepersister type and wild type, depending on the environmental conditions. The number of persistersremain at low levels and act as a staging ground for the bacteria to repopulate after the antimicro-bial is removed. Pulsed protocols of antibiotics, however, can disrupt this process and lead to theelimination of the bacteria. Experimental studies have shown that pulsed protocols can be effectivein controlling such a system [46].Although pulsed protocols can eliminate non-persister and persister colonies, they have moredifficulty in eliminating colonies with resistant phenotypes that can grow when antibiotics arepresent. However, these protocols have been shown to be effective in containing an infectionboth theoretically and experimentally [5, 22]. In such cases, antimicrobials can act as ecologicaldisturbances and can be approximately as effective as a constant application of the antimicrobial incontrolling the bacterial load while also diminishing the probability of the emergence of resistance[5]. With short durations of high concentrations of drugs, the period under which resistance isselected for can be minimized.The above studies have explored pulsed protocols in different ways: controlling persisters, andcontrolling emergence of resistance. The resistant strain we consider here does grow in the presenceof antibiotics, and thus are not persisters. Our scenario is thus more similar to and an extensionof [5]. Our main contribution here is to show how leveraging competition can not only suppressthe emergence of resistance as in [5], but also reduce the overall bacterial load. Additionally, weexplore the impact of other important mechanisms on pulsed protocols including the evolvabilityof the bacteria and the lethality of the antibiotic.The key mechanism of our models in suppressing the population is competition between thetwo phenotypes. Competition can be modelled in several ways such as through the competitiveLotka-Volterra equation or resource competition models. Competition, however, can be low whenthe total size of the population is small (e.g. the population is well below the carrying capacityor there is a high amount of resource relative to the number of bacteria). Yet, if competition is2ufficiently strong even when the population is small, then the resistant type can be eliminatedduring the no-antibiotic phase of the pulse. We explore the scenarios where competition is and isnot substantial at low population levels. Our models also feature stochasticity, which we developin a stochastic kinetic framework [55]. We show that only when selection against the resistant typeis high when the antibiotic is off can pulsed protocols effectively eliminate the population. Stochastic birth-death processes are widely used in biological modelling [37], and, in particular,stochastic modelling of the Lotka-Volterra system [24]. Our stochastic model features microscopicprocesses of birth, death, competition, and mutation, as detailed in Box 1. These processes operateon two phenotypes 𝑋 and 𝑌 , which represent a wild-type strain, which is susceptible to the antibiotic,and a mutant strain, which is resistant, respectively.Consider first the dynamics of the birth and death processes without the presence of antibiotics,i.e. in the antibiotic-free regime. Reaction set 1 represents these processes, where 𝑏 and 𝑑 are thebirth and death rates, respectively, for the wild-type strain with 𝑏 > 𝑑 >
0. We assume that thedeath rate for the resistant strain is also 𝑑 . However, assuming a cost 𝑐 > 𝑏 − 𝑐 (costs applied to birth rather than death rates have also beenapplied similarly in ecological games [23]).In the presence of the antibiotic, the above processes still occur, but with an additional set ofreactions involving the antibiotic. Since we are considering a concentration-independent or time-dependent antibiotic, we will assume that the amount of antibiotic, ¯ 𝐴 , remains unchanged while weare in the antibiotic regime. At the maximum dose ¯ 𝐴 =
1, normalized. The antibiotic is bactericidal,and kills both types of bacteria. Reaction set 3 represents death from the antibiotic with rates 𝛼 and 𝛼 (cid:48) for the susceptible and resistant strains, respectively. Note that for 𝑋 to be susceptible and 𝑌 tobe resistant, we must have 𝛼 > 𝑏 − 𝑑 > 𝑐 + 𝛼 (cid:48) .The species also die from competition as represented in Reaction set 2. Let 𝛾 be the rate atwhich they come into contact and compete where one dies. We assume that the rate of death fromcompetition is proportional to the difference between the mean growth rate (not including the effectof antibiotics) for the focal type vs. its competitor. Thus, for same type competition, we set the rateof death to 1. However, the rate of death of a mutant in competition with a wild-type is 𝜅 > / 𝜅 .Box 1: Stochastic Lotka-Volterra processesBirth/death: 𝑋 𝑏 −→ 𝑋 𝑌 𝑏 − 𝑐 −−→ 𝑌 𝑋 𝑑 −→ ∅ 𝑌 𝑑 −→ ∅ (1)Competition: 2 𝑋 𝛾 −→ 𝑋 𝑌 𝛾 −→ 𝑌 𝑋𝑌 𝛾 / 𝜅 −−→ 𝑌 𝑋𝑌 𝛾𝜅 −−→ 𝑋 (2)Death via antibiotic: ¯ 𝐴𝑋 𝛼 −→ ∅ ¯ 𝐴𝑌 𝛼 (cid:48) −→ ∅ (3)Mutation: 𝑋 𝜇 ( − ¯ 𝐴 )+ 𝜇 (cid:48) ¯ 𝐴 −−−−−−−−−→ 𝑌 𝑌 𝜇 −→ 𝑋 (4)The bacteria can turn into the other type via mutation at rate 𝜇 , which occurs in both regimes(Reaction set 4). However, under the stress of the antibiotic, the susceptible type 𝑋 will mutate to 𝑌 at a higher rate 𝜇 (cid:48) > 𝜇 .The environmental switching is controlled by a choice of the on and off durations of the drug.We assume 100% bioavailability of the drug at application, e.g. intravenous application. Thus,when the antibiotic is “turned on,” its effects are immediate. Further, when it is “turned off,” thedissipation of the antibiotic — i.e. the rate at which it breaks into ineffective material, is metabolized,3tc. — is rapid, which is common for concentration-independent and time-dependent antibiotics.In the antibiotic regime, we apply the maximum effective dose. Thus, the set of pulsed protocolswe consider are sequences of durations of the regimes, where ¯ 𝐴 = ¯ 𝐴 = 𝑏 = .
35 Growth rate. 𝑐 = .
05 Cost for resistance. 𝑑 = . 𝛼 = . 𝛼 (cid:48) = . 𝛼 = .
04 Death rate of resistant via antibiotic bacteria. 𝛾 = − Bacteria contact rate. 𝜅 > 𝜇 = − Mutation rate. 𝜇 (cid:48) = 𝜇 = − Stress induced mutation rate of susceptible to resistant bacteria.
Table 1: Summary definitions of parameters and variables.
We conducted numerical simulations of the model to test the effects of various protocols. Weaverage realizations for each parameter combination. We use the Gillespie algorithm [18] fromJulia’s DifferentialEquations and Catalyst packages [43]. Table 1 lists the default parameter valuesused with rates per 15 minutes. We vary these values to explore nearby parameter space. Weassume that a new generation occurs after 1 /( 𝑏 − 𝑑 ) ≈ ( 𝑏 − 𝑐 − 𝑑 )/( 𝑏 − 𝑑 ) ≈ . = ⇒ 𝑐 = . E. Coli has an average rate of mutation pergenome per generation on the order of 10 − [40]. Further, under stress from the antibiotic, themutation rate can be even larger, up to ten times the non-stressed rate [28]. Thus, we consider 𝜇 (cid:48) = 𝜇 . We assume that resistant bacteria die from the antibiotic at 1 / th the rate susceptiblesdo [11], i.e. 𝛼 (cid:48) = . 𝛼 . The initial condition is a population of 100% susceptible bacteria, 𝑋 = .We explore a variety of competition parameters 𝜅 . We consider fixed on/off durations, where werepeat switching until the population is extinct or 𝑡 =
14 days.
We first consider the mean field (i.e. non-stochastic) behaviour of the system with average values (cid:104) 𝑋 (cid:105) and (cid:104) 𝑌 (cid:105) (for details of the derivation see Appendix A.) (cid:164)(cid:104) 𝑋 (cid:105) = ( 𝑏 − 𝑑 − 𝜇 ( − ¯ 𝐴 ) − ( 𝛼 + 𝜇 (cid:48) ) ¯ 𝐴 )(cid:104) 𝑋 (cid:105) + 𝜇 (cid:104) 𝑌 (cid:105) − 𝛾 (cid:104) 𝑋 (cid:105) − 𝛾𝜅 (cid:104) 𝑋 (cid:105)(cid:104) 𝑌 (cid:105) , (5) (cid:164)(cid:104) 𝑌 (cid:105) = ( 𝑏 − 𝑐 − 𝑑 − 𝜇 − 𝛼 (cid:48) ¯ 𝐴 )(cid:104) 𝑌 (cid:105) + ( 𝜇 ( − ¯ 𝐴 ) + 𝜇 (cid:48) ¯ 𝐴 )(cid:104) 𝑋 (cid:105) − 𝛾𝜅 (cid:104) 𝑋 (cid:105)(cid:104) 𝑌 (cid:105) − 𝛾 (cid:104) 𝑌 (cid:105) . (6)In either environment, both types will coexist at equilibrium due to mutations, though the lessadapted type will remain at low frequency. When the population size is sufficiently low (i.e.below the zero isoclines), competition will also be low, which will allow both types to increase inabundance. However, the higher the competition term 𝜅 , the smaller this region is. A large 𝜅 willcause the mutant strain to be suppressed in the antibiotic-free environment, and in the stochasticscenario, this effect increases the chance that the mutant strain be eliminated. In the remainder ofthe results, we detail the effects of switching the drug on and off, competition, and stochasticity inthe stochastic scenario.Figure 1 depicts a representative time series for a switching protocol vs. a constant application ofthe antibiotic. With a sufficient competitive disadvantage for resistance, i.e. high 𝜅 , we can effectively4 a)(b)(c)(d)Figure 1: Representative times series for resistance emerging under constant application, panel (a), andresults for a pulsed protocol (on and off for hrs each), panels (b)-(d). Blue are wild strain and red resistant.Pulsed protocols can suppress the bacterial load (b with 𝛼 = . and 𝜅 = ). However, if we increase theantibiotic kill rate to 𝛼 = . , the pulsed protocol fails (c). We may also observe oscillations (d with 𝛼 = . , 𝛾 = . × − , 𝜅 = , and 𝜇 = − ). 𝑋 = wild-type strain bacteria, and the remaining parameters, ifnot mentioned here, are from Table 1. control the bacterial load and antibiotic resistance. However, due to stochastic effects, resistancecan briefly rise as seen in Figure 1c before it is brought under control again. This phenomenonemerges from the interactions between the two types at different time scales. Oscillations between5he types emerge at a time scale longer than the period of the protocol. When the mutants are attheir peak, the wild-type strain is not, and vice verse. However, peaks of the wild-type strain arerelatively low as can be seen in panel d. Between day 10 until sometime during day 12, the mutantsare suppressed and the wild-type are at a relative high. This high, however, is not much higherthan when they are suppressed. As such, the total bacterial load is relatively low. Regardless ofthese waves of resistant bacteria, the time average of the total bacterial load is still less than whenthe mutant becomes established under constant application.Figures 1b and 1d are specific instances where the mutant becomes established under theconstant application protocol. However, the mutations and fluctuatios in abundances are stochastic,so it is also possible that we are fortunate and constant application drives the population extinctbefore resistance emerges. Thus, to better understand the effectiveness of therapies we must evaluatethe statistics of the bacterial load as a function of system parameters.Averaged over 50 realizations, we calculate the average total bacterial load over time for pulsedprotocols with various on and off durations and compare these to the bacterial load for constantapplication of the antibiotics. We plot these results in heat maps, where the colour indicates thelong-term bacterial load relative to the outcome from constant antibiotic application. Red indicatesthat the therapy is on average worse than constant application, yellow is on average equal, and blueindicates that it is on average better.We observe that pulsed protocols along a diagonal do best. One reason for this is that theswitching times explored here are much less than the time to reach carrying capacity in eitherregime. For example, even a day-period protocol will not reach carrying capacity (the expectedtime to reach carrying capacity is between one and two days). In such a case, the population canswing from predominantly one type to the other (see Appendix B for an example time series).However, this behaviour can still be beneficial, since each application of the antibiotic is anotherchance of eliminating the population, since switching environments drives the dominant type downpotentially to extinction before the other type can become established. In addition to plotting heatmaps, we plot the total bacterial load over time for constant and pulsed (2hrs on and off each) forvarious values of each parameter averaged over 50 realizations.Figure 2 shows that the higher the competition, the lower the diagonal (i.e. the best results comefrom protocols where the duration on is greater than the duration off). The increased competitionsuppresses the emergence of resistance even in the antibiotic environment, and thus the duration ofapplication can be longer. However, we note that the effectiveness of the optimal pulsed protocolsdrastically falls if the duration on relative to off is too high. We do not see this sensitivity when re-ducing the time on. We also observe an intermediate level of competition is best for pulsed protocolsrelative to constant application. We can see this effect in Figure 2g. Increasing 𝜅 decreases the meanbacterial load for the constant application as we would expect. Since, high competition betweenthe types will suppress the emergence of resistant mutants (which is true in both environments).However, increasing competition has an initially steeper effect upon pulsed protocols before it levelsoff. A sufficient amount of competition is required for pulsed protocols to work. As 𝜅 is increased,the difference between the outcomes of the two approaches decreases.Another reason for angle of the optimal diagonal of successful protocols in Figure 2, is due tothe relationships between the mean growth rates and the antibiotic kill rates. Figure 3 depicts theresults for various values of 𝛼 . The higher the antibiotic kill rate, the shorter the duration on forthe most successful protocols. Like in the case of 𝜅 , 𝛼 impacts the effectiveness of pulsed protocolsnonlinearly. Figure 3g depicts the mean results for various 𝛼 . The higher the 𝛼 , the better constantapplication does. However, this is not true for pulsed protocols. An intermediate value is best. Thisresult is due to the impact of 𝛼 on competitiveness. If 𝛼 is too high, then the wild-type is suppressedtoo much, and thus cannot be used to suppress the mutant strain through competition.To explore how robust our results are to mutation rates, we considered various values of 𝜇 and 𝜇 (cid:48) . We can see the effects of various 𝜇 in the rows of Figure 4, which show that the pulsed protocolsare more effective under a higher mutation rate. The first row depicts the case where there is a stress6 a) (b) (c)(d) (e) (f)(g)Figure 2: Heatmaps of the average bacterial load from pulsed protocols relative to that of constant appli-cation of the antibiotic for 𝜅 = , , , , , . Protocols that matched the average outcome of the constantapplication therapy are coloured in yellow. Those protocols that did worse are in red, and those that didbetter are in blue. Panel g depicts the total bacterial load over time for constant and pulsed ( hrs on andoff each) for various 𝜅 . The points are the results for individual realizations and the curves the average. induced mutation rate from wild-type to mutant ( 𝜇 (cid:48) = 𝜇 ). The second row depicts the resultswhere stress does not increase the mutation rate ( 𝜇 (cid:48) = 𝜇 ). The stress induced mutation makes theantibiotic environment more conducive to generating resistance, and thus makes it harder to controlthe emergence of resistance. The impact is a small relative effect on the pulsed protocols vs. constantapplication protocols. Figure 4g shows that the mutation rate impacts the constant application moreso than the pulsed protocol. This result matches intuition; the higher the mutation rate, the less likelya constant application can eliminate the colony before a mutant arises and becomes established. Insummary, the more evolvable the system, the better switching environments works.7 a) (b) (c)(d) (e) (f)(g)Figure 3: Heatmaps of the average bacterial load from pulsed protocols relative to that of constant applica-tion of the antibiotic for 𝛼 = . , . , . , . , . , . and 𝜅 = (panels a-f). Protocols that matched theaverage outcome of the constant application therapy are coloured in yellow. Those protocols that did worseare in red, and those that did better are in blue. Panel g depicts the total bacterial load over time for constantand pulsed ( hrs on and off each) for various 𝛼 . The points are the results for individual realizations andthe curves the average. In Figure 5 we explored the effect of varying the contact rate 𝛾 , and observed that switchingwas more effective for an intermediate 𝛾 . For high 𝛾 , the region where both types can grow issmall, which magnifies the impact of stochastic effects leading to elimination of emerging mutants.Further, same-type competitive interactions are also more intense, and thus the population isdriven to extinction more quickly. Therefore, a constant application is best. If 𝛾 is too low, however,there is not sufficient levels of cooperation for pulsed protocols to outperform constant ones. Thisobservation has broader implications to theory such as the paradox of the plankton. In a switching8 a) (b) (c)(d) (e) (f)(g)Figure 4: Heatmaps of the average bacterial load from pulsed protocols relative to that of constant appli-cation of the antibiotic for 𝜇 = − , − , − , 𝜇 (cid:48) = 𝜇 , 𝜇 , and 𝜅 = . Protocols that matched the averageoutcome of the constant application therapy are coloured in yellow. Those protocols that did worse are inred, and those that did better are in blue. Panel g depicts the total bacterial load over time for constant andpulsed ( hrs on and off each) for various 𝜇 . The points are the results for individual realizations and thecurves the average. environment, intermediate levels of interaction between types can result in more diversity. Herewhat we mean is that switching environments which do poorly either lead to extinction or theresistant strain dominating (i.e less diversity than if both types are controlled at low levels). Figure5d shows that the effectiveness of both constant application and pulsed protocols increases as 𝛾 increases (and thus the carrying capacity is lower). However, the gap between the two shrinks, andthus for a high 𝛾 system, the absolute bacterial load of the constant application is similarly to thelow level of the pulsed protocol. 9 a) (b) (c)(d)Figure 5: Heatmaps of the average bacterial load from pulsed protocols relative to that of constant applica-tion of the antibiotic for 𝛾 = . × − , − , × − and 𝜅 = . Protocols that matched the average outcomeof the constant application therapy are coloured in yellow. Those protocols that did worse are in red, andthose that did better are in blue. Panel d depicts the total bacterial load over time for constant and pulsed( hrs on and off each) for various 𝛾 . The points are the results for individual realizations and the curvesthe average. Our model furthers the contribution of mathematical modelling to the fight against resistance [38].We found a control strategy that can eliminate the infection while avoiding resistant mutants.Our study differs from [5] in that their study showed how resistance could be greatly reducedwhile keeping bacterial load steady. However, this only works when competition is high, else at lowpopulation levels, resistance can be maintained. Previous empirical and theoretical research has alsofound the importance of high competition in containing an infection [22]. Further, pulsed protocolsand competition can be effective in containing an infection even in a well-mixed population [22].Spatial effects, such as those found in biofilms, could heighten the degree of competition and thusthe effectiveness of the pulsed protocols. Since, spatial heterogeneity due to clumping could keepthe competition level between different types high even when the population is small relative to thecarrying capacity. The competition effects we consider can be interpreted crudely as arising fromsuch effects.Bactericides with significant post-antibiotic effects (PAE), such as fluoroquinolones, may hamperour control strategy. Such antimicrobials can impact the bacteria at sub-MIC levels long afterthey have been removed from the system [2], and as such can select for resistant bacteria after theantibiotic has been turned off. Additionally, sub MIC levels can lead to multidrug resistance throughradical-induced mutagensis [26]. Therefore, we must have a rapid dissipation of the antibiotic once10elow the MIC to prevent selection for the mutant (the range between the MIC and the point atwhich the susceptible strain is selected for) [20]. PAE is frequently caused by antibiotics that impairDNA functioning. Hence, 𝛽 -lactams, which inhibit cell wall production, are a good choice for ourtherapies. Antibiotics with a short half life would also be effective with our protocols.Intermediate drug concentration is a dimension we did not consider here. The effects of con-centration on the spread of the disease due to within host and between host dynamics have beenfound to be an important factor [45]. It is not always beneficial to have a high concentration due to aU-shaped probability of resistance emerging vs. drug concentration [15]. Here we only considereda specific concentration for all applications.We restrict ourselves to the case where the system cannot be well monitored. Clearly, the beststrategy is to alter the durations dependent on the state of the system, which can both maximizethe duration of the antibiotic regime and prevent the emergence of resistance. The more repeatedapplications of the treatment, the less likely it will work. However, such observations may not befeasible, especially when the bacterial load is small and heterogeneous.Future models could incorporate other biological factors. For example, the setting and source ofresistance (source-sink dynamics) are important factors in controlling antibiotic resistance [41]. Themethod by which resistance is spread is another important factor such as where plasmids conferresistance. In such a scenario, resistance can reemerge rapidly. Since, plasmids can remain in thepopulation due to horizontal transfer even when the plasmid confers a cost [32]. Although, thetransfer rate has been show to dramatically fall once the population is low [21]. Future modelscould also incorporate more of the complexity of interactions between the bacteria and the patients’natural flora [17, 52].Our work has applications to cancer [25, 27] and which families of chemotherapies will workbest. More generally, it applies to species conflict in varying environments and the impact of suchvariations on extinction and coexistence. Though our aim here has been to lower the overall bacterialload, in some simulations (see Figure 1d), alternating environments prevents one species dominatingthereby sustaining coexistence. This temporal heterogeneity in competitiveness can thus act as astabilizing mechanism that promotes diversity (as measured by the relative proportions of eachtype over time). Acknowledgements
We would like to thank Plamen Kamenov and Giuseppe Forte for assistance with earlier versionsof this project. This material is based upon work supported by the Defense Advanced ResearchProjects Agency under Contract No. HR0011-16-C-0062, and the University of Pennsylvania.
Author contributions
Both authors contributed to the conception of the study and the final manuscript. B.M. developedthe code for and analyzed the numerical simulations, and wrote the first draft.
Data availability
The code to run the numerical simulations and make figures and the data for the figures are availableat https://github.com/bmorsky/antibioticresistance . References [1] Acar, N., and Cogan, N. G. Enhanced disinfection of bacterial populations by nutrient andantibiotic challenge timing.
Mathematical biosciences 313 (2019), 12–32.112] AliAbadi, F. S., and Lees, P. Antibiotic treatment for animals: effect on bacterial population anddosage regimen optimisation.
International Journal of Antimicrobial Agents 14 , 4 (2000), 307–313.[3] Andersson, D. I., and Hughes, D. Antibiotic resistance and its cost: is it possible to reverseresistance?
Nature Reviews Microbiology 8 , 4 (2010), 260.[4] Austin, D., and Anderson, R. Studies of antibiotic resistance within the patient, hospitals andthe community using simple mathematical models.
Philosophical Transactions of the Royal Societyof London B: Biological Sciences 354 , 1384 (1999), 721–738.[5] Baker, C. M., Ferrari, M. J., and Shea, K. Beyond dose: Pulsed antibiotic treatment schedulescan maintain individual benefit while reducing resistance.
Scientific Reports 8 , 1 (2018), 5866.[6] Balaban, N. Q., Merrin, J., Chait, R., Kowalik, L., and Leibler, S. Bacterial persistence as aphenotypic switch.
Science 305 , 5690 (2004), 1622–1625.[7] Basra, P., Alsaadi, A., Bernal-Astrain, G., O’Sullivan, M. L., Hazlett, B., Clarke, L. M.,Schoenrock, A., Pitre, S., and Wong, A. Fitness tradeoffs of antibiotic resistance in extrain-testinal pathogenic escherichia coli.
Genome Biology and Evolution 10 , 2 (2018), 667–679.[8] Bhagunde, P., Singh, R., Ledesma, K. R., Chang, K.-T., Nikolaou, M., and Tam, V. H. Mod-elling biphasic killing of fluoroquinolones: guiding optimal dosing regimen design.
Journal ofAntimicrobial Chemotherapy 66 , 5 (2011), 1079–1086.[9] Bonhoeffer, S., Lipsitch, M., and Levin, B. R. Evaluating treatment protocols to preventantibiotic resistance.
Proceedings of the National Academy of Sciences 94 , 22 (1997), 12106–12111.[10] Borges-Walmsley, M. I., McKeegan, K. S., and Walmsley, A. R. Structure and function of effluxpumps that confer resistance to drugs.
Biochemical Journal 376 , 2 (2003), 313–338.[11] Coates, J., Park, B. R., Le, D., Şimşek, E., Chaudhry, W., and Kim, M. Antibiotic-inducedpopulation fluctuations and stochastic clearance of bacteria. eLife 7 (2018), e32976.[12] Cogan, N. Effects of persister formation on bacterial response to dosing.
Journal of TheoreticalBiology 238 , 3 (2006), 694–703.[13] Cogan, N., Brown, J., Darres, K., and Petty, K. Optimal control strategies for disinfectionof bacterial populations with persister and susceptible dynamics.
Antimicrobial Agents andChemotherapy 56 , 9 (2012), 4816–4826.[14] Czock, D., and Keller, F. Mechanism-based pharmacokinetic–pharmacodynamic modelingof antimicrobial drug effects.
Journal of Pharmacokinetics and Pharmacodynamics 34 , 6 (2007),727–751.[15] Day, T., and Read, A. F. Does high-dose antimicrobial chemotherapy prevent the evolution ofresistance?
PLoS Computational Biology 12 , 1 (2016), e1004689.[16] Ender, M., McCallum, N., Adhikari, R., and Berger-Bächi, B. Fitness cost of sccmec andmethicillin resistance levels in staphylococcus aureus.
Antimicrobial Agents and Chemotherapy48 , 6 (2004), 2295–2297.[17] Estrela, S., and Brown, S. P. Community interactions and spatial structure shape selection onantibiotic resistant lineages.
PLoS Computational Biology 14 , 6 (2018), e1006179.[18] Gillespie, D. T. A general method for numerically simulating the stochastic time evolution ofcoupled chemical reactions.
Journal of Computational Physics 22 , 4 (1976), 403–434.1219] Gloede, J., Scheerans, C., Derendorf, H., and Kloft, C. In vitro pharmacodynamic models todetermine the effect of antibacterial drugs.
Journal of Antimicrobial Chemotherapy 65 , 2 (2009),186–201.[20] Gullberg, E., Cao, S., Berg, O. G., Ilbäck, C., Sandegren, L., Hughes, D., and Andersson, D. I.Selection of resistant bacteria at very low antibiotic concentrations.
PLoS Pathogens 7 , 7 (2011),e1002158.[21] Händel, N., Otte, S., Jonker, M., Brul, S., and ter Kuile, B. H. Factors that affect transfer ofthe inci1 𝛽 -lactam resistance plasmid pesbl-283 between e. coli strains. PloS One 10 , 4 (2015),e0123039.[22] Hansen, E., Karslake, J., Woods, R. J., Read, A. F., and Wood, K. B. Antibiotics can be usedto contain drug-resistant bacteria by maintaining sufficiently large sensitive populations.
PLoSBiology 18 , 5 (2020), e3000713.[23] Hauert, C., Wakano, J. Y., and Doebeli, M. Ecological public goods games: cooperation andbifurcation.
Theoretical Population Biology 73 , 2 (2008), 257–263.[24] Huang, W., Hauert, C., and Traulsen, A. Stochastic game dynamics under demographicfluctuations.
Proceedings of the National Academy of Sciences 112 , 29 (2015), 9064–9069.[25] Katouli, A. A., and Komarova, N. L. The worst drug rule revisited: mathematical modelingof cyclic cancer treatments.
Bulletin of Mathematical Biology 73 , 3 (2011), 549–584.[26] Kohanski, M. A., DePristo, M. A., and Collins, J. J. Sublethal antibiotic treatment leads tomultidrug resistance via radical-induced mutagenesis.
Molecular Cell 37 , 3 (2010), 311–320.[27] Komarova, N. Stochastic modeling of drug resistance in cancer.
Journal of Theoretical Biology239 , 3 (2006), 351–366.[28] Kuban, W., Jonczyk, P., Gawel, D., Malanowska, K., Schaaper, R. M., and Fijalkowska, I. J.Role of escherichia coli dna polymerase iv in in vivo replication fidelity.
Journal of Bacteriology186 , 14 (2004), 4802–4807.[29] Kussell, E., Kishony, R., Balaban, N. Q., and Leibler, S. Bacterial persistence: a model ofsurvival in changing environments.
Genetics 169 , 4 (2005), 1807–1814.[30] Laxminarayan, R., Duse, A., Wattal, C., Zaidi, A. K., Wertheim, H. F., Sumpradit, N., Vlieghe,E., Hara, G. L., Gould, I. M., Goossens, H., et al. Antibiotic resistance—the need for globalsolutions.
The Lancet Infectious Diseases 13 , 12 (2013), 1057–1098.[31] Lipsitch, M., and Levin, B. R. The population dynamics of antimicrobial chemotherapy.
An-timicrobial Agents and Chemotherapy 41 , 2 (1997), 363–373.[32] Lopatkin, A. J., Meredith, H. R., Srimani, J. K., Pfeiffer, C., Durrett, R., and You, L. Persistenceand reversal of plasmid-mediated antibiotic resistance.
Nature Communications 8 , 1 (2017), 1689.[33] Martínez, J. L., and Baquero, F. Interactions among strategies associated with bacterial infec-tion: pathogenicity, epidemicity, and antibiotic resistance.
Clinical Microbiology Reviews 15 , 4(2002), 647–679.[34] Melnyk, A. H., Wong, A., and Kassen, R. The fitness costs of antibiotic resistance mutations.
Evolutionary Applications 8 , 3 (2015), 273–283.1335] Nielsen, E. I., Cars, O., and Friberg, L. E. Pk/pd indices of antibiotics predicted by a semi-mechanistic pkpd model–a step towards model-based dose optimization.
Antimicrobial Agentsand Chemotherapy (2011), AAC–00182.[36] Nielsen, E. I., and Friberg, L. E. Pharmacokinetic-pharmacodynamic modeling of antibacterialdrugs.
Pharmacological Reviews 65 , 3 (2013), 1053–1090.[37] Novozhilov, A. S., Karev, G. P., and Koonin, E. V. Biological applications of the theory ofbirth-and-death processes.
Briefings in Bioinformatics 7 , 1 (2006), 70–85.[38] Opatowski, L., Guillemot, D., Boelle, P.-Y., and Temime, L. Contribution of mathematicalmodeling to the fight against bacterial antibiotic resistance.
Current Opinion in Infectious Diseases24 , 3 (2011), 279–287.[39] Pacheco, J. O., Alvarez-Ortega, C., Rico, M. A., and Martínez, J. L. Metabolic compensationof fitness costs is a general outcome for antibiotic-resistant pseudomonas aeruginosa mutantsoverexpressing efflux pumps. mBio 8 , 4 (2017), e00500–17.[40] Perfeito, L., Fernandes, L., Mota, C., and Gordo, I. Adaptive mutations in bacteria: high rateand small effects.
Science 317 , 5839 (2007), 813–815.[41] Perron, G. G., Gonzalez, A., and Buckling, A. Source–sink dynamics shape the evolution ofantibiotic resistance and its pleiotropic fitness cost.
Proceedings of the Royal Society of London B:Biological Sciences 274 , 1623 (2007), 2351–2356.[42] Poole, K. Mechanisms of bacterial biocide and antibiotic resistance.
Journal of Applied Microbi-ology 92 (2002), 55S–64S.[43] Rackauckas, C., and Nie, Q. Differentialequations. jl – a performant and feature-rich ecosys-tem for solving differential equations in julia.
Journal of Open Research Software 5 , 1 (2017).DOI:http://doi.org/10.5334/jors.151.[44] Schmidt, S., Sabarinath, S. N., Barbour, A., Abbanat, D., Manitpisitkul, P., Sha, S., and Deren-dorf, H. Pharmacokinetic-pharmacodynamic modeling of the in vitro activities of oxazolidi-none antimicrobial agents against methicillin-resistant staphylococcus aureus.
AntimicrobialAgents and Chemotherapy 53 , 12 (2009), 5039–5045.[45] Scire, J., Hozé, N., and Uecker, H. Aggressive or moderate drug therapy for infectiousdiseases? trade-offs between different treatment goals at the individual and population levels.
PLoS Computational Biology 15 , 8 (2019), e1007223.[46] Sharma, B., Brown, A. V., Matluck, N. E., Hu, L. T., and Lewis, K. Borrelia burgdorferi, thecausative agent of lyme disease, forms drug-tolerant persister cells.
Antimicrobial Agents andChemotherapy 59 , 8 (2015), 4616–4624.[47] Sun, J., Deng, Z., and Yan, A. Bacterial multidrug efflux pumps: mechanisms, physiologyand pharmacological exploitations.
Biochemical and Biophysical Research Communications 453 , 2(2014), 254–267.[48] Tam, V. H., Ledesma, K. R., Vo, G., Kabbara, S., Lim, T.-P., and Nikolaou, M. Pharmacodynamicmodeling of aminoglycosides against pseudomonas aeruginosa and acinetobacter baumannii:identifying dosing regimens to suppress resistance development.
Antimicrobial Agents andChemotherapy 52 , 11 (2008), 3987–3993.[49] Tam, V. H., Louie, A., Deziel, M. R., Liu, W., Leary, R., and Drusano, G. L. Bacterial-populationresponses to drug-selective pressure: examination of garenoxacin’s effect on pseudomonasaeruginosa.
The Journal of Infectious Diseases 192 , 3 (2005), 420–428.1450] Tam, V. H., Schilling, A. N., Poole, K., and Nikolaou, M. Mathematical modelling responseof pseudomonas aeruginosa to meropenem.
Journal of Antimicrobial Chemotherapy 60 , 6 (2007),1302–1309.[51] Van Bambeke, F., Balzi, E., and Tulkens, P. M. Antibiotic efflux pumps.
Biochemical Pharmacology60 , 4 (2000), 457–470.[52] Wade, M. J., Harmand, J., Benyahia, B., Bouchez, T., Chaillou, S., Cloez, B., Godon, J.-J.,Boudjemaa, B. M., Rapaport, A., Sari, T., et al. Perspectives in mathematical modelling formicrobial ecology.
Ecological Modelling 321 (2016), 64–74.[53] Wang-Kan, X., Blair, J. M., Chirullo, B., Betts, J., La Ragione, R. M., Ivens, A., Ricci, V.,Opperman, T. J., and Piddock, L. J. Lack of acrb efflux function confers loss of virulence onsalmonella enterica serovar typhimurium. mBio 8 , 4 (2017), e00968–17.[54] Webber, M., and Piddock, L. The importance of efflux pumps in bacterial antibiotic resistance.
Journal of Antimicrobial Chemotherapy 51 , 1 (2003), 9–11.[55] Wilkinson, D. J.
Stochastic modelling for systems biology . CRC press, 2011.[56] Zhang, Y., Yew, W. W., and Barer, M. R. Targeting persisters for tuberculosis control.
Antimi-crobial Agents and Chemotherapy 56 , 5 (2012), 2223–2230.
A Derivation of the Fokker-Planck and mean field equations
Here we derive the mean field equations. This model includes a stress induced mutation rate 𝜇 (cid:48) (see Appendix for details). Let 𝑃 ( 𝑋 , 𝑌, 𝑡 ) be the probability of 𝑋 and 𝑌 numbers of each type attime 𝑡 . Then, given that we have a two dimensional Markov process in continuous time, the masterequation is 𝜕 𝑃 ( 𝑋 , 𝑌, 𝑡 ) 𝜕 𝑡 = 𝑏 ( 𝑋 − ) 𝑃 ( 𝑋 − , 𝑌, 𝑡 ) + ( 𝑏 − 𝑐 )( 𝑌 − ) 𝑃 ( 𝑋 , 𝑌 − , 𝑡 )+ (cid:104) 𝑑 + 𝛼 ¯ 𝐴 + 𝛾 ( 𝑋 + ) + 𝛾𝜅 𝑌 (cid:105) ( 𝑋 + ) 𝑃 ( 𝑋 + , 𝑌, 𝑡 )+ (cid:104) 𝑑 + 𝛼 (cid:48) ¯ 𝐴 + 𝛾𝜅 𝑋 + 𝛾 ( 𝑌 + ) (cid:105) ( 𝑌 + ) 𝑃 ( 𝑋 , 𝑌 + , 𝑡 )+ 𝜇 ( 𝑌 + ) 𝑃 ( 𝑋 − , 𝑌 + , 𝑡 ) + ( 𝜇 ( − ¯ 𝐴 ) + 𝜇 (cid:48) ¯ 𝐴 )( 𝑋 + ) 𝑃 ( 𝑋 + , 𝑌 − , 𝑡 )− (cid:104) ( 𝑏 + 𝑑 + ( − ¯ 𝐴 ) 𝜇 + ( 𝛼 + 𝜇 (cid:48) ) ¯ 𝐴 ) 𝑋 + ( 𝑏 − 𝑐 + 𝑑 + 𝜇 + 𝛼 (cid:48) ¯ 𝐴 ) 𝑌 + 𝛾 𝑋 + (cid:16) 𝜅 + 𝜅 (cid:17) 𝛾 𝑋𝑌 + 𝛾 𝑌 (cid:105) 𝑃 ( 𝑋 , 𝑌, 𝑡 ) . (7)Taking the time derivative of the average values (cid:104) 𝑋 ( 𝑡 )(cid:105) = (cid:205) 𝑋,𝑌 𝑃 ( 𝑋 , 𝑌, 𝑡 ) and (cid:104) 𝑌 ( 𝑡 )(cid:105) = (cid:205) 𝑋,𝑌 𝑃 ( 𝑋 , 𝑌, 𝑡 ) ,we can find the mean field model: (cid:164)(cid:104) 𝑋 (cid:105) = ( 𝑏 − 𝑑 − 𝜇 ( − ¯ 𝐴 ) − ( 𝛼 + 𝜇 (cid:48) ) ¯ 𝐴 )(cid:104) 𝑋 (cid:105) + 𝜇 (cid:104) 𝑌 (cid:105) − 𝛾 (cid:104) 𝑋 (cid:105) − 𝛾𝜅 (cid:104) 𝑋 (cid:105)(cid:104) 𝑌 (cid:105) , (8) (cid:164)(cid:104) 𝑌 (cid:105) = ( 𝑏 − 𝑐 − 𝑑 − 𝜇 − 𝛼 (cid:48) ¯ 𝐴 )(cid:104) 𝑌 (cid:105) + ( 𝜇 ( − ¯ 𝐴 ) + 𝜇 (cid:48) ¯ 𝐴 )(cid:104) 𝑋 (cid:105) − 𝛾𝜅 (cid:104) 𝑋 (cid:105)(cid:104) 𝑌 (cid:105) − 𝛾 (cid:104) 𝑌 (cid:105) . (9)Here we are assuming that (cid:104) 𝑋𝑋 (cid:105) = (cid:104) 𝑋 (cid:105) , (cid:104) 𝑋𝑌 (cid:105) = (cid:104) 𝑋 (cid:105)(cid:104) 𝑌 (cid:105) , and (cid:104) 𝑌 (cid:105) = (cid:104) 𝑌 (cid:105) .15 Further simulations results
Here we depicts other simulations results. In Figure 6 we depict an example time series for a oneday on and off protocol where the bacteria is not eliminated. Longer durations on or off can lead tothe population reaching the carrying capacity, after which they will be markedly less successful atsuppressing the bacteria than the constant application therapy.