Drug-target binding quantitatively predicts optimal antibiotic dose levels in quinolones

Antibiotic resistance is rising and we urgently need to gain a better quantitative understanding of how antibiotics act, which in turn would also speed up the development of new antibiotics. Here, we describe a computational model (COMBAT-COmputational Model of Bacterial Antibiotic Target-binding) that can quantitatively predict antibiotic dose-response relationships. Our goal is dual: We address a fundamental biological question and investigate how drug-target binding shapes antibiotic action. We also create a tool that can predict antibiotic efficacy a priori. COMBAT requires measurable biochemical parameters of drug-target interaction and can be directly fitted to time-kill curves. As a proof-of-concept, we first investigate the utility of COMBAT with antibiotics belonging to the widely used quinolone class. COMBAT can predict antibiotic efficacy in clinical isolates for quinolones from drug affinity (R2>0.9). To further challenge our approach, we also do the reverse: estimate the magnitude of changes in drug-target binding based on antibiotic dose-response curves. We overexpress target molecules to infer changes in antibiotic-target binding from changes in antimicrobial efficacy of ciprofloxacin with 92–94% accuracy. To test the generality of our approach, we use the beta-lactam ampicillin to predict target molecule occupancy at MIC from antimicrobial action with 90% accuracy. Finally, we apply COMBAT to predict antibiotic concentrations that can select for resistance due to novel resistance mutations. Using ciprofloxacin and ampicillin as well defined test cases, our work demonstrates that drug-target binding is a major predictor of bacterial responses to antibiotics. This is surprising because antibiotic action involves many additional effects downstream of drug-target binding. In addition, COMBAT provides a framework to inform optimal antibiotic dose levels that maximize efficacy and minimize the rise of resistant mutants.


Introduction
The rise of antibiotic resistance represents an urgent public health threat. In order to effectively combat the spread of antibiotic resistance, we must optimize the use of existing drugs and develop new drugs that are effective against drug-resistant strains. Accordingly, methods to improve antibiotic dose levels to i) maximize efficacy against susceptible strains and ii) minimize resistance evolution play a key role in our defense against antibiotic resistant pathogens.
It is noteworthy that dosing strategies for treatment of susceptible strains (e.g., dosing level [1], dosing frequency [2], and treatment duration [3][4][5]) have recently been substantially improved, even for antibiotic treatments that have been standard of care for decades. This suggests that there likely remains significant room for optimization in our antibiotic treatment regimens. It also highlights the difficulty in identifying optimal dosing levels for new antibiotics. Indeed, optimizing dosing is one of the biggest challenges in drug development. Traditionally, antibiotic efficacy was mainly described by a single value, the minimal inhibitory concentration (MIC). While correlations between treatment success and MIC have been demonstrated there is limited predictive power [6,7]. When susceptibility is assessed by MIC, not all patients infected with "susceptible" bacteria are successfully treated with antibiotics. Additionally, a large majority of patients with a "resistant" infection can be successfully treated with antibiotics even when the underlying infection is serious and untreated patients are unlikely to recover [8]. This was e.g. shown for patients with complicated intraabdominal infections [7]. Reasons for this mismatch may include that the MIC only gives the minimal concentration to suppress bacterial growth and contains no information on antibiotic efficacy above or below MIC [9]. This makes the MIC ill-suited to describe efficacy of the strongly fluctuating antibiotic concentrations in patients. This has led to an increase in more sophisticated dose-response measurements where bacteria are exposed to multiple antibiotic concentrations and the kill rate is assessed at each concentration individually (pharmacodynamic profiles). However, these approaches require orders of magnitude more experimental effort than simple MIC measurements because they involve a multitude of antibiotic concentrations and time points. This process is too time-consuming when testing new drug candidates.
It is even more challenging to optimize dose levels to minimize the emergence of antibiotic resistance, both for existing and novel antibiotics. Typically, not only the MIC changes when a strain acquires resistance, also other properties such as the steepness of the dose-response curve and the maximal kill rate at very high concentrations change [10]. Predicting the changes in the dose-response curve is therefore not trivial. Thus, a full pharmacodynamic profile should be assessed for each potential resistant strain. To this end, resistant strains must be isolated and due to the amount of different resistance mechanisms, a good saturation of the mutational target must be achieved. This requires substantial and lengthy evolutionary experiments. In addition, there remains substantial debate about which dosing strategies best prevent the emergence of resistance during treatment [11][12][13]. In this context, a useful concept that links antibiotic concentrations with resistance evolution is the resistance selection window (mutant selection window) that ranges from the lowest concentration at which the resistant strain grows faster than the wild-type, usually well below the wild-type MIC, to the MIC of the resistant strain [14][15][16]. Antibiotic concentrations above the resistance selection window safeguard against de novo resistance emergence. Antibiotic concentrations below the resistance selection window do not kill the susceptible strain, but also do not favor the resistant strain and therefore do not promote emergence of resistance. To limit resistance, it is therefore important to identify the resistance selection window and optimize dosing accordingly. However, this again requires obtaining a full pharmacodynamic profile of a majority of the expected resistant mutants and is therefore not feasible as a standard assessment in drug development.
The next challenge to successfully designing antibiotic treatment arises when the experimental information is integrated into mathematical pharmacodynamic models that then predict efficacy under realistic, fluctuating concentrations in patients. Pharmacodynamic models from 1910 (E max or Hill-models) [17] are still widely used despite assuming instantaneous equilibria of antibiotic-target binding and therefore being often inaccurate when antibiotic concentrations fluctuate. Recently described models that relax these assumptions have been useful in gaining a better qualitative understanding of realistic dosing and complicated drug effects, such as post-antibiotic effects, inoculum effects, and bacterial persistence [18][19][20][21]. However, to speed the development of new antibiotics or to inform practices which minimize resistance, we require quantitative predictions for antibiotics or resistant bacterial strains that do not exist yet. Models which permit quantitative predictions of changes in drug efficacy as a function of modification of antibiotic molecules (i.e. new drugs) or novel resistance mutations would be invaluable. Such tools would advance our general mechanistic understanding of antibiotic action, could guide dosing trials of new drugs, and suggest better dosing of existing drugs.
In this report, we describe a mechanistic computational modeling framework (COMBAT-COmputational Model of Bacterial Antibiotic Target-binding) that allows us to predict full pharmacodynamic profiles based solely on accessible biochemical parameters describing drugtarget interaction. These parameters can be determined early in drug development. We use this framework to investigate how changes in drug target binding, either due to improvements in existing antibiotics or due to resistance mutations in bacteria, affect antibiotic efficacy. We first show that COMBAT accurately predicts bacterial susceptibility as a function of drug-target binding and, conversely, allows inference of these biochemical parameters on the basis of observed patterns of bacterial growth suppression or killing. We then use COMBAT to predict the susceptibility of newly arising resistant variants based on the molecular mechanism of resistance and determine the resistance selection window.

Quinolone target affinities correlate with antibiotic efficacy
To investigate how biochemical changes in antibiotic action modifies bacterial susceptibility, we explored how the affinity of antibiotics to their target affects the MIC. We compared the MICs of quinolones, an antibiotic class in which individual antibiotics have a wide range of affinities to one of their targets, gyrase (K D~1 0 −4 -10 −7 M) but are of similar molecular sizes and have a similar mode of action [22]. This choice allowed us to isolate the effects of differences in drug-target affinity on the MIC.
We obtained binding affinities of quinolones to their gyrase target in Escherichia coli from previous studies [23][24][25][26][27]. We then retrieved MIC data for several quinolones from clinical Enterobacteriaceae isolates collected before 1990 [28], i.e., before the widespread emergence of quinolone resistance [22]. We assume that quinolone affinities obtained from clinical Enterobacteriaceae isolates collected before the emergence of resistance correspond to those measured in wild-type E. coli.
To make qualitative predictions of MICs, we employed a simplified model based on the assumptions that i) drug-target binding occurs much more quickly than bacterial replication, ii) the antibiotic concentration remains constant and iii) that during the 18 hours of an MIC assay, the concentration gradient of the drug inside and outside the cell has equilibrated. Under these assumptions, the MIC can be expressed as where K D represents the affinity constant and f c the fraction of the target bound at the MIC [29]. Accordingly, this model predicts that the MIC is linearly correlated with K D . Fig 1 shows the correlations between drug-target affinities and MICs for seven quinolones and clinical isolates of 11 different Enterobacteriaceae species. We observed a significant (p < 0.018) linear correlation between MIC and K D in all species, confirming the qualitative model prediction.

A quantitative model to predict antibiotic efficacy
While it was encouraging that our model can qualitatively predict MIC changes, our aim was to quantitatively predict antibiotic treatment performance. The simplified model assumes that the binding kinetics are much faster than bacterial replication, which may not be true in all cases. To expand the generalizability of the model, we extended the modeling framework to allow that bacterial replication may occur in a similar time frame as drug-target binding events.
The full model (COMBAT-COmputational Model of Bacterial Antibiotic Target-binding) describes the binding and unbinding of antibiotics to their targets and predicts how such binding dynamics affects bacterial replication and death (Fig 2A). In previous work linking drugtarget binding kinetics with bacterial replication [21], we described a population of bacteria with θ target molecules per cell with a system of θ + 1 (bacteria with 0, 1, . . ., θ bound target molecules) ordinary differential equations (ODEs). This system increases in complexity with the number of target molecules and makes fitting the model to data computationally too demanding for most settings. To simplify this prior approach, we developed new mathematical models based on partial differential equations (PDEs), where a single equation describes all bacteria simultaneously. The sum of bacteria within all target occupancy states over time can be described by a time kill curve (Fig 2B), during which the bacterial population is characterized by the distribution of bacterial cells with different levels of target occupancies at each time-step ( Fig 2C). This curve can be visualized as a two-dimensional surface in a threedimensional coordinate system where the number of bacteria is represented on the z-axis, the percent of bacteria with the fraction of bound target molecules on the x-axis, and time on the y-axis (Fig 2D).
Antibiotic action is described by rates of binding (k f ) and unbinding (k r ) to bacterial target molecules (Fig 2A and 2E). The binding of an antibiotic to a target results in the formation of an antibiotic-target molecule complex x, where x ranges between 0 and θ.
COMBAT consists of two mass balance equations: Eq 2 describing bacterial numbers as a function of bound targets and time and Eq 3 describing antibiotic concentration as a function of time (see Methods).  [24,[26][27][28], and the y-axes show the MICs, both in mol/L. The adjusted R 2 and p-value of each correlation are given. In cases where there was more than one K D value reported in the literature, we used the mean for this analysis. The tested MIC values are the median of several clinical isolates described previously [28]. https://doi.org/10.1371/journal.pcbi.1008106.g001

PLOS COMPUTATIONAL BIOLOGY
Drug-target binding quantitatively predicts optimal antibiotic dose levels in quinolones and v r can be seen as a generalized velocity v ¼ dx dt . Eq 4 (part of the replication term in Eq 2) describes how daughter cells inherit bound target molecules from the mother cell during replication: Eq 5 (part of the replication term in Eq 2) is a logistic growth model describing reduced bacterial replication as the carrying capacity is approached:

Model fit to ciprofloxacin time-kill data
We used the quinolone ciprofloxacin to quantitatively fit bacterial time-kill curves, since this is a commonly used antibiotic for which binding parameters have been directly measured. S1 Table gives an overview of the known parameters used for fitting; S2 Table gives the parameters resulting from our fit.
The functional relationship between the levels of bacterial replication and death on the fraction of bound target molecules is extremely hard to obtain experimentally. We therefore treated the relationships between the fraction of bound target and bacterial replication and death as free parameters in our model fitting. Ciprofloxacin is considered to have both bacteriostatic and bactericidal action (mixed action) [30,31], and we fitted functions for a monotonically decreasing replication and a monotonically increasing killing with each successively bound target molecule (see Methods & S1 Fig).
Overall, we found that COMBAT could fit the time-kill curves well (R 2 = 0.93, Fig 3A). Fig  3B shows the predicted bacterial replication r(x) and death as a function of target occupancy δ (x) based on the fit obtained in Fig 3A. After model calibration, we simulated bacterial replication during exposure to different antibiotic concentrations for 18 h. For this simulation, positive values indicate an increase in the number of bacteria, and negative values indicate a decrease in the number of bacteria. We estimated a MIC of 0.0139 mg/L (Fig 3C), a value that is within the range of MIC determinations for wt E. coli (0.01 mg/L, 0.015 mg/L, 0.017 mg/L and 0.023 mg/L [15,[32][33][34]). unbinding events of the antibiotic to its target molecule in the cell.k f is the adjusted forward reaction rate, k r is the reverse reaction rate, A is the concentration of antibiotics inside the bacterium, x is the number of bound targets, θ is the number of targets and B x is the number of bacteria with x bound targets. b, Modeled sample time-kill curve, in which the sum of bacteria in all binding states (i.e., the entire population of living bacteria) is followed over time after exposure to antibiotics. The vertical dotted lines indicate the time points depicted in (c); 1 min (grey), 14 min (yellow), and 80 min (purple). c, The percentage of bound antibiotic targets in the bacterial population at indicated time points. d, Illustration of how the partial differential equation describes the bacterial population as a surface in a three-dimensional coordinate system, the dimensions of which represent percent bound target (x-axis), time (y-axis), and number of bacteria (z-axis). The three time points shown in (c) represent two-dimensional cross-sections at different points of the y-axis. e, Overview of used parameters and functions. https://doi.org/10.1371/journal.pcbi.1008106.g002

Accurate prediction of target overexpression from time-kill data
Having shown that COMBAT can quantitatively fit experimental data on antibiotic action within biologically plausible parameters, we continued to test the predictive ability of the model. Given our hypothesis that modifications in antibiotic-target interactions lead to predictable changes in bacterial susceptibility, we experimentally induced changes in the antibiotic-target interaction of ciprofloxacin in E. coli. We then quantified these biochemical changes by fitting COMBAT to corresponding time-kill curves and compared them to the experimental results. Ciprofloxacin acts on gyrase A 2 B 2 tetramers [22]. We used an E. coli strain for which both gyrase A and gyrase B are under the control of a single inducible

PLOS COMPUTATIONAL BIOLOGY
Drug-target binding quantitatively predicts optimal antibiotic dose levels in quinolones promoter (P lacZ ), such that the amount of gyrase A 2 B 2 tetramer can be experimentally manipulated [35]. We measured net growth rates for this strain at different ciprofloxacin concentrations in the presence of 10 μM isopropyl β-D-1-thiogalactopyranoside (IPTG; mild overexpression) and 100 μM IPTG (strong overexpression) and compared it to the wild-type in the absence of the inducer (Fig 4A).
Like previously reported, we find that increasing gyrase content makes E. coli more susceptible to ciprofloxacin [35]. We fitted net growth rates allowing the target molecule content, i.e. gyrase A 2 B 2 , to vary. We assumed that the only change between the different conditions was the amount of target. We further assumed that the relationship between bound target and bacterial replication or death did not differ between the control strain containing a mock plasmid (no IPTG) and the experiments with overexpression ( Fig 4B, between 0% and 100%). Finally, we assumed that the maximal kill rate at very high antibiotic concentrations was accurately measured in our experiments and forced the function describing bacterial death through the measured value when all target molecules are bound. We found the best fit for a 1.31x increase in GyrA 2 B 2 target molecule content for bacteria grown in the presence of 10 μM IPTG and a 2.02x increase in GyrA 2 B 2 target molecule content for those grown in the presence of 100 μM IPTG.
We subsequently tested these predictions experimentally by analyzing Gyrase A and B content by western blot Figs 4C and S2). Using realistic association and dissociation rates for biological complexes [36], we predicted a range of functional tetramers based on the relative amount of Gyrase A and B proteins ( Fig 4D). S3 Table details the individual measurements, and the procedure to estimate tetramers is provided in the methods section. We found that the observed overexpression was very close to our theoretical prediction, with 1.

Accurate prediction of target occupancy at MIC from time-kill data
Next, we tested whether COMBAT can be applied to the action of the beta-lactam ampicillin, a very different antibiotic with a distinct mode of action from quinolones. Using published pharmacodynamic data of E. coli exposed to ampicillin [34] also allowed us to compare COMBAT predictions to established pharmacodynamic approaches. Most of the biochemical parameters for ampicillin binding to its target, penicillin-binding proteins (PBPs), have been determined experimentally (S1 Table). Ampicillin is believed to act as a bactericidal drug [37], and this mode of action is supported by findings from single-cell microscopy [29]. We therefore assume that ampicillin binding does not affect bacterial replication. In order to model the consumption of beta-lactams at target inhibition and eventual target recovery, we made small adjustments to Eq 13 (see Methods, description of beta-lactam action).
We fitted COMBAT to published time-kill curves of E. coli exposed to ampicillin ( Fig 5A). Again, COMBAT provides a good fit to the experimental data between 0 min and 40-60 min. After that time, observed bacterial killing showed a characteristic slowdown at high ampicillin concentrations which is often attributed to persistence [21] (Fig 5A). For the sake of simplicity, we chose to omit bacterial population heterogeneity in this work and therefore cannot describe persistence, even though COMBAT can be adapted to capture this phenomenon [21]. Because ampicillin acts in an entirely bactericidal manner, we assume a constant replication rate (see  Table). Fig 5C shows the predicted net growth rate over a range of drug concentrations. We estimated a MIC of 2.6 mg/L. This MIC is based on the Clinical & Laboratory Standards Institute definition of the MIC determined at 18 h. The original source of the MIC, which was based on experimental data and a pharmacodynamic model [34] determined an MIC of 3.4 mg/L at 1 h. If we change our prediction to 1 h, our estimated MIC is 3.32 mg/L, which is within 2.5% of the reported value [34].
Having established that COMBAT can also adequately capture the pharmacodynamics of ampicillin, we next tested whether we can estimate experimentally determined target occupancy at the MIC. Our estimated mean occupancy considering both living and dead bacteria is 89% (Fig 5B), a value within previously reported experimental estimates from Staphylococcus aureus (84-99%) [38].

Sensitivity of antibiotic efficacy to parameters of drug-target binding
It is possible to vary all parameters in COMBAT and explore their effect. We used this to test how hypothetical chemical changes to ampicillin or ciprofloxacin would affect antibiotic efficacy (S3-S11 Fig). These changes could reflect either bacterial resistance mutations or modifications of the antibiotics themselves. We predict that changes in drug-target affinity, K D , have more profound effects than changes in target molecule content, bacterial reaction to increasingly bound target (i.e. δ(x) and r(x)), or changes in target molecule content. We also predict that the individual binding rates k r and k f , and not just the ratio of these terms, the K D , are important factors in efficiency. The faster a drug binds, the more efficient we predicted it will be. One intuitive explanation for the observation that k f drives efficacy is that a slow binding fails to rapidly interfere with bacterial replication, which may allow for the production of additional target molecules and thereby reduce the ratio of free antibiotic to target molecules.

Forecasting the resistance selection window
Finally, we illustrate how COMBAT can be used to explore how the molecular mechanisms of resistance mutations affect antibiotic concentrations at which resistance can emerge, i.e., the resistance selection window. We compared predicted net growth rates as a function of ciprofloxacin concentrations for a wild-type strain and an archetypal resistant strain. For this analysis, we assumed that the resistant strain has a 100x slower drug-target binding rate (i.e.~100x  [34]. The points represent experimental data, and the lines represent the fit of the model. Each color indicates a single ampicillin concentration, as described in the legend. b, Replication (blue) and death (red) rates as a function of the number of bound targets predicted by the model fit in (a). The black line indicates the predicted distribution of target occupancies in a bacterial population (both living and dead cells) exposed to ampicillin at the MIC for 18 h. c, The net growth rate, as determined by the slope of a line connecting the initial bacterial density and the bacterial density at 18 h on a logarithmic scale predicted from the model fit in (a), is shown as function of the drug concentration (blue). The dotted horizontal line indicates zero net growth, and the intersection with the blue line predicts the MIC (2.6 mg/mL). https://doi.org/10.1371/journal.pcbi.1008106.g005

PLOS COMPUTATIONAL BIOLOGY
Drug-target binding quantitatively predicts optimal antibiotic dose levels in quinolones increased MIC, realistic for novel point mutations [39]) and that the maximum replication rate of the resistant strain is 85% of the wild type strain [40]. We then predicted the antibiotic concentrations at which resistance would be selected. Interestingly, when comparing COM-BAT to previous pharmacodynamics models (Fig 5), we observed that estimates of replication rates depend on the selected time frame (Fig 6A). When the timeframe for MIC determination is set to 18 h as defined by CLSI [41], the "competitive resistance selection window", i.e., the concentration range below the MIC of both strains where the resistant strain is fitter than the wild type, ranges from 0.002 mg/L to 0.014 mg/L for ciprofloxacin ( Fig 6A) and 1 mg/L to 2.6 mg/L for ampicillin (S12 Fig), respectively. This corresponds well with previous observations that ciprofloxacin resistance is selected for well below MIC [15]. However, when measuring after 15 min or 45 min, the results are substantially different. The reason for this is illustrated in Fig 6B. COMBAT reproduces non-linear time kill curves where bacterial replication continues until sufficient target is bound to result in a negative net growth rate. This compares well with experimental data around the MIC in Figs 3A and 5A. In Fig 6B,  Thus, we estimate that sustained levels of 1.27-7mg/L would safeguard against resistance. While ciprofloxacin plasma concentrations typically reach concentrations of 2mg/L after oral uptake and 6mg/L after intravenous administration [42], levels of 2.6 mg/L and above were shown to be chondrotoxic in young animals [43] and concentrations of 40 mg/L are toxic to mitochondria [44]. Clearly, toxicity and risk of resistance must be carefully weighed when deciding on dosing.

Discussion
Optimizing dosing levels of antibiotics is important for maximizing drug efficacy against wildtype strains as well as for minimizing the rise of resistant mutants. Antibiotic efficacy is traditionally described by a single value, the minimal inhibitory concentration (MIC), which has limited predictive power [7,8]. In more sophisticated dose-response measurements, bacteria are exposed to multiple antibiotic concentrations and the kill rate is assessed at each concentration individually in dose-response curves (pharmacodynamic profiles). However, this approach requires substantial experimental effort and is too time-consuming when testing large libraries of new drug candidates. Limited predictive power of standard measures of pharmacodynamics is not only a problem for antibiotic development, drug attrition in general is mainly due to insufficient predictions of pharmacodynamics rather than pharmacokinetics [45].
Because of the experimental effort, pharmacodynamic profiles for either novel drug candidates or novel resistant strains are often not obtained. Thus, we need a transferrable framework that allows quantitative predictions based on parameters that can be determined a priori. Recent studies have reported methods to predict MICs from whole genome sequencing data [46,47]. However, these methods require transfer of prior knowledge on how the resistance mutations affect MICs in other organisms. There are no methods that could predict a priori how chemical changes to an antibiotic structure or novel resistance mutations affect bacterial growth at a given antibiotic concentration.
Here, we accurately predict antibiotic action on the basis of accessible biochemical parameters of drug-target interaction. Our computational model, COMBAT provides a framework to predict the efficacy of compounds based on drug-target affinity, target number, and target occupancy. These parameters may change both when improving antibiotic lead structures as well as when bacteria evolve resistance. Importantly, they can be measured early in drug development and may even be a by-product of target-based drug discovery [48]. When these data are available, COMBAT makes only one assumption: that the rate of bacterial replication decreases and/or the rate of killing increases with successive target binding. While fitting, we allow this relationship to be gradual or abrupt and select the best fit. This means we do not model specific molecular mechanisms down-stream of drug-target binding, but their effects are subsumed in the functions that connect the kinetics of drug-target binding to bacterial replication and death.
In previous work, for example on antipsychotics [18], antivirals [19] and antibiotics [20,21], models of drug-target binding kinetics have been used to improve our qualitative understanding of pharmacodynamics. Our study substantially advances this work by making quantitative predictions across antibiotics and bacterial strains when measurable biochemical characteristics change with extremely high accuracy. This is possible because COMBAT employs an efficient and versatile mathematical approach, based on partial differential equations, that makes it computationally feasible to fit the model to a large range of data. Importantly, we are not only able to predict antibiotic action from biochemical parameters, but can also vice versa use COMBAT to accurately predict biochemical changes from observed patterns of antibiotic action. We have confirmed the excellent predictive power of COMBAT with clinical data as well as experiments with antibiotics with very different mechanisms of action. The high predictive power makes it possible to use modeling to guide dosing. This gives us confidence that biochemical parameters are major determinants of antibiotic action in bacteria and that COMBAT helps to make rational decisions about antibiotic dosing.
In drug development, our mechanistic modeling approach provides insight into which chemical characteristics of drugs may be useful targets for modification. For example, our sensitivity analyses indicate that antibiotics with a similar affinity but faster binding inactivate bacteria more quickly and therefore prevent replication and production of more target molecules, which would change the ratio of antibiotic to target. Furthermore, because e.g. antibiotic binding and unbinding rates can be determined early in the drug development process, such insight can help the transition to preclinical and clinical dosing trials. This may contribute to reducing bottlenecks between these phases of drug development and thereby save money and time.
Avoiding antibiotic concentrations that select for resistance is challenging for two reasons. First, the differences in the pharmacodynamic curves of wild-type and resistant strains are not trivial. Resistance can affect not only the MIC, but also the maximal kill rate at high drug concentrations and the steepness of the dose-response curve [10]. Therefore, one would need to record full pharmacodynamic profiles rather than just MICs to assess the mutation selection window for resistant mutants. Second, this process would have to be repeated for all (or at least a representative set of) potential emerging resistant mutants. This makes it extremely time-and resource-consuming to safeguard against resistance by determining the resistance selection window.
COMBAT offers insight into determinants of the resistance selection window and builds transferrable knowledge that allows estimating useful dose ranges. In concordance with a recent meta-analysis of experimental data [49], our sensitivity analyses predict that changes in drug target binding and unbinding have a greater impact on susceptibility than changes in target molecule content or down-stream processes. Thus, a more comprehensive characterization of the binding parameters of spontaneous resistant mutants would allow an overview of the maximal biologically plausible levels of resistance that can arise with one mutation. Dosing above this level should then safeguard against resistance. This is especially useful for compounds for which it is difficult to saturate the mutational target for resistance, or for safeguarding against resistance to newly introduced antibiotics for which we do not yet have a good overview of resistance conferring mutations.
Good quantitative estimates on the dose-response relationship of new drugs would also help defining the therapeutic window, i.e. the range of drug concentrations at which the drug is effective but not yet toxic. For example in ciprofloxacin, the doses found necessary to prevent resistance after marketing were found to be toxic [50], and an early assessment of doses that might become necessary after resistance is wide-spread might preserve antibiotic utility. If toxicity, solubility or other constraints do not allow dosing above the MIC of expected resistant strains, COMBAT can also predict the concentration range at which resistance is less strongly selected. This could guide decisions on treating with low versus high doses, which is currently controversially debated [11,12]. COMBAT therefore offers new promise to reduce the failure rates of candidate compounds late in the drug development process when resistance is observed in patients and substantial resources have been invested.
Our quantitative work can help to identify optimal dosing strategies at constant antibiotic concentrations for homogeneous bacterial populations. These measures are commonly used to assess antibiotic efficacy. In addition, previous work has demonstrated that drug-target binding models outperform traditional pharmacodynamic models for the fluctuating concentrations that actually occur in patients [29,51]. To illustrate this, we coupled a pharmacokinetic model describing different modes ampicillin administration in patients and predict the pathogen load in infected tissues based on the realistic, fluctuating antibiotic concentrations (S1 Text, S13 and S14 Figs and S5 Table). They can also explain complicated phenomena such as biphasic kill curves, the post-antibiotic effect, or the inoculum effect [20,21,52] that often complicate the clinical phase of drug development. COMBAT has similar characteristics that allow capturing these complex phenomena. Therefore, employing COMBAT may be useful for guiding drug development to maximize antibiotic efficacy and minimize de novo resistance evolution.

Mathematical model
COMBAT incorporates the binding and unbinding of antibiotics to their targets and describes how target binding affects bacterial replication and death. This work extends the model developed in [21]. COMBAT consists of a system of two mass balance equations: one PDE for bacteria (describing replication and death as a function of both time and target binding) and one ODE for antibiotic molecules (describing the concentrations as function of time).
In the most basic version of COMBAT, we ignored differences between extracellular and intracellular antibiotic concentrations and only followed the total antibiotic concentration A, assuming that the time needed for drug molecules to enter bacterial cells is negligible. We model ciprofloxacin (to which there is a limited diffusion barrier [53]) and ampicillin (where the target is not in the cytosol, even though the external membrane in gram negatives has to be crossed to reach PBPs). We therefore believe that this assumption is justified in wild-type E. coli. This basic version of COMBAT is therefore more accurate for describing antibiotic action where the diffusion barrier to the target is weak.

Binding kinetics
We describe the action of antibiotics as a binding and unbinding process to bacterial target molecules [21]. For simplicity, we assume a constant number of available target molecules θ. The binding process is defined by the formula A + T Ð x, where the intracellular antibiotic molecules A react with target molecules T at a rate k f and form an antibiotic-target molecule complex x, where values for x range between 0 and θ. If the reaction is reversible, the complex dissociates with a rate k r .
In [21], the association and dissociation terms are described by the following terms dB i ðtÞ dt ¼k f AðtÞððy À i þ 1ÞB iÀ 1 ðtÞ À ðy À iÞB i ðtÞÞ zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {

Association term
À k r ðiB i ðtÞ À ði þ 1ÞB iþ1 ðtÞÞ zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {

Dissociation term
; i�½0; y� ð6Þ wherek f ¼ k f V tot n A , k f is the association rate, V tot is the volume in which the experiment is performed, n A is Avogadro's number, k r is the dissociation rate, B i is the number of bacteria with i bound targets, and θ is the total number of targets.
This approach requires the use of a large number of ordinary differential equations, (0 + 1) for the bacterial population and one for the antibiotic concentration. To generalize this approach, we assume that the variable of bound targets is a real number x 2 R. Under this continuity assumption, we consider the bacterial cells as a function of x and the time t, thereby reducing the total number of equations to two.
Under the continuity approximation (x 2 R), we can rewrite the binding kinetics in the form @Bðx; tÞ @t ¼ @ @x �k f AðtÞðy À xÞBðx; tÞ � zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl { Association term À @ @x � k r xBðx; tÞ � zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {

Dissociation term ð7Þ
or simply @Bðx; tÞ where v f ¼k f A t ð Þ y À x ð Þ and v r = k r x can be considered as two velocities, i.e., the derivative of the bound targets with respect to the time dx dt . Replication rate. We assume that the replication rate of bacteria, r(x), is dependent on the number of bound target molecules x. The function r(x) is a monotonically decreasing function of x, such that fewer bacteria replicate as more target is bound. r(0) is the maximum replication rate, corresponding to the replication rate of bacteria in absence of antibiotics. Thus, r(x) describes the bacteriostatic action of the antibiotics, i.e., the effect of the antibiotic on bacterial replication.
Carrying capacity. Replication ceases as the total bacterial population approaches the carrying capacity K. At that point, the replication term of the equation is @Bðx; tÞ @t ¼ rðxÞBðx; tÞ where is the replication-limiting term due to the carrying capacity K, and 0 � F lim � 1.
Distribution of target molecules upon division. We assume that the total number of target molecules doubles at replication, such that each daughter cell has the same number as the mother cell. We also assume that the total number of drug-target complexes is preserved in the replication and that the distribution of x bound target molecules of the mother cell to its progeny is described by a hypergeometric sampling of n molecules from x bound and 2θ−x unbound molecules. Under the continuity assumption, we generalize the concept of hypergeometric distribution. Because the hypergeometric distribution is a function of combinations and because a combination is defined as function of factorials, we can use Γ functions in place of factorials and redefine a continuous hypergeometric distribution as a function of Γ functions. A Γ function is where z is a complex number. In this way, the distribution can be expressed as a probability density function of continuous variables. The amount of newborn bacteria is given by the term r(x)B(x, t)F lim (t). We assume that bound target molecules are distributed randomly between mother and daughter cells, with each of them inheriting 50% upon division on average. This means that twice the amount of newborn cells must be redistributed along x to account for the random distribution process. For example, if a mother cell with 4 bound targets divides, we have two daughter cells, each with a number of bound targets between 0 and 4 (their sum has to be 4), following the generalized hypergeometric distribution. For simplicity, we define S(x,t) to be a function related to the replication rate that depends on the number of bacteria with a number of bound target molecules ranging between x and θ, their specific replication rate r(x), and the fraction of their daughter cells expected to inherit x antibiotic-target complexes h(x,z): Death rate. The death rate function δ(x) depends on the number of bound target molecules. The function δ(x) is assumed to be a monotonically increasing function of x, where δ(θ) is the maximum death rate, when all targets in the bacteria have been bound by antibiotics. The shape of this function describes the bactericidal action of the antibiotic.
Bacteriostatic and bactericidal effects. We consider several potential functional forms of the relationship between the percentage of bound targets and replication and death rates, because the exact mechanisms how target occupancy affects bacteria is unknown (S1 Fig). We use a sigmoidal function that can cover cases ranging from a linear relationship to a step function. When the inflection point of a sigmoidal function is at 0% or 100% target occupancy, the relationship can also be described by an exponential function. We assume that replication in bactericidal and death in bacteriostatic drugs is independent of the amount of bound target. With sufficient experimental data, the replication rate r(x) and/or the death rate δ(x) can be obtained by fitting COMBAT to time-kill curves of bacterial populations after antibiotic exposure. The sigmoidal shape of r(x) and δ(x) can be written as: where x rth is the replication rate threshold, x dth is the death rate threshold, and both represent the point where the sigmoidal function reaches ½ of its maximum. γ r and γ d represent the shape parameters of the replication and death rate functions, respectively. These factors determine the steepness around the inflection point. When they are extreme, the relationship approaches a linear or a step function. Full equation describing bacterial population. Putting these components together, the full equation describing a bacterial population is: @Bðx; tÞ @t þ @ @x ðv f ðx; tÞBðx; tÞ À v r ðx; tÞBðx; tÞÞ zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {

Description of beta-lactam action.
Beta-lactams acetylate their target molecules (PBPs) and thereby inhibit cell wall synthesis. The acetylation of PBPs consumes beta-lactams. However, PBPs can recover through deacetylation. We modified the term of drug-target dissociation in the equation describing antibiotic concentrations (Eq 3), and set the unbinding rate k r = 0. To reflect the recovery of target molecules, we substituted the dissociation rate k r in the equation describing the bacterial population with the deacetylation rate k a , as described in [29].
Initial and boundary conditions. At t = 0, we assume that all bacteria have zero bound targets (x = 0), and the initial concentration of bacteria is B(x, 0) = 0, x > 0, and B(0,0) At the boundaries of the partial differential equation (x = 0, x = θ), we specify that the outgoing velocities are zero. For x = 0, i.e. no bound target molecules, the unbinding velocity v r (0, t) = 0, and in x = θ, i.e. all targets are bound, the binding velocity v f (θ, t)=0. When the replication term at x = 0 and the death term at x = θ are known, we can solve the partial differential equation with two ordinary differential equations at the boundaries. They are similar to the equations at x = 0 and at x = θ described by Abel zur Wiesch et al. [21], but taking into account that x is a continuous variable instead of a natural number.
Numerical schemes. To solve our system of differential equations, we used a firstorder upwind scheme. Specifically, we used the spatial approximation u f À ¼ u i ð ÞÀ uðiÀ 1Þ Dx for the binding term (v f > 0) and the spatial approximation for the unbinding term (v r < 0). For the time approximation of both the PDEs and the ODEs, we used the forward approximation DB Dt ¼ B nþ1 À B n Dt [54]. We also verified that the Courant-Friedrichs-Lewy condition is satisfied. For fitting the experimental data of bacteria exposed to ciprofloxacin and ampicillin, we used the particle swarm method ("particleswarm" function in Matlab, MathWorks software).
Time-kill curves. Overnight cultures of BW25113 or SoA3329 and SoA3330 were diluted 1:1000 in pre-warmed LB or LB-Cm and LB-Cm-IPTG, respectively, and grown with shaking to OD 600~0 .5. A 1:3 dilution series of ciprofloxacin was made and added to the cultures at indicated concentrations. Additional cultures without antibiotics and with a very high concentration of ciprofloxacin (2.187 mg/L) were used to determine the minimal and maximal kill rate, respectively. Samples were taken immediately prior to addition of the antibiotic and iñ 20 min intervals or after 45 min, respectively. Samples were washed once in phosphate buffered saline (PBS) before colony forming units (CFUs) were determined for each sample by plating a 1:10 dilution series in PBS on LB agar plates.
GyrAB quantification. To quantify the relative amount of GyrAB, samples of SoA3329 and SoA3330 were collected after 45 min of drug treatment as described above. An equal number of cells corresponding to 1 mL culture at OD 600 = 1 were harvested by centrifugation. Pelleted cells were lysed at room temperature for 20 min using B-PER bacterial protein extraction reagent (90078, Thermo Scientific) supplemented with 100 μg/mL lysozyme, 5 units/mL DNaseI (all part of B-PER with Enzymes Bacterial Protein Extraction Kit, 90078, Thermo Scientifc) and 100 μM/ mL PMSF (52332, Calbiochem). Samples were stored at -80˚C until further use.
Band intensities were quantified from unmodified images using the record measurement tool of Photoshop CS6, normalized to the CRP loading control after background subtraction, and reported relative to SoA3330. For clarity, the "levels" tool of Photoshop CS6 was used to enhance the contrast of shown Western blot images. We use the model fitted to experimental data to explore the sensitivity of our results to changes in the turnover rate of the drug-target complex. We changed k r and k f (0.01x, 0.1x, 1x, 10x, 100x, and 1000x original value) while keeping the ratio between k f and k r , the affinity K D , constant. a, shows the net growth rate (log 10 (bacterial number at 18 h)-log 10 target r(x). We use the model fitted to experimental data to explore the sensitivity of our results to changes in the replication rate with increasingly bound target r(x). We change the value of bound target at which we obtain a half-maximal replication rate, x 1/2 . a, Functions connecting bacterial replication rates r(x) to percentage of bound target molecules with different half-maximal replication rates. b, Net growth rate (log 10 (bacterial number at 18 h)-log 10 bound target δ(x). We use the model fitted to experimental data to explore the sensitivity of our results to changes in the death rate with increasingly bound target δ(x). We change the value of bound target at which we obtain a half-maximal death rate, x 1/2 . a, Functions connecting bacterial death rates δ(x) to percentage of bound target molecules with different half-maximal death rates. b, Net growth rate (log 10 (bacterial number at 18 h)-log 10 5) to explore the sensitivity of our results to changes in k a (0.01x, 0.1x, 1x, 10x, and 100x original value). a, Net growth rate (log 10 (bacterial number at 18 h)-log 10 (bacterial number at 0 h))/18 h) as function of drug concentration for different values of the binding rate k a (see legend). The dotted horizontal line indicates zero net growth. The intersections of the simulated dose-response curves with this line indicate the respective MICs. b, Sensitivity of the MIC to k a obtained from simulations in (a). The color code indicates the MIC corresponding to the simulation with the same color in (a). (TIF) S10 Fig. Sensitivity analysis of ampicillin fit: Changes in the drug-target turnover rate. We use the model fitted to experimental data (Fig 5) to explore the sensitivity of our results to changes in the turnover rate of the drug-target complex. We changed values for k a and k f (0.01x, 0.1x, 1x, 10x, and 100x original value) while keeping the ration of k a /k f constant. a, Net growth rate (log 10 (bacterial number at 18 h)-log 10  . We use the model fitted to experimental data ( Fig 5) to explore the sensitivity of our results to changes in the death rate with increasingly bound target δ(x). We change the value of bound target at which we obtain a half-maximal death rate, x 1/2 (see legend). a, Functions connecting bacterial death rates δ(x) to percentage of bound target molecules with different half-maximal replication rates x 1/2 (see legend). b, Net growth rate (log 10 (bacterial number at 18 h)-log 10 (bacterial number at 0 h))/18 h) as function of drug concentration for different values of δ(x) (see legend). The dotted horizontal line indicates zero net growth. The intersections of the simulated dose-response curves with this line indicate the respective MICs. c, Sensitivity of the MIC to δ(x) obtained from simulations in (b). The color code indicates the MIC corresponding to the simulation with the same color in (a&b). (TIF) S12 Fig. Predicted mutation selection windows for E. coli exposed to ampicillin. The drug concentration of ampicillin is shown on the x-axes, and the average bacterial net growth rate over 18 h is given on the y-axes. The blue line represents the wild-type strain based on the fits shown in Fig 5, and the red line represents a strain with a theoretical resistance mutation that decreases the binding rate (k f ) 100-fold and imparts a 15% fitness cost. The dotted horizontal line represents no net growth. The first vertical dotted line indicates where the resistant strain becomes fitter than the wild-type (the start of the competitive resistance selection window), the solid vertical line indicates the MIC of the wild-type (the start of the classical resistance selection window), and the dashed vertical line indicates the MIC of the resistant strain, above which selection for resistance should be minimal because both growth of the wild-type and the resistant strain is inhibited. (TIF) S13 Fig. Schematic of pharmacokinetic model. We simulate plasma and tissue concentrations of ampicillin with a two-compartment pharmacokinetic model. This model described intravenous drug input into the "plasma" compartment, which has an apparent volume of V 1 . From there, it can enter the peripheral "tissue" compartment, characterized by the apparent volume V 2 , with a rate k 12 . Conversely, the drug can also re-enter the plasma compartment with a rate k 21 . From the plasma compartment, the drug is eliminated with a rate k 10 . (TIF) S14 Fig. Coupling a pharmacokinetic model to COMBAT. Two modes of drug administrations, 2 g of drug per day given as single, 5 min i.v. infusions (blue line) and 2 g of drug per day given as continuous i.v. infusion (red line), are simulated. a, Simulated drug concentrations in the tissue (i.e., infected) compartment of a two compartment pharmacokinetic model over two days. The dotted black line indicates the MIC of the pathogen (2.6 mg/L). b, Pathogen load in the tissue compartment in response to fluctuating drug concentrations predicted by COMBAT over the same timeframe. (TIF) S1 Table. Kinetic parameters for used antibiotics. To our knowledge, the association rate for ciprofloxacin has not been determined directly. Because the values of the ratio of dissociation rate k r and association rate k f , K D , diverge by more than an order of magnitude in the literature, we chose to fit the association rate k f as a free parameter in our model while constraining K D to remain within the published range (the resulting value of k f is given together with other fitted parameters in S2 Table).  Table. Results of gyrase level determination and estimated GyrA 2 B 2 tetramer levels.

Supporting information
(m) indicates mild overexpression, (s) indicates strong overexpression. Columns headed GyrA and GyrB show the experimentally determined overexpression as fold expression compared to the wild type. The columns headed GyrA 2 B 2 show the estimated tetramer levels resulting from each measurement. For the GyrA 2 B 2 tetramer estimation, we sampled 10 4 sets association and dissociation rates from a uniform distribution within their reported limits (Latin hypercube approach). We report the standard deviation for each estimate. We give summary estimates in the last row of the table. (DOCX) S4 Table. Parameters resulting from model fit to experimental time-kill curves of ampicillin in wild-type E. coli. Fig 5B shows the resulting death rate δ(x) as an exponential function of the number of bound targets δ(x) = a 3 e b3x + c 3 . (DOCX) S5 Table. Parameters used in pharmacokinetic model. The parameters were obtained from [60]. (DOCX) S1 Text. Predicting pathogen load in patient by coupling COMBAT with a pharmacokinetic model. (DOCX)