Mechanisms of antibiotic action shape the fitness landscapes of resistance mutations.

Antibiotic-resistant pathogens are a major public health threat. A deeper understanding of how an antibiotic's mechanism of action influences the emergence of resistance would aid in the design of new drugs and help to preserve the effectiveness of existing ones. To this end, we developed a model that links bacterial population dynamics with antibiotic-target binding kinetics. Our approach allows us to derive mechanistic insights on drug activity from population-scale experimental data and to quantify the interplay between drug mechanism and resistance selection. We find that both bacteriostatic and bactericidal agents can be equally effective at suppressing the selection of resistant mutants, but that key determinants of resistance selection are the relationships between the number of drug-inactivated targets within a cell and the rates of cellular growth and death. We also show that heterogeneous drug-target binding within a population enables resistant bacteria to evolve fitness-improving secondary mutations even when drug doses remain above the resistant strain's minimum inhibitory concentration. Our work suggests that antibiotic doses beyond this "secondary mutation selection window" could safeguard against the emergence of high-fitness resistant strains during treatment.

Introduction

The emergence and spread of antibiotic-resistant bacterial pathogens is an urgent global problem that threatens to undermine one of the most essential components of modern medicine [1]. Antibiotic resistance is also expensive, adding an average of US $1400 to the costs of treatment for each of the 2.8 million patients who become infected with a drug-resistant bacterium in the United States annually [2], [3], [4]. The scarcity of promising new antimicrobial drugs with novel mechanisms of action further exacerbates the challenges associated with managing the spread of resistance [5], [6]. Given the increasing incidence of resistant bacterial infections and the lack of new drugs on the horizon, clinicians, researchers, and global leaders must act to preserve the effectiveness of the world’s existing antibiotic drug arsenal [1].

Antibiotic treatment induces a strong selective pressure on bacterial populations to evolve resistance [7], [8]. Resistance mutations raise the minimum inhibitory concentration (MIC) of an antibiotic, the amount of drug needed to suppress the growth of a bacterial culture [9]. However, alleles that confer drug resistance also frequently carry fitness costs [10], [11], [12], predominantly because antibiotics target vital cellular functions (such as DNA replication and protein synthesis). Resistance mechanisms reduce the ability of a drug to disrupt its target, but do so at the expense of optimal physiological function [13].

With few exceptions [14], resistance-causing alleles induce fitness impairments in both drug-free and drug-containing environments, though resistant strains may only suffer a strict competitive disadvantage (i.e. a slower growth rate) against sensitive strains in drug-free conditions. A range of antibiotic concentrations therefore exists within which drug-resistant strains have a selective advantage over their drug-susceptible counterparts. Drugs dosed within this “resistance selection window” (also called the “mutant selection window”) favor the proliferation of drug-resistant subpopulations [15], [16], [17]. Recent advances in antimicrobial pharmacodynamics have leveraged resistance selection windows to design dosing strategies that minimize the selection of resistant pathogens without sacrificing treatment efficacy [17], [18], [19].

The existence of resistance mutations that confer fitness impairments in both drug-free and drug-containing environments implies that resistant strains face selective pressures to evolve secondary mutations that alleviate these impairments, and that these selective pressures exist even under continuous drug exposure [20], [21]. Secondary mutations can increase bacterial fitness (through faster growth rates) in the absence of drugs, or they can confer elevated levels of drug tolerance to preexisting resistant subpopulations (through attenuated drug-target interactions, faster growth rates in the presence of drugs, or both). In the case of increased bacterial fitness, secondary mutations enable drug-resistant mutants to compete against drug-susceptible strains in resource-limited, antibiotic-free environments [10], [22], [23], and are implicated in the spread of drug resistance across a wide range of timescales and clinical settings [24]. In the case of increased drug tolerance, secondary mutations can be the underlying cause of treatment failure [25], [26]. Elucidating the dynamics of secondary mutation emergence during treatment is thus crucial for managing the spread of resistance.

Since resistance mutations are frequently associated with fitness costs [11], [12] both in vivo [27] and in vitro [28], studies on the resistance selection window and on secondary adaptation have yielded valuable insights into the emergence of drug-resistant bacteria during treatment. However, the design of optimal resistance-mitigating drug dosing strategies remains challenging for two reasons. One obstacle is that bacteria may acquire resistance through a multitude of mechanisms that reduce antibiotic efficacy [29]. These molecular mechanisms may themselves influence the fitness landscape of resistance mutations (that is, the relationship between the fitness cost of resistance and the selective advantage conferred by the resistance mutation in drug-containing environments) [30]. A second challenge is that an antibiotic’s mechanism of action may affect the strength of selection for resistant strains over drug-susceptible strains during treatment. One important feature of an antibiotic’s cellular-level mechanism of action is whether the drug controls bacterial populations by increasing the rate of bacterial killing (i.e. bactericidal action) or by decreasing the rate of bacterial replication (i.e. bacteriostatic action). Clinicians and researchers alike have argued that these modes of antimicrobial action influence the dynamics of resistance selection [31], [32].

The design of resistance-mitigating antibiotic usage therefore depends on an understanding of how a drug’s mechanism of action, a pathogen’s mechanism of resistance, and the fitness landscape of resistance affect selection pressures during treatment. Tractable and quantitative strategies for systematically exploring all of these factors have so far been lacking. To address this gap, we developed a dynamical model that simulates the growth and death of bacterial populations under antibiotic exposure using molecular-scale descriptions of drug-target binding kinetics and cellular-scale descriptions of a drug’s mechanism of action. In our model, higher numbers of inactivated drug-target complexes within a cell lead to increases in antibiotic effect (either bacteriostatic, bactericidal, or a combination of the two). The relationship between drug-target inactivation and antibiotic effect can take the shape of a linear (i.e. gradual) or stepwise (i.e. sudden) function, as well as other intermediate forms (Supplementary Fig. S1). The model enables us to estimate critical pharmacodynamic parameters from experimental datasets as effectively as with classical approaches [33], to simulate the fitness landscapes of resistance mutations against drugs with diverse mechanisms of action, and to quantify the probability of secondary mutation emergence within resistant subpopulations of bacteria during treatment.

The mathematical model described here is a linear case of a nonlinear formulation (COMBAT) reported previously to study the influence of drug-target binding kinetics on optimal antibiotic dosing [34]. Linearization results in a > 102-fold computational speed-up that enables us to robustly fit experimental kill-curve data and to simulate antibiotic dose–response relationships at high resolution. Our linear formulation also allows us to calculate an antibiotic’s MIC directly from experimentally measurable molecular parameters. We leverage the mathematical tractability and computational efficiency of our model to investigate the selective pressures that antibiotics with diverse mechanisms of action induce on growing bacterial populations, a task that would be impractical with previous approaches.

We find that bacteria with resistance mechanisms that confer even modest reductions in drug-target binding affinity can incur strikingly high (80–99 %) fitness costs while still maintaining higher drug tolerances than their susceptible counterparts, regardless of the antibiotic’s mechanism of action. We also find that drugs with stepwise effects on bacterial growth and death as a function of target inactivation have narrower resistance selection windows than do drugs with linear effects. However, our model suggests that whether a drug acts primarily through bactericidal or bacteriostatic action has comparatively little influence on the strength of resistance selection during treatment. We further demonstrate that, even with aggressive treatment regimens, heterogeneous drug-target occupancy within a population enables fitness-impaired resistant strains to develop secondary mutations that can lead to treatment failure. Our work cautions that fitness costs may not limit the emergence of resistant strains that evolve through reductions in drug-target binding affinity. We propose the “secondary mutant selection window” as a novel pharmacodynamic characteristic of a drug that should be assessed alongside other classic parameters such as the MIC and the resistance selection window when designing robust resistance-mitigating antibiotic dosing strategies.

Results

A model that links bacterial population dynamics with molecular mechanisms of antibiotic action

We developed a linear dynamical model to describe the effect of antibiotic exposure on the growth and death rates of a bacterial population (Fig. 1A) (see Methods, Model formulation and analysis for a mathematical description of the model). We assume that each bacterial cell in the population carries an identical number N of intracellular proteins that the drug targets for inactivation. Drug molecules inactivate target proteins by binding to them with a rate k and can dissociate from the target with a rate k. The affinity KD of the drug is thus the ratio of off-rate to on-rate, KD = k/k. The model assumes that the growth and death rates of a bacterial cell depend on its drug-target occupancy (that is, the number of inactivated drug-target complexes it contains) [34], [35]. We denote drug-target occupancy with the index i, which ranges from 0 to N. The model is a system of N + 1 ordinary differential equations; the ith equation of the system describes the change in the size of the bacterial subpopulation with i inactivated drug-target complexes as a function of time. Cells harboring successively larger numbers of inactivated drug-target complexes have successively faster death rates and/or slower growth rates, depending on the mechanism of action of the drug (see Results, Classification of drug action). We thus define the growth rate (G[i]) and death rate (D[i]) of each subpopulation as discrete monotonic functions of drug-target occupancy. In practice, G[i] and D[i] take the form of constrained logistic functions each controlled by a steepness and inflection point parameter, allowing us to define quasi-linear, quasi-stepwise, quasi-exponential, and sigmoid curves (Supplementary Fig. S1).

Fig. 1

Features of a model that links bacterial population dynamics with the molecular mechanisms of antibiotic drug action. (A) Illustration of the model. We consider a population B of bacterial cells harboring i inactive drug-target complexes. The change in the size of B is a function of cellular growth and death rates (each of which is determined by the value of i, Supplementary Fig. S1), and of the molecular kinetics of the drug binding and unbinding to its protein target. The total bacterial population is given by the sum B + B + … + B-1 + B, where N is the number of drug targets per cell. (B) Dynamics of a bacterial population exposed to a drug dose above the minimum inhibitory concentration (MIC). The black line represents the total bacterial population; shaded lines represent subpopulations with X and fewer inactivated drug-target complexes. Population dynamics as a function of drug concentration are shown in Supplementary Fig. S2. (C) Proportion of the bacterial subpopulation B as a share of total population for the first three hours of the curve shown in panel (B). (D) Pharmacodynamic curves derived from the model for a wild-type (light purple) and drug-resistant (dark purple) bacterial strain. The MIC is denoted as the drug concentration at which the net bacterial growth rate is zero. Inset: the resistance selection window (purple shading) is given by the drug concentration range within which the drug-resistant strain exhibits a higher—but still positive—net growth rate compared to the wild-type strain. G denotes the growth rate of the wild-type strain in the absence of antibiotic (i.e. the growth rate for subpopulation B). D denotes the maximum death rate of bacterial strains when all N cellular targets are inactivated (i.e. the death rate of subpopulation B). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The model tracks the growth and death of all N + 1 bacterial subpopulations, each denoted B, over time (Fig. 1B). Drug concentration determines the net growth rate of the entire bacterial population (Supplementary Fig. S2). In the absence of drug, the population grows exponentially at a rate equal to the difference between the drug-free growth and death rates (G and D, respectively). When drug is present, the composition of bacterial subpopulations asymptotes towards a steady state after a transient phase during which drug molecules bind to their targets (Fig. 1C). At steady state, the relative composition of bacterial subpopulations does not depend on the total size of the population.

We can calculate the MIC of a drug directly from model parameters when assuming a drug concentration that is constant in time (see Methods, Calculation of the minimum inhibitory concentration). We can also simulate clinically observed drug resistance mutations by modulating the parameters of the model that influence the value of the MIC. Changes in the binding kinetics of the drug (i.e. k and k) simulate target modification mutations that decrease the affinity of an antibiotic molecule to a cellular target [36], [37], [38]. Changes to the value of N represent changes in the number of protein targets per cell, equivalent to target up- or downregulation [39], [40], [41]. We assume that fitness costs associated with resistance alleles take the form of reduced growth rates, and we simulate this cost by reducing the drug-free growth rate of resistant strains by a factor c such that the maximum growth rate of a resistant strain (G) relative to that of a wild-type strain is G = G(1–c). When c ranges from 0 (no cost) to 1 (no growth), the resistant strain exhibits a slower growth rate relative to that of the wild-type. If c is negative, the resistant strain exhibits a faster drug-free growth rate than does the wild-type strain, as has been observed in rare cases with some fluoroquinolone-resistant Escherichia coli isolates [42]. The model also enables us to generate pharmacodynamic curves by calculating the net growth rates of simulated bacterial populations over a range of drug concentrations (Fig. 1D). The resistance selection window constitutes the range of drug concentrations over which a drug-resistant mutant strain has a higher but strictly positive net growth rate relative to that of its wild-type counterpart (Fig. 1D, inset). By leveraging biochemical descriptions of drug-target kinetics to simulate the growth and death of bacterial populations under antibiotic exposure, our approach enables us to model a diversity of antibiotic mechanisms of action, bacterial mechanisms of resistance, and clinically relevant pharmacodynamic parameters.

Inferring antibiotic mechanisms of action from population-scale data

To test the utility of our biochemical model for gaining cellular-scale insights into antimicrobial drug mechanisms from population-scale experiments, we calibrated our model to experimental time-kill curves of the gram-negative bacterium Escherichia coli challenged to ciprofloxacin, a fluoroquinolone first brought to market in 1987, and ampicillin, a β-lactam introduced in 1961. Ciprofloxacin has two known molecular targets in bacteria, both of which are heterotetrameric type-II topoisomerases: the DNA gyrase complex (GyrA2B2) and DNA topoisomerase IV (ParC2E2). However, ciprofloxacin preferentially binds to the GyrA2B2 complex in gram-negative bacteria [43]. We used a mass-spectrometry based estimate for the number of GyrA2B2 complexes per E. coli cell (N ∼ 183) as the number of drug targets within each bacterium [44]. The targets of ampicillin are the penicillin binding proteins (PBPs), which play critical roles in peptidoglycan synthesis [45]. Ampicillin and other β-lactams inactivate PBPs by acylating a catalytic serine residue. Of the > 10 PBPs that have been described in E. coli, the high-molecular mass (HMM) PBPs (1a/1b/2/3) are known to play essential physiological roles [46]. When fitting our model to ampicillin time-kill curve data, we used a measurement for the number of HMM-PBPs per E. coli cell (N ∼ 600) as the number of drug targets within each bacterium [47].

We implemented an adaptive simulated annealing algorithm to calibrate the parameters of our model to time-kill curves (Methods, Model calibration via simulated annealing). Simulated annealing is a global optimization method that seeks to minimize the value of an objective function (in our case, the difference between experimental data and model predictions) by modifying model parameters according to a probability distribution that resembles a Boltzmann distribution. At each iteration, the algorithm accepts a new set of model parameter values with a probability defined by the objective function value given by the parameter set from the previous iteration, the objective function value given by the parameter set from the current iteration, and a characteristic system temperature that controls the probability with which the algorithm accepts a new parameter set. As the temperature of the system is reduced with each iteration, the algorithm converges on a parameter set that minimizes the difference between experimental data and model predictions (see Supplementary text for a more detailed description of the algorithm) [48]. For the ciprofloxacin dataset, we performed 249 independent parameterizations using the algorithm and selected the parameter set that yielded the lowest objective function value (Fig. 2A, Table 1, Supplementary Fig. S3). Bacterial persistence [49], [50] likely plays a role in the slower-than-expected population decline that we observe experimentally at high drug concentrations. At antibiotic doses below those that elicit persistence, the calibrated model accurately recapitulates the pharmacodynamic curve derived from experimental data (Supplementary Fig. S4).

Fig. 2

Calibrating the model to experimental data reveals underlying mechanisms of drug action. (A) Comparison between calibrated biochemical model (solid lines) and experimental data (shaded points). The experimental data (Supporting Data File S1) represent time-kill curves of Escherichia coli exposed to ciprofloxacin. A summary of all independent model calibrations is shown in Supplementary Fig. S3. (B) Comparison of the calibrated biochemical model with the E pharmacodynamic model [33]. We fit the E model to the same experimental dataset shown in panel (A) and compared Pearson correlation coefficients (R2) and MICs. Red points in the MIC panel denote experimentally-measured ciprofloxacin MICs for E. coli strains isolated prior to the widespread emergence of quinolone resistance (Supporting Data File S2). The solid horizontal line represents the mean of experimental measurements, and the dashed lines indicate the 95% confidence interval. A comparison of the pharmacodynamic curves obtained from the models is shown in Supplementary Fig. S4. (C) Cellular growth and death rates as a function of ciprofloxacin-GyrA2B2 complex number (i) for the model calibrated to the experimental data shown in panel (A). (D) Four extreme schemes of drug action resulting from two characteristics (activity and steepness) of a drug’s effect on growth and death rates as a function of drug-target occupancy. Supplementary Fig. S5 shows the simulated bacterial kill curves for these schemes at 4x MIC. Model fits for drug-free growth rate (G) and drug-saturated death rate (D) are shown in Supplementary Fig. S6. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1

Model parameterization to ciprofloxacin time-kill curves. We obtained the values of k, k, α, α, γ, and γ by calibrating the model to experimental data (Fig. 2). We inferred antibiotic-free growth rate and antibiotic-saturated death rate (G and D) by fitting an exponential curve to ciprofloxacin kill curves using 0 and 2.19 µg/ml of drug, respectively (Supplementary Fig. S6). We use a constrained logistic function to model the growth and death rates of bacterial cells as a function of inactivated target number, where α controls the steepness of the logistic function and γ controls the inflection point of the logistic function (Supplementary Fig. S1). Parameters not obtained from model calibrations to experimental data were retrieved from the literature. For the bacterial death rate in the absence of drug (D), we used the mean of death rates reported in Wang et al., 2010 [51].

Model parameters
Name	Description	Value	Units	Source
N	Number of target proteins per cell (i.e. GyrA2B2 copy number)	183	cell−1	[44]
G0	Bacterial growth rate in the absence of drug	0.526	hr−1	Fit by model
GN	Bacterial growth rate in saturating concentrations of drug	0	hr−1	Fit by model
D0	Bacterial death rate in the absence of drug	5.40 × 10−3	hr −1	[51]
DN	Bacterial death rate in saturating concentrations of drug	7.53	hr −1	Fit by model
kF	Drug-target binding rate	5.23 × 103	M −1sec−1	Fit by model
kR	Drug-target unbinding rate	3.17 × 10−4	sec−1	Fit by model
αG	Steepness of growth rate function G[i]	16.8	# drug-target complexes−1	Fit by model
αD	Steepness of death rate function D[i]	7.29	# drug-target complexes−1	Fit by model
γG	Inflection point of growth rate function G[i]	24.9	# drug-target complexes	Fit by model
γD	Inflection point of death rate function D[i]	359	# drug-target complexes	Fit by model
B	Initial size of bacterial population at the start of drug treatment	Varies	cell ml−1	n/a
µR	Mutation rate for drug resistance emergence	2.00 × 10−7	cell−1 division−1	[52], [53]
µC	Mutation rate for emergence of secondary mutations in resistant strains	2.00 × 10−6	cell−1 division−1	[52], [53]
cR	Cost of resistance mutation, such that the antibiotic-free growth rate of a resistant mutant is G0 (1 - cR)	0.25	Non-dimensional	[54]

We compared our biochemical model’s ability to capture critical pharmacodynamic characteristics of a drug against that of an E model [33]. The E approach describes net bacterial growth rate directly as a function of drug concentration and does not accommodate molecular descriptions of drug-target interactions. Such models have been used extensively to estimate pharmacodynamic parameters, to design drug dosing regimens, and to predict the strength of resistance selection at nonzero drug concentrations [17], [55]. Our formulation delivers performance comparable to that of the E model for fitting experimental time-kill curves (Fig. 2B, left panel) and more accurately estimates MIC (which we calculated to be 8.8 × 10−3 µg/ml for ciprofloxacin) from these data (Fig. 2B, right panel). We note that the eigenvalue-based method for estimating MIC using our model seeks to calculate the drug concentration at which a bacterial population undergoes zero net growth at infinite time and can thus be interpreted as a lower bound on the experimental MIC (see Methods, Calculation of minimum inhibitory concentration). If we use the Clinical & Laboratory Standards definition of MIC as the concentration of drug that yields zero net growth at 18 h, our model predicts an MIC of 9.4 × 10−3 µg/ml. This is within 20 % of the mean of experimental MIC measurements.

Our model furthermore offers capabilities that the E approach lacks, including the ability to estimate molecular kinetic parameters of drug-target binding from population-scale data. We estimated the gyrase-ciprofloxacin unbinding rate (k) to be 3.17 × 10−4 sec−1, near the experimentally measured value of 3 × 10−4 sec−1 [56]. We also analyzed the KD values for ciprofloxacin binding to E. coli GyrA2B2 generated for the 249 independent parameterizations described above. As our fitting method is stochastic, not all model calibrations reach local minima. However, the best 90 % of all calibrations (that is, the 224 fits with the lowest objective function values) consistently converged upon a narrow range of affinity values (95 % confidence interval: 7.2 × 10−8 to 1.6 × 10−7 M) (Supporting Data File S3). Our estimates lie within the range of KD values of ciprofloxacin for E. coli GyrA2B2 reported from experimental measurements, which span from 3.2 × 10−8 to 3.0 × 10−6 M [57], [58], [59], [60] (see Supplementary Text for a discussion of the convergence of other model parameters). These results suggest that our model can estimate molecular kinetic parameters from population-scale data, but that numerous simulated annealing runs are required to robustly identify a global minimum within parameter space.

We also fit our model to time-kill curves of E. coli exposed to ampicillin (Supplementary Table S1, Supplementary Fig. S7). Ampicillin acts predominantly as a bactericidal agent [61]; we thus assumed that bacterial growth rate is invariant to the number of acylated HMM-PBPs per cell (G[i] = G for all values of i) and fit the death rate function (D[i]) to experimental data. Because β-lactams acylate HMM-PBPs at different rates, we also used experimental measurements for the acylation and deacylation rates of E. coli PBP1b exposed to ampicillin [62]. As with the ciprofloxacin time-kill curves, we observed persistence at high antibiotic concentrations, but our calibrated model accurately recapitulates experimental data at lower (≤96 µg/ml) drug concentrations. The death rate function inferred by the model is similar to that inferred by COMBAT on separate experimental replicates of ampicillin time-kill curves [34], suggesting that our linear model can fit parameters to experimental data as robustly as COMBAT. These results indicate that our model can fit experimental time-kill curves of bacterial populations exposed to antibiotics with distinct mechanisms of action, predict MICs in close agreement with experimental measurements, and recapitulate results inferred by more computationally intensive nonlinear models.

Classification of antibiotic action

Another unique feature of our approach is the ability to describe bacterial growth and death rates as a function of drug-target occupancy. For ciprofloxacin, the calibrated model predicts three regimes of bacterial subpopulation dynamics in relation to GyrA2B2 inactivation: a growth regime in which bacterial replication dominates among subpopulations with low numbers of inactivated targets, a stalling regime for intermediate numbers of drug-target complexes in which neither growth nor death is appreciable, and a killing regime at high numbers of inactivated targets in which bacterial death increases quasi-exponentially (Fig. 2C). The forms of G[i] and D[i] that we obtain here suggest that ciprofloxacin has a multimodal mechanism of action, a result consistent with prior experimental studies [43], [63], [64] and with COMBAT [34]. The drug stalls cellular replication at intermediate target occupancies and induces killing only at higher doses. Like many antibiotics, ciprofloxacin thus exhibits both bactericidal and bacteriostatic effects on microbial populations [64], [65]. Our biochemical model represents this explicitly.

Most drugs nonetheless demonstrate a greater degree of bactericidal or bacteriostatic activity at clinically relevant doses [66], and we hypothesized that the ability of a drug to stall growth or to accelerate death may affect the selection for resistant strains and the emergence of secondary mutations. We also suspected that the relationship between drug-target occupancy and antibiotic effect—reflected in the steepness of the G[i] and D[i] functions—could further shape the dynamics of resistance selection.

These two characteristics (bactericidal versus bacteriostatic activity and drug effect steepness) represent two general dimensions along which a drug’s mechanism of action can affect the growth and death of bacterial populations. Four extreme cases of drug action thus exist (Fig. 2D). In the case of a purely bacteriostatic antibiotic, death rates are a constant function of inactivated drug-target complex number (that is, D[i] = D for all values of i). For a purely bactericidal antibiotic, the growth rate of all bacterial subpopulations remains constant (G[i] = G for all values of i). The steepness of the drug effect is reflected in the form of the function D[i] for bactericidal antibiotics and G[i] for bacteriostatic antibiotics (Supplementary Fig. S1). We defined linear and stepwise onset of action as our two extremes, as other monotonic forms are intermediate cases of these curves.

The opposing effects of increased drug resistance and decreased cellular fitness

Strains with mutations that confer resistance against antibiotics often have reduced growth rates compared to those of drug-susceptible strains [10], [11]. The balance of replication and death determines bacterial net growth both in the absence and in the presence of antibiotics, and very high fitness costs associated with resistance can prevent bacterial viability at any drug concentration [67]. We sought to elucidate the quantitative basis for the trade-off between drug resistance and cellular growth and to investigate how the drug mechanisms defined above influence the range of permissible fitness costs that a drug-resistant mutant can incur while still maintaining a drug susceptibility that is lower than that of a wild-type strain. In the simplest case of the model, where the number of target molecules per cell is 1, the expression for the MIC captures the opposing effects of drug resistance and cellular growth (see Methods, Calculation of minimum inhibitory concentration for derivation):

The MIC increases with reductions of the on-rate kinetics of drug-target binding (k) and with increases in the off-rate kinetics of drug-target binding (k), but decreases with fitness costs that manifest as reductions in the drug-free growth rate (G). These proportionalities hold for any number N of drug targets.

We modeled the opposing effects of biochemical changes that reduce drug susceptibility (i.e. altered drug-target binding kinetics or target overexpression) and the fitness costs of these biochemical changes. To model different drug mechanisms, we considered a set of six antibiotics (Supplementary File S2, Supplementary Fig. S5). Two antibiotics in the set feature growth and death dynamics derived from the model calibration to ciprofloxacin and ampicillin time-kill curves. The other four antibiotics are hypothetical and feature growth and death dynamics that represent four extremes of antibiotic action (Fig. 2D). These hypothetical drugs use molecular kinetic parameters (N, k, k) identical to those used for ciprofloxacin simulations. We simulated mutant strains of E. coli that acquire drug resistance phenotypes either through changes in the molecular kinetics of drug binding (k or k) or by increasing the copy number N of the drug’s cellular protein target. Each of these resistance mechanisms has been observed in clinical isolates of drug-resistant, gram-negative bacteria [11], [29], [68]. We then simulated fitness costs associated with the resistance mutation and calculated the mutant strain’s MIC relative to that of the wild-type strain.

For resistance acquired through changes in the kinetics of drug-target binding (k and k), we found that mutants can tolerate strikingly high (80–99 %) fitness costs while still maintaining an MIC that is greater than that of the drug-susceptible wild-type (Fig. 3, top and middle rows). This permissibility of fitness costs exists for all six of the drug mechanisms we simulated, although ampicillin and drugs that act with linear effects (Bacteriostatic/Linear and Bactericidal/Linear) have a narrower range of permissible fitness costs than do ciprofloxacin and drugs that act with stepwise effects. For all drug mechanisms, mutant strains make larger gains in MIC by decreasing the on-rate kinetics of drug-target binding (k) than they do by increasing the off-rate kinetics of drug-target binding (k) by the same amount (Supplementary Fig. S8). That is, mutations that lead to the same change in drug-target affinity (as quantified by the dissociation constant KD = k/k) through different changes in the on– and off-rate binding kinetics do not necessarily have the same range of permissible fitness costs. This has biological significance—limiting the opportunity for a drug to bind to its target, thereby preventing the drug from actuating its effects on cellular growth and death, should lead to lower drug susceptibilities than would accelerating the rate at which an already-formed drug-target complex disassociates. The difference in the fitness effects of mutations that modify k and k is especially pronounced for ampicillin, where a larger number of target molecules results in smaller incremental changes to the death rate with each successive drug deacylation event. We conjecture that affinity-reducing mutations to target proteins within resistant bacterial strains may occur preferentially through reductions in on-rate kinetics, especially for bactericidal agents that target abundant cellular components (such as PBPs or ribosomes), although further experimental studies will be required to validate this hypothesis.

Fig. 3

Drug mechanism influences the fitness landscapes of resistance mutations. We calculated the MIC, expressed as a fold-change relative to the MIC of the wild-type, for mutant strains carrying (top row) drug targets with reduced binding kinetics (k), (middle row) drug targets with accelerated unbinding kinetics (k), or (bottom row) increased numbers of drug target molecules (N). Mutant strains also carry fitness costs, expressed as a fractional reduction in drug-free growth rate relative to wild-type. Cost-free MIC as a function of k and k for all mechanisms of action are shown in Supplementary Fig. S8.

Biochemical experiments have shown that ciprofloxacin exhibits a bactericidal effect by permitting GyrA2B2-mediated cleavage of DNA but preventing DNA re-ligation, resulting in the formation of toxic DNA double-stranded breaks [43], [69]. When simulating the overexpression of target proteins in resistant cells (Fig. 3, bottom row) we therefore assumed that bacterial killing increases when a fixed number of inactivated drug-target molecules form within a cell (that is, we assume a toxicity threshold whereby γ remains constant with changing N). This assumption is consistent with experiments that have observed a positive correlation between the number of quinolone-stabilized DNA-gyrase complexes and the rate of bacterial killing [70]. Conversely, we assumed that a resistant cell requires a fixed number of active, non-complexed target proteins in order to maintain its maximum growth rate (that is, a survival threshold). γ thus changes in step with N such that N-γ remains constant. We made these same assumptions for the four hypothetical antibiotics. Contrary to ciprofloxacin, ampicillin does not stabilize a toxic molecular intermediate and instead exerts a bactericidal effect by inhibiting the transpeptidase activity of PBPs [46], [71]. When assessing the effects of PBP overexpression on ampicillin resistance, we therefore simulated a survival threshold on the death rate function whereby a constant number of uninhibited targets is required to maintain a minimum death rate.

We found that target overexpression has a diversity of effects on resistance that depend on the mechanism of action of the drug. Overexpression leads to substantial gains in resistance against bacteriostatic drugs that exhibit stepwise effects, even at very high fitness costs. The effect of target overexpression on drug resistance is negligible for bactericidal drugs and for bacteriostatic drugs with a linear effect on growth stalling. For ciprofloxacin, small (2–3-fold) increases in target number can lead to modest increases in MIC. However, larger increases in target number lead to increased antibiotic susceptibility. This result is consistent with experimental studies on target amplification, in which the overexpression of gyrAB in E. coli resulted in increased susceptibility to ciprofloxacin but did not induce reductions in growth rate [40].

Our model also suggests that, in the absence of fitness costs, PBP overexpression confers increased resistance to ampicillin. Increased resistance due to target overexpression is an expected—and clinically observed—phenomenon for antibiotics whose primary mode of action is simple protein inhibition without toxic product formation [71], [72]. However, experimentation has shown that E. coli exhibits reduced growth rates when overexpressing PBPs and that these fitness costs lead to reduced β-lactam resistance when PBPs are overexpressed beyond certain thresholds [40]. Indeed, our model suggests that the MIC of a strain with increased PBP expression can be lower than that of the wild-type if fitness costs are sufficiently high. Our model thus recapitulates the experimentally observed effects of target overexpression on resistance to two clinically relevant antibiotics with distinct mechanisms of action: Overexpression confers reduced resistance to drugs that induce toxic intermediates even in the absence of fitness costs (as is the case for ciprofloxacin), and overexpression confers reduced resistance to drugs that inhibit protein function only in the presence of fitness costs (as is the case for ampicillin).

Drug mechanism shapes the resistance selection window

To understand how a drug’s mechanism of action affects the propensity to select for resistance during treatment, we simulated the pharmacodynamics of wild-type and drug-resistant strains challenged to each of the six drugs in the set outlined above. MICs for clinical isolates of ciprofloxacin-resistant E. coli strains with single point mutations in GyrA, which may reduce the affinity of ciprofloxacin to GyrA2B2, range from 10 to 16 times greater than the MIC of a drug-susceptible wild-type [36], [68], [73], [74]. Data on the fitness costs associated with mutant GyrA-mediated ciprofloxacin resistance in E. coli are sparse, but studies of rifampicin-resistant clinical isolates of Mycobacterium tuberculosis with point mutations in the rpoB gene have suggested that a 20–30 % reduction in growth rate is approximately the maximum fitness cost that drug-resistant mutants can incur before facing extinction in competitive drug-free environments [54]. To model drug-resistant strains, we therefore scaled k and k such that the MIC of the resistant strain is 12 times that of its drug-susceptible counterpart given a 25 % fitness cost (c = 0.25) (Fig. 4A). We applied the same strategy to ampicillin-resistant strains.

Fig. 4

The propensity to select for resistant mutants depends on drug mechanism. (A) We modeled wild-type strains using the parameters obtained from the model fits to ciprofloxacin (Fig. 2) and ampicillin (Supplementary Fig. S7) time-kill curves. (B) Relationship between MICs of resistant strains (expressed as multiples of MICWT) and fitness cost of resistance. Horizontal dashed lines indicate the MICs of the wild-type and resistant strains described in panel (A); the vertical dashed line indicates the fitness cost at which all resistant strains have the same fold-increase in MIC relative to that of wild-type (c = 0.25). (C) Pharmacodynamic curves for the wild-type and resistant strains described in panel (A). (D) Resistance selection windows for drug-resistant strains. The fitness advantage of resistant strains over wild-type strains is shown within the drug concentration range in which the resistant strain has a positive net growth rate that is larger than that of the wild-type. The fitness advantage is expressed as a proportion of the resistant strain’s growth rate in the absence of drug (G). Supplementary Fig. S9 illustrates the relationship between the size of the resistance selection window and the steepness of a drug’s pharmacodynamic curve. CIP: ciprofloxacin; AMP: ampicillin; S/S: bacteriostatic/stepwise effect; S/L: bacteriostatic/linear effect; C/S: bactericidal/stepwise effect; C/L: bactericidal/linear effect; MICWT: MIC of the wild-type strain; MICRES: MIC of the resistant strain.

A nearly linear relationship exists between drug resistance and fitness cost for strains resistant to drugs with a linear effect on growth or death (Fig. 4B, Bacteriostatic/Linear and Bactericidal/Linear). By contrast, drugs with stepwise effects on growth and killing (Bacteriostatic/Stepwise and Bactericidal/Stepwise) exhibit only modest reductions in MIC until they incur very high (>90 %) fitness costs. For ciprofloxacin and ampicillin, the relationships between fitness cost and drug resistance lie in between these extremes. We determined resistance selection windows for strains resistant to the six drugs in our set by simulating pharmacodynamic curves for wild-type and resistant strains (Fig. 4C). To quantify the magnitudes of selection for resistant strains, we calculated the difference in net growth rates between wild-type and susceptible strains over the concentration range that defines the resistance selection window for each drug (Fig. 4D). For linear-effect bacteriostatic drugs (Bacteriostatic/Linear), we found that the resistance selection window begins at drug concentrations as low as 200x below the MIC of the susceptible strain. Drugs with stepwise effects on growth or killing (Bacteriostatic/Stepwise and Bactericidal/Stepwise) have narrower resistance selection windows than their counterparts with more linear activity profiles.

Consistent with prior studies on the pharmacodynamic profiles of antimicrobial agents [17], [19], [75], we find that the size of the resistance selection window is associated with the steepness of a drug’s pharmacodynamic curve. Given a cellular effect (i.e. bacteriostatic or bactericidal), drugs with steeper pharmacodynamic curves tend to have narrower selection windows. This holds true for the antibiotics we investigated experimentally as well—compared to ciprofloxacin, ampicillin has both a steeper pharmacodynamic curve (as measured by a Hill coefficient) and a narrower resistance selection window (Supplementary Fig. S9). However, we also find that strains resistant to drugs with narrower resistance selection windows have higher net growth rates within the resistance selection regime than do strains resistant to drugs with wider resistance selection windows (Fig. 4D). This finding has clear clinical significance: drugs with steeper pharmacodynamic profiles feature relatively small concentration ranges that select for resistance, but the negative consequences of dosing within the resistance selection window are higher for these drugs.

The secondary mutant selection window is narrower for bacteriostatic antibiotics

The genotypic space for mutations that confer resistance to antibiotics by modifying the binding kinetics of a drug to its target, such as those described in Fig. 4, is typically highly constrained [22], [76]. Accordingly, a return to a drug-susceptible state requires reversion of the specific genetic changes that conferred resistance in a bacterial population. In contrast to resistance reversion, secondary mutation accumulation can involve a wider range of genetic changes that affect a cell’s metabolic, transcriptional, and/or translational states [23]. Therefore, the probability that a bacterial population evolves secondary mutations that compensate for the fitness costs of a resistance mutation is often higher than the probability that a bacterial population will revert to susceptibility in drug-free environments [20], [77]. During treatment, resistant bacterial populations may also accumulate secondary mutations that further raise MIC. In order to understand how drug mechanism influences such secondary adaptation, we simulated the emergence of secondary mutants from drug-resistant subpopulations of a bacterial population faced with antibiotic challenge (Supplementary Fig. S10; Methods, Simulating the emergence of secondary mutations).

At a given antibiotic concentration, the probability of secondary mutation emergence is substantially higher for drugs with linear effects on cellular growth and death than it is for drugs with stepwise effects (Fig. 5A). This holds true for both bactericidal and bacteriostatic agents. Counterintuitively, then, the suppression of secondary mutation emergence is not necessarily guaranteed by rapid killing as suggested by earlier studies [78]. Likewise, rapid attenuation of cell division does not halt the emergence of secondary mutations. We studied the basis for this result by investigating the steady-state target occupancy distributions of cells under antibiotic exposure. By accounting for the kinetics of drug-target binding, our biochemical model shows that target occupancy among cells follows a distribution and is not a single value even in otherwise clonal bacterial subpopulations (Fig. 5B). This results in heterogeneous replication rates within the drug-resistant subpopulation (Supplementary Fig. S11) that allow some bacteria to mutate. Classical population-dynamic models of antibiotic action [33], [78], which assume that a drug affects the net growth rate of all cells equally, overlook this phenomenon.

Fig. 5

Emergence of secondary mutations among resistant subpopulations of infecting bacteria. (A) Probability of a drug-resistant strain with secondary mutations emerging from an infecting bacterial population before the infection is cleared (i.e. before the total bacterial population decreases to less than 1, Supplementary Fig. S10). The initial population size for this simulation is 109 cells. Inset shows probabilities of secondary mutation emergence before infection clearance when the drug concentration used is 2x MICRES. (B) Frequency distributions of inactive drug-target complexes for drug-resistant subpopulations undergoing steady-state exponential decline at 2x MICRES. Growth and death rate distributions for these populations are shown in Supplementary Fig. S11. (C) Probability of secondary mutant emergence from bacterial subpopulations with i inactivated drug-target complexes, shown for ciprofloxacin and ampicillin dosed at 2x MICRES. (D) Probability of secondary mutant emergence from bacterial subpopulations as a function of drug dose, shown for ciprofloxacin. Probabilities are shown as absolute values (left panel) and as values normalized to the total probability of compensation for the entire bacterial population over the course of treatment (right panel). (E) Resistance and secondary mutant selection windows for different drug action mechanisms. The resistance selection window (middle purple) is defined as the drug concentration range over which a drug-resistant strain has a growth advantage over the wild-type. The secondary mutant selection window (dark purple) is defined as the drug concentration range over which the probability of a resistant strain with secondary mutations emerging before infection clearance exceeds 10−4 (see Supplementary Fig. S12 and Methods, Simulating the emergence of secondary mutations). Dashed lines indicate the MICs of the wild-type and resistant strains. CIP: ciprofloxacin; AMP: ampicillin; S/S: bacteriostatic/stepwise effect; S/L: bacteriostatic/linear effect; C/S: bactericidal/stepwise effect; C/L: bactericidal/linear effect; MICWT: MIC of the wild-type strain; MICRES: MIC of the resistant strain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

For ciprofloxacin and ampicillin doses only slightly above the MIC of the resistant strain ([Drug] = 2x MICRES), we found that secondary mutations are most likely to emerge once the bacterial population has reached a steady-state target occupancy distribution (Fig. 5C). For ciprofloxacin, a considerable probability of secondary mutation emergence nonetheless exists among bacterial subpopulations with low numbers of inactivated drug-target complexes. These low-occupancy subpopulations have faster growth rates and thus higher mutation rates (for ampicillin, where we assume only bactericidal action, all subpopulations have equal growth rates and thus equal mutation rates). Low-occupancy subpopulations are also present in very large numbers during the initial stages of treatment, when drug molecules are binding to their cellular targets and before the overall population begins to decline (Fig. 1C). We found that drug concentration influences the likelihood of a secondary mutant arising from a steady-state or a low-occupancy subpopulation (Fig. 5D). While the overall probability of secondary mutation emergence decreases with higher drug dose (Fig. 5D, left panel), the relative probability that a secondary mutation arises from a low-occupancy population is greater for higher drug doses (Fig. 5D, right panel). This implies that secondary mutations are more likely to emerge very early during treatment (i.e. within the first few hours) when high doses of drugs with bacteriostatic effects are used.

Prior studies have estimated that the probability of the existence of a fitness cost-free bacterial pathogen prior to treatment ranges from 5 × 10−5 to 3 × 10−4 per infection [79]. We sought to determine the range of drug concentrations over which the likelihood of secondary mutation emergence during treatment is at least as high as the likelihood for preexisting secondary adaptation. We therefore determined the drug concentration at which the probability of secondary mutation emergence before population extinction equals 10−4 (that is, each treatment course has a 1 in 10,000 chance of giving rise to a resistant strain bearing secondary mutations). We used this value as an upper boundary for the “secondary mutant selection window,” the range of drug concentrations over which the probability of the emergence of a drug-resistant bacterial strain with secondary mutations is substantial (Supplementary Fig. S12). The secondary mutant selection window extends the range of drug concentrations defined by the resistance selection window over which drug-resistant strains may be selected (Fig. 5E).

As with the resistance selection window, we found that the size of the secondary mutant selection window varies dramatically depending on a drug’s mechanism of action. The secondary mutant selection window for ampicillin is nearly-two orders of magnitude larger than that for ciprofloxacin. Drugs that fully suppress cellular replication above MIC (i.e. Bacteriostatic/Stepwise) have small secondary mutant selection windows, as the probability that additional mutations emerge over the course of treatment is equal to the probability that a resistant strain with secondary mutations emerges during the transient phase of drug-target binding immediately after treatment onset, which lasts on the order of a few hours (Fig. 1C).

We find that, for a given effect steepness (linear or stepwise), bacteriostatic drugs have narrower secondary mutant selection windows than do bactericidal drugs. This is expected given that secondary mutations arise from actively replicating bacterial populations. This result nonetheless contrasts with the characteristics of the resistance selection window, where we observed that drug effect steepness plays a stronger role in determining the size of the window than does bacteriostatic or bactericidal effect (Fig. 4D). We therefore conclude that drugs with narrow resistance selection windows are not necessarily effective at suppressing secondary mutation. This is indeed the case for the antibiotics we investigated experimentally—ampicillin has a slightly narrower resistance selection window than does ciprofloxacin, but a considerably wider secondary mutation selection window.

Comparison with other modeling approaches

The model we describe in this study is a linear case of COMBAT, a previously reported formulation for simulating the growth and death of bacterial populations under antibiotic exposure using chemical reaction kinetics [34]. To assess the precision and performance of our model relative to COMBAT, we simulated pharmacodynamic curves with both COMBAT and our linear model using identical parameter sets and drug concentrations (Supplementary Fig. S13A). With both models, we simulated bacterial populations for two hours and calculated net growth rates. We found that both models predicted similar net growth rates across a range of drug concentrations, and that they predicted MIC values within 0.01 % of one another. However, our linear model simulated its pharmacodynamic curve > 1000 times more rapidly than did COMBAT (mean computation time of 0.017 sec for the linear model vs 46 sec for COMBAT) (Supplementary Fig. S13B). Additionally, COMBAT relies on numerical schemes to solve a set of partial differential equations. Its computation time therefore increases with increasing simulation time. Because our model relies on an eigenvector representation of a system of ordinary differential equations to simulate bacterial populations under a constant drug concentration, computation time does not increase with increasing simulation time (Supplementary Fig. S13C). Because of these attributes, we recommend the linear model reported in this study for fitting experimental data where drug concentration can be assumed to be constant (such as laboratory-derived time-kill curves) and for simulation tasks that require the exploration of parameter space at high resolution (such as those reported in Fig. 3).

Discussion

The increasing prevalence of first line- and multi-drug resistant bacteria [1], [2] signals the need for new antibiotics and robust drug dosing strategies that minimize the emergence and spread of resistance [4]. Despite this need, little is known about the role that a drug’s mechanism of action plays in the evolution of antibiotic resistance. We studied the relationship between drug mechanism and drug resistance using a mathematical model that connects bacterial population dynamics with molecular-scale descriptions of drug-target inactivation kinetics (Fig. 1A). Our biochemical model allows us to describe bacterial replication and death as functions of drug-target inactivation, enables us to estimate molecular kinetic parameters from population-scale data, and delivers performance on par with that of classical pharmacodynamic models (Fig. 2B).

We calibrate the model to experimental datasets of ciprofloxacin and ampicillin time-kill curves, and we show that drug-resistant strains can incur strikingly high fitness costs associated with mutations that reduce drug-target binding kinetics (Fig. 3). We find that the relationship between drug-target inactivation and antibiotic effect (i.e. bacterial killing, growth stalling, or both) exerts a strong influence on the strength of selection for resistant strains during treatment, regardless of whether the drug is bactericidal or bacteriostatic (Fig. 4D). We also show that the molecular kinetics of drug-target binding within cells results in heterogeneous replication rates among members of an otherwise homogenous population (Fig. 5B). This enables some drug-resistant strains to develop secondary mutations that can further reduce drug susceptibility, increase resilience in drug-free environments, and ultimately lead to treatment failure.

The clinical consequence of the frequently-observed trade-off between bacterial fitness and drug resistance [10] is the existence of a resistance selection window—a range of drug concentrations that selects for the propagation of drug-resistant strains over their drug-susceptible counterparts [5], [15]. It is important to note that numerous factors not captured by the resistance selection window can contribute to resistance selection in clinical settings, most notably ecological interactions between drug-susceptible strains, drug-resistant strains, and host physiology [80]. Our approach nonetheless enables us to isolate the roles that a drug’s mechanism of action play in driving the emergence of resistance.

We show that the resistance selection window is narrower for drugs that exert their effects on growth or death in a stepwise (i.e. sudden) manner, resulting in a steeper pharmacodynamic curve (Fig. 4C-4D, Supplementary Fig. S9). This result is consistent with other studies on the pharmacodynamics of antimicrobial agents, which have found that the size of the resistance selection window is associated with the steepness of the pharmacodynamic curve [17], [19], [75]. The characteristics of antimicrobial agents that enable steeper pharmacodynamic curves nonetheless remain poorly described. Models that capture the effects of antibiotic drugs on multiple scales, such as that described in this study and elsewhere [34], [35], could serve as helpful tools for studying the interplay between a drug’s molecular mechanism and its effect on bacterial population dynamics, enabling the design of new antimicrobial agents with narrow resistance selection windows.

Mutations that alleviate the fitness costs associated with drug resistance and/or that further raise a strain’s MIC play a major role in driving the spread of antimicrobial resistance across bacterial populations and clinical settings [24]. Our study sheds quantitative light on the mechanistic factors that govern the emergence of these secondary mutations during treatment. We propose the use of the secondary mutant selection window (Supplementary Fig. S12) as a tool for assessing the likelihood of further mutation acquisition at nonzero drug concentrations. As with the size of the resistance selection window, the size of the secondary mutant selection window varies greatly depending on the mechanism of action of the antibiotic (Fig. 5E). We stress that the secondary mutant selection window does not necessarily indicate a region on the pharmacodynamic profile of a drug over which the selection of a resistant strain with secondary mutations is favored. The strength of selection depends on the physiological effect of the secondary mutation itself—that is, whether the mutation accelerates growth rate, slows drug-target binding, or exerts a multitude of other possible effects. Indeed, secondary mutations that act strictly by restoring growth rates to wild-type levels lead only to modest (usually sublinear) increases in MIC (Fig. 4B), implying that strains with cost-free resistance phenotypes would still have MICs well below the upper boundary for the secondary mutant selection windows shown in Fig. 5E. Rather, the secondary mutant selection window defines the drug concentration range within which the accumulation of secondary mutations is substantial and therefore clinically significant.

Suppressing secondary mutation is nonetheless crucial for reducing the survival of drug-resistant mutants in antibiotic-free environments, where drug-resistant strains enter into direct competition with other microbial organisms for limited resources [10], [23]. We demonstrate that dosing drugs at or slightly above the MIC of a resistant strain may not be sufficient for preventing the spread of resistance, and that—for purely bactericidal drugs and/or drugs with low bacteriostatic potency—there exist appreciable risks of selecting for secondary mutations even at doses substantially above the MIC of the resistant strain. Reassessing the range of drug concentrations that selects for resistant mutants as a composite of the resistance selection window and the secondary mutant selection window (Fig. 5E, Supplementary Fig. S12) could facilitate the design of drug dosing strategies that holistically mitigate the emergence and spread of resistance.

Our study shows that both bactericidal and bacteriostatic drugs can exhibit narrow resistance selection windows and low probabilities of secondary mutation emergence in bacterial populations subjected to antibiotic treatment. This finding challenges the long-accepted notion that bactericidal agents are superior to bacteriostatic agents in suppressing the emergence of resistance during treatment [31], and signals the need to look beyond a drug’s ability to kill or stall bacterial replication to assess the risks of resistance emergence. The relationship between drug-target inactivation and overall antibiotic effect has a much stronger influence on the strength of resistance selection than does the drug’s bacteriostatic or bactericidal activity (Fig. 4D). The processes that may dictate such a relationship for any given antibiotic nonetheless remain enigmatic. This underscores the need for deeper experimental and theoretical research on the molecular processes that govern the pharmacodynamics of antibiotic drugs.

Several similarities and differences between the model described in this study and the simpler E pharmacodynamic model [33] are worth discussing. For both ciprofloxacin and ampicillin time-kill curves, we observed that both our model and the E model were capable of fitting experimental data with high fidelity. Furthermore, our model predicts that drugs with steeper pharmacodynamic curves have narrower mutant selection windows, a pattern that has also been observed using E models [17]. However, the E model describes the steepness of the pharmacodynamic curve explicitly using a Hill coefficient, whereas our formulation models the steepness of the pharmacodynamic curve indirectly via descriptions of growth and death as a function of drug-target inactivation. Despite these phenomenological similarities, our model offers three key advantages over the E approach. First, our biochemical model can estimate molecular kinetic parameters from population-scale data. We use this capability to infer the k, k, and KD values of gyrase-ciprofloxacin interactions with good agreement to experimental measurements. Second, whereas the E model fits MIC directly to data, our model provides a strategy for calculating MIC from model parameters. This enables us to study the fitness landscapes of specific drug resistance mechanisms (such as changes in drug-target binding kinetics and target overexpression) in a manner that would be impossible to do using an E model (Fig. 3). Third, our biochemical model describes heterogenous growth and death rates within bacterial populations as a function of drug-target occupancy. This feature enables us to quantify the probabilities of mutation and fitness compensation among bacterial populations exposed to drugs with different mechanisms of action (Fig. 5).

We note that the model reported here makes several simplifying assumptions that limit its scope and generalizability. One key assumption made is that growth and death rates are monotonically decreasing and increasing functions, respectively, of drug-target inactivation. Non-monotonic dose–response curves have been described for numerous drugs since the early years of the antibiotic era [81], and these imply the existence of non-monotonic drug-target occupancy schemes or of drug-induced cellular responses (such as reduced outer membrane permeability) that lower drug-target occupancy at high antibiotic concentrations. Our model also has limitations on the scope of resistance mechanisms that it can recapitulate. While some classes of antibiotics (particularly fluoroquinolones and rifamycins) frequently elicit resistance through altered drug-target affinity, other classes elicit resistance through additional mechanisms (including drug efflux, enzymatic degradation of drug molecules, and off-target binding) not captured in the linear model presented here. Other models have been devised that link additional mechanisms of resistance (such as efflux pump activity, membrane permeability, and cellular metabolic states) with critical pharmacologic parameters (i.e. MIC) [30], [82], but do not accommodate explicit descriptions of an antibiotic’s mechanism of action. Still other models have provided valuable insights into the genotypic determinants of antimicrobial resistance fitness landscapes [83]. Finally, we note that the origins of the fitness costs of resistance mutations remain poorly understood, and models that link resistance mechanisms with mechanistic descriptions of impaired growth may yield valuable insights into the evolution of resistance traits in bacterial populations. Adapting existing models to study the relationship between antibiotic mechanism, fitness cost, and other mechanisms of resistance constitutes a promising direction for future research.

Conclusions

The proper use of antibiotics in clinical and non-clinical settings constitutes a core action for addressing the worldwide threat of antibiotic resistance [4]. The quantitative approach we present in this study may prove useful for identifying strategies that manage the emergence of resistance to existing and future antimicrobial agents. We argue that dosing regimens should account for a drug’s resistance and secondary mutant selection windows if they are to minimize the selection of resistance phenotypes during treatment. Our findings suggest that even drugs with seemingly straightforward pharmacodynamic classifications (i.e. bacteriostatic versus bactericidal action) can set bacterial populations on complex and sometimes counterintuitive evolutionary trajectories with respect to resistance selection. In the clinic, there exists little evidence that bactericidal antibiotics lead to more favorable outcomes than do bacteriostatic antibiotics, especially for combatting uncomplicated infections [65], [84]. Yet it is precisely in the treatment of uncomplicated, drug-susceptible infections that the greatest gains are to be made in mitigating the emergence of resistance. Mechanistic models such as that presented in this study can help to uncover clinically useful drug characteristics that classical models may overlook. We envision a coupling of our quantitative approach with high-throughput experimental platforms [85], [86] to aid in the development of new drugs with optimal pharmacodynamic profiles and to accelerate the discovery of drug- and pathogen-specific dosing regimens that minimize resistance emergence.

Methods

Bacterial time-kill curve experiments

We obtained time-kill curves using Escherichia coli strain BW25113 (Coli Genetic Stock Center #7636) for ciprofloxacin time-kill curve experiments and E. coli strain MG1655 (Coli Genetic Stock Center #7740) for ampicillin time-kill curve experiments [87]. We diluted liquid overnight cultures of E. coli 1:1000 into pre-warmed lysogeny broth (LB) and grew cells to an optical density at 600 nm (OD600) of 0.50. We then prepared dilution series of ciprofloxacin (highest concentration: 2.19 µg/ml) and ampicillin (highest concentration: 256 µg/ml) and added the antibiotics to bacterial cultures. We quantified bacterial population sizes at regular (10–30 min) time intervals by plating a 1:10 dilution series of liquid culture onto LB agar plates and counting colony forming units. We performed colony counting blind to plate condition, and we did not exclude any plates from the analysis. To keep shot noise below 15 % during colony counting, we quantified plates with 50 or greater colony forming units.

To further assess the biological reproducibility of our experiment, we repeated ciprofloxacin cytotoxicity assays on different days, once with a fixed timepoint measurement at 90 min post-ciprofloxacin exposure, and another with a timecourse (i.e. that presented in Fig. 2A and Supporting Data File S1). When compared at matching timepoints of drug exposure (90 min), dose–response data from these biological replicates collected on different days were highly reproducible, with Pearson correlation of 0.987, p less than 10−5. Each time the experiment was performed, counts of colony forming units before drug treatment were conducted in technical triplicate.

The time-kill curve obtained at the highest antibiotic concentration was used to determine the maximum death rate (D) of bacterial cells, and a growth curve obtained using the same protocol with the omission of antibiotic was used to determine the maximum growth rate (G) of cells in an antibiotic free environment (Supplementary Fig. S6).

Model formulation and analysis

Our biochemical model constitutes a system of linear ordinary differential equations that describe how successive numbers of inactivated drug-target complexes affect bacterial replication and death. We consider a population of initial size B of phenotypically homogenous bacteria exposed to an antibiotic. When no drug is present, bacterial cells replicate at a rate G and die at a rate D. All cells have an identical number N of proteins that drug molecules target for inactivation. We assume first-order kinetics for drug-target binding: drug molecules bind to cellular protein targets within cells, thereby inactivating the protein, at a rate k. Inactivated drug-protein targets dissociate at a rate k. The first-order affinity of the drug to its protein target (KD) is therefore the ratio of the molecular dissociation rate to the molecular on-rate (KD = k/k).

We stratify the entire bacterial population into N + 1 subpopulations according to the number i of inactivated drug-target complexes within each cell (i.e. the drug-target occupancy), and we assume that growth and death rates of each bacterial cell depend on the drug-target occupancy. That is, bacterial subpopulations with a higher drug-target occupancy have slower growth rates and/or higher death rates than do bacterial subpopulations with a lower drug-target occupancy. Growth rate is therefore a monotonically decreasing discrete function of i (G[i]), and death rate is a monotonically increasing discrete function (D[i]). We use generalized logistic equations (Supplementary Fig. S1) to describe overall growth and death rates as a function of drug-target occupancy, allowing these functions to take the form of a line, a sigmoidal curve, an exponential curve, or a step function. We assume that when a drug inactivates all N protein targets in a cell, growth rate falls to zero (for bacteriostatic drugs), death rate attains a maximal value D (for bactericidal drugs), or growth and death rates are both affected (for drugs with mixed bactericidal and bacteriostatic action). In all of these cases, the maximal rate of killing or growth attenuation can occur before all N target proteins are inactivated if, for instance, G[i] and/or D[i] are step functions with inflection points between 0 and N. During replication, a bacterial cell partitions its inactivated drug-target complexes to two daughter cells according to a binomial distribution.

If drug is assumed to remain at a constant concentration C, the change over time in the number of bacterial cells with exactly i inactivated drug-target complexes (B) depends on the growth rate G, the death rate D, and the binding kinetics of the drug to its protein target:

The first four terms on the right side of Eq. (2) describe changes in B due to drug-target binding and unbinding. The fifth term describes bacterial death, the sixth term describes bacterial growth, and the seventh term describes the partitioning of drug-target complexes upon replication according to a binomial distribution. Eq. (2) is a linear form of a model we have described previously that treats drug-target complex number as a continuous variable rather than as a natural number [34]. Linearization allows us to define B(t) as a vector whose elements comprise the set of all bacterial subpopulations (B, B, …, B, …, B, B) at a given time t. We can then describe the temporal dynamics of the entire bacterial population as a system of linear differential equations:

In the equation above, denotes the matrix of coefficients describing the system of equations for the vector B(t). The values for the coefficients in depend on the concentration C of drug, on the drug’s binding kinetics, and on the growth and death rate functions G[i] and D[i].

Eq.(3) represents an initial value problem. This system of linear differential equations with a constant coefficient matrix has a unique solution given by

where the vector denotes the initial composition of bacterial subpopulations at t = 0. The solution can also be written as a linear superposition of a product of complex exponentials (with arguments determined by eigenvalues) and polynomials (whose degree is determined by the geometric multiplicity of these eigenvalues and whose coefficients are uniquely determined by the initial conditions). In practice, B(t) describes a family of exponential growth and decay curves that represent the replication and death of all N + 1 bacterial subpopulations over time (Fig. 1B). We solve for B(t) numerically by calculating the matrix exponential of using a scaling and squaring algorithm implemented in MATLAB (MathWorks, Newton, MA) [88].

Calculation of minimum inhibitory concentration

We define the MIC as the concentration C of an antibiotic such that any concentration of drug at or above C is guaranteed to cause the eventual extinction of the bacterial population. This occurs precisely when one eigenvalue of matrix (from Eq. (3)) is zero and all other eigenvalues have a negative real component. For this calculation, we assume that C is constant in time. We express the MIC as

With this formulation, finding the MIC amounts to finding the value of C such that the greatest real component of the eigenvalues of is zero. Deriving the expression for the MIC in the simplest case of the model, when N = 1, serves to illustrate this approach. For the purposes of this derivation, we consider a drug that elicits both a bactericidal and a bacteriostatic effect, so G[i = 1] = 0 and D[i = 1] = D. However, the approach for finding the MIC is identical for any mechanism of drug action. The matrix describing all bacterial subpopulations (B and B) in this simple case is

We wish to find the concentration C of antibiotic that yields negative real components of all but one eigenvalues λ of matrix . For the 2-by-2 matrix given by Eq. (6), the characteristic polynomial is given by λ2 - tr()λ + det(), and the Routh-Hurwitz stability criterion needed to satisfy the negative value constraints on λ is tr() ≤ 0 and det() ≥ 0. For the matrix described in Eq. (6), these expressions correspond to

and

Solving for the concentration C in both of these cases yields.

in the case of Eq. (7) and.

in the case of Eq. (8). We expect the value of the death rate at saturating drug concentrations (D) to be nonzero, positive, and larger than G. Therefore, Eq. (9) is guaranteed to be satisfied if Eq. (10) is also satisfied. We thus find the expression for the MIC to be

From this expression, we can infer the following proportionalities for the value of the MIC relative to the values of other model parameters:

Polynomial expressions for the MIC, as shown in Eq. (11), become exceedingly complex beyond N = 3. However, we conjecture (although we have not been able to prove) that the structure of the linear system shown in Eq. (3) guarantees the existence of the MIC for any N. For larger values of N, we leverage numerical schemes to calculate the eigenvalues of matrix . We use MATLAB’s eig() function, which calculates eigenvalues using the QZ algorithm [89].

Model calibration via simulated annealing

Numerical values for the model parameters N, D, µ, and µ were obtained from the literature (Table 1). For ciprofloxacin, the values for G and D were obtained by fitting experimental kill curves at drug concentrations of zero and 2.19 µg/ml, respectively, to exponential functions (Supplementary Fig. S6). For ampicillin, the values for G and D were obtained by fitting experimental kill curves at drug concentrations of zero and 256 µg/ml, respectively, to exponential functions. We leveraged an adaptive simulated annealing algorithm coupled with local gradient descent to obtain the remaining parameters (k, k, α, α, γ, and γ). Detailed descriptions of the adaptive simulated annealing algorithm are available elsewhere [48], [90]; in brief, simulated annealing is a global optimization algorithm capable of escaping local minima. It is therefore well suited to applications involving the optimization of many parameters. Adaptive simulated annealing is a variant on the classical simulated annealing algorithm that probes global parameter space with greater efficiency by accounting for each parameter’s magnitude when formulating a new parameter set at every iteration of the algorithm. We used adaptive simulated annealing to minimize the difference between experimental time-kill curves and model simulations of bacterial populations challenged to the same antibiotic doses. The difference between experimental observation and simulation is expressed through the objective function, whose value ψ the algorithm seeks to minimize:

denotes an m-by-n matrix of experimentally-measured population sizes at m drug concentrations and n timepoints, denotes simulated population sizes at the same drug concentrations and timepoints, and denotes an m-by-n weighting matrix (for our application, simply a matrix of ones). is a function of the parameters being optimized (that is,  = f(k, k, α, α, γ, γ)).

Coupling the adaptive simulated annealing optimization with a local gradient descent assures that our calibration procedure always converges on a local minimum. We used an exponential cooling schedule for the simulated annealing algorithm, which allows the optimization to run ergodically [90]. That is, repeating the optimization many times from random initial starting conditions in parallel yields roughly the same results as running the optimization once for a very long time. This allowed us to parallelize the optimization procedure by running the algorithm repeatedly across several cores of a computer and to characterize the distributions of parameter values obtained from these calibrations (Supplementary Fig. S3). After performing 249 independent model calibrations, we selected the parameter set with the lowest objective function value to use in subsequent simulations. The parameter values for this set are shown in Table 1. Parameter sets for all model optimizations performed are available in Supporting Data File S3.

Simulating the emergence of secondary mutations

We assumed that drug-resistant bacterial strains with secondary mutations that compensate for fitness costs and/or that further increase MIC emerge from preexisting drug-resistant subpopulations present in the initial population at the start of treatment (Supplementary Fig. S10). The size of the drug-resistant subpopulation in the absence of antibiotic (B) is given by the mutation-selection balance, which expresses the equilibrium at which the rate of emergence of drug resistance alleles by spontaneous mutation equals the rate of elimination of those alleles due to competitive fitness costs [91]:

Here, µ denotes the mutation rate for drug resistance emergence per unit time.

In order to quantify the probability of secondary mutation emergence from this drug-resistant subpopulation, we adapted a formulation that Lipsitch and Levin developed to study the evolution of drug-resistant bacterial strains during antibiotic treatment [78]. We assumed that secondary mutations emerge exclusively due to errors in DNA replication during bacterial growth. The expected number of resistant cells with secondary mutations that emerge from a bacterial population with i inactivated drug-target complexes (E(M)) is proportional to the total number of replications that the subpopulation undergoes before extinction and the rate of secondary mutation emergence:

In this equation, µ denotes the secondary mutation rate, G represents the growth rate of a resistant strain with exactly i inactivated drug-target complexes, B(t) describes the population dynamics of the ith drug-resistant bacterial subpopulation, and t describes the amount of time elapsed from treatment onset until the bacterial subpopulation is eliminated (B = 1 when t = t). The total number E(M) of resistant mutants with secondary mutations that we expect to observe over the course of treatment is thus the sum of Eq. (15) over all values of i, and the probability P that a compensated resistant mutant will emerge over the course of treatment follows from the Poisson assumption that secondary mutations arise stochastically and independently of other mutations:

The summation term in Eq. (16) describes the total number of resistant strains with secondary mutations expected to emerge before extinction. This equation thus quantifies the Poisson probability that at least one resistant strain with a secondary mutation will emerge over the course of treatment.

Code and data

We wrote all code in MATLAB. All of the code used to implement and solve our model, to analyze experimental data, and to generate simulation data shown in all figures is available as a software package in Supplementary File S1. Experimental data represented in Fig. 2A & 2B and in Supplementary Fig. S4 are available within Supporting Data Files S1, S2 & S4, respectively, and the parameter values for all iterations of model optimization are available in Supporting Data File S3.

Funding

This work was funded by Bill and Melinda Gates Foundation Grant OPP1111658 (to T.C. & P. AzW.); Research Council of Norway (NFR) Grant 262686 (to P.AzW.); and the Yale School of Public Health Ralph Skolnik Summer Internship Fund (to C.H.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

CRediT authorship contribution statement

Colin Hemez: Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. Fabrizio Clarelli: Formal analysis, Investigation, Writing – review & editing. Adam C. Palmer: Investigation, Methodology, Resources, Writing – review & editing. Christina Bleis: Investigation. Sören Abel: Methodology, Project administration. Leonid Chindelevitch: Formal analysis, Methodology, Validation, Writing – review & editing. Theodore Cohen: Funding acquisition, Project administration, Resources, Writing – review & editing. Pia Abel zur Wiesch: Conceptualization, Formal analysis, Funding acquisition, Methodology, Project administration, Resources, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.