How molecular mechanisms of resistance affect resistance evolution

Antibiotic treatment may be compromised by the sporadic appearance and selection of drug resistant mutants during therapy. Identifying optimal dosing strategies for treating bacterial infections that minimize the risk of resistance is difficult, and improving dosing guidelines usually requires long and costly in vitro and in vivo experimentation. Previously, we proposed a mathematical model that links bacterial population biology with chemical reaction kinetics and demonstrated that this model had high predictive and explanatory power for antibiotic pharmacodynamics. Here, we extend our model to incorporate several distinct molecular mechanisms of resistance, e.g. altered drug-target affinity or target molecule content, to explore the how these mechanisms affect the risk of acquired resistance during treatment. Our model is based on a system of partial and ordinary differential equations. The central assumption encoded in the model is that bacterial growth declines (and/or bacterial kill rate increases) as the fraction of bound antibiotic target molecules within the bacterium increases. We use experimental data from E. coli exposed to ciprofloxacin and ampicillin to estimate the dependence of bacterial growth and death on the fraction of bound target molecules and to calibrate our model. We then predict how several alternative molecular mechanisms of resistance should affect bacterial growth and compare these predictions with data from experimentally-manipulated bacteria with varying target molecule content. Finally, we use our model to estimate antibiotic susceptibility for different antibiotics that share similar mechanisms of action but differ in their drug-target affinity. Our extended model offers a novel approach for predicting how resistance will evolve during therapy and provides a framework for generating new testable hypotheses about how the risk of resistance depends on the molecular mechanism of resistance.


Introduction
The emergence and spread of antibiotic resistance is a major threat for global public health. Despite the magnitude of the problem, we have limited knowledge how to employ antibiotics in a way that maximizes treatment efficacy while minimizing resistance [1]. Whether optimal treatment outcomes are best achieved under higher or lower antibiotic exposure is a matter of intense debate [2][3][4]. As early as 1913, Paul Ehrlich recommended strategies that "hit hard" in order to maximize efficacy of antibiotic treatment [5]. The underlying intuition supporting such approaches is that high doses of antibiotics provide a obstacle to the emergence of resistance by suppressing the appearance of mutants that possess intermediate resistance; this "hit hard" recommendation has been supported by many experimental and theoretical studies. Conversely, high doses of antibiotics create a greater selective pressure in favor of extant resistant mutants, and in some situations "hitting hard" may result in a greater risk of antibiotic resistance [6]. The competing consequences of high dose antibiotic treatment (i.e. suppression of the rise of intermediately resistant mutants and selection of existing resistant mutants) have been consolidated in the concept of a resistanceselection window [7,8]. This window is the worrisome range of antibiotic concentrations in which growth of the susceptible (wild-type) strain is sufficiently impaired that resistant strains with a single resistance mutation can be expected to have a higher replication rate, i.e. higher fitness. At the same time, the antibiotic concentration is not high enough to eradicate all possible strains with a single resistance mutation. Antibiotic concentrations above the mutation-selection window are expected to safeguard against de novo resistance emergence during treatment. Antibiotic concentrations below the mutation selection window do not kill the susceptible strain, but also do not sufficiently favor the resistant strain to allow for out competition and emergence of resistance. While these concepts are useful, we note that it is rarely known the precise concentrations that select for resistance.
One challenge for predicting optimal antibiotic dosing to suppress resistance is that there are often several different mutations that can lead to resistance to a single drug [9], and each of these mutations may affect bacterial susceptibility and bacterial fitness in different ways. Tools that facilitate predictions of the susceptibility of the bacterial population from the possible molecular mechanisms of resistance would thus be useful. In this manuscript, we use a mechanistic mathematical model that describes the effect of intracellular reaction kinetics on bacterial populations [9] to investigate how different mechanisms of antibiotic resistance affect antibiotic susceptibility. We also consider how modifications to existing antibiotics, for example, small chemical modifications that increase their target affinity (as it has been done with each generation of quinolones [10]), may be leveraged to improve treatment outcomes.

Methods
We introduce a mathematical model based on a system of partial and ordinary differential equations. Similar to [1], the model describes the binding and unbinding of antibiotics to their target and how this affects bacterial growth and death. The model includes two mass balance equations, one for antibiotics and one for bacteria.

Binding kinetics
The total number of available target molecules θ is assumed to be a constant, the association rate is kf, the dissociation rate is kr and the number of bound targets x ranges between 0 and θ. In a first model, we ignore differences between extracellular and intracellular antibiotic concentrations, and only follow the total antibiotic concentration A and assume that the time needed for drug molecules to enter bacterial cells is negligible. This first model is therefore most appropriate for describing antibiotic action where the diffusion barrier to the target is weak (e.g. because the target is in the cell envelope as is the case for beta-lactams).
In [1], the association and dissociation terms are described by the following terms , and kf is the association rate, Vtot is the volume containing all bacteria and drug molecules, nA is the Avogadro number, , is the dissociation rate, . 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 differential equations ( + 2). To simplify and generalize this approach, we now assume that the variable of bound targets is a real number ∈ ℛ. 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.

Replication rate
We assume that living bacteria replicate with a rate r(x) depending on the number of bound target molecules. The function r(x) is assumed to be a monotonically decreasing function of x, i.e. the more target is bound, the less bacteria replicate. r(0) is the maximum growth 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 when the total bacterial population approaches the carrying capacity K. Thus, we can write the replication term as where Flim is the growth limiting term due to the carrying capacity K, and 0 ≤ B.C ≤ 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 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 − unbound molecules. Under the continuity assumption, we must generalize the concept of hypergeometric distribution. Since the hypergeometric distribution is a function of combinations, and 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 now be expressed as a probability density function of continuous variables.
In other words, the quantity of new bacteria is given by the term ( ) ( , ) B.C ( ), but twice this quantity must be redistributed along x. For illustration, if a mother cell with 4 bound targets dives, 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 which 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
Here, we define the death rate function ( ), which depends on the number of bound target

Bacteriostatic and bactericidal effects
We consider several functional forms of the relationship between the percentage of bound targets and growth and death rates ( Figure 1: step function, sigmoidal, linear). With sufficient experimental data, the growth rate r(x) and/or the death rate ( ) can be obtained by fitting our model to timekill curves of bacterial populations after antibiotic exposure. The sigmoidal shape of r(x) and d(x) can be written as: where xrth is the growth rate threshold, xdth is the death rate threshold, and both of them represent the point where the sigmoidal function reaches ½ of its maximum. gr and gd are the shape parameters respectively of the sigmoidal growth and death rate functions.

Full equation describing bacterial population
Finally, the full equation of bacteria is where B(x,t) is the number of bacteria.

Equation describing antibiotic concentration
The equation describing antibiotic concentration is given simply by the mass conservation, i.e. subtracting all the molecules associated in the complexes drug-target and adding all the molecules dissociated. The antibiotic equation is Where A(t) represents the number of drug molecules.

Initial and boundary conditions
To close the system of our equations, we must specify the initial and boundary conditions. At t=0, we assume that all bacteria have zero bound targets, i.e. = 0, and the initial drug concentration is C(0) = C0.
On the boundaries ( = 0, = ) of the partial differential equation, we specify that the outgoing velocities are zero. Thus, for = 0, we have , (0, ) = 0, and in = we have that # ( , ) = 0. Specifying this structure of our velocities permits us to solve the partial differential equation (using a finite differences scheme at the first order in space and time approximations) simplifying our needs on the boundary conditions; in other words we have the same conditions in = 0 and = present in [1]. To solve this system, we use a finite differences numerical scheme at the first order in space and time. Additional information and how this numerical scheme is developed can be found in [11,12].

Description of beta-lactam action
Beta-lactams acetylate penicillin-binding proteins (PBPs, the target molecules), and thereby inhibit cell wall synthesis. The acetylation of PBPs consumes beta-lactams, and therefore the drug-target reaction is not reversible. To include this effect in our model, we modify the term of drug-target dissociation. In the equation describing antibiotic concentration, the unbinding rate kr=0, while, in the bacteria equation, we substitute the dissociation rate kr with the deacetylation rate ka, see [13].

Model validation
First, we use time-kill curves with the beta-lactam ampicillin, for which most parameters have been experimentally determined (see table 1) to explore how well our modeling approach describes experimental data. It is difficult to determine r(x) and d(x) (the dependence of bacterial replication and death rates on the number of bound target molecules, respectively) experimentally. In single cell analyses, ampicillin was found to have no effect on bacterial replication rates at levels up to 80% of its MIC [14] and ampicillin is usually thought to be a bactericidal antibiotic [15]. We therefore assume that bacterial replication rates are not affected by ampicillin binding and that ampicillin action is entirely due to an increase in bacterial death rates with successively bound target molecules ( Figure 1C & D). We fitted our model parametrized with previously determined experimental data (Table 1) to time-kill curves of E. coli exposed to ampicillin (data from [16]) to obtain d(x). We used the following drug concentrations: 6 mg/l, 48 mg/l, 256 mg/l to fit the model and 0mg/l to estimate r0, the maximal growth rate in absence of antibiotics. We used simulated annealing to identify parameters of our model, finding the following optimal set of parameters (  The estimated mean occupancy considering both living and dead bacteria is 89%. This value is very close to experimental estimates in Staphylococcus aureus (86-99%) [17]. Thus, our model is able to describe bacterial time-kill curves from first principles with measurable parameters of drugtarget binding.
We next use our model to predict antibiotic efficacy (as measured by the effective growth rate [1/h] divided log10(e) to match the experimental measure of log10 change per hour) for other concentrations. We find the same characteristic sigmoidal curve for the dependence of effective growth rate on drug concentration that has been described earlier [16].

Changing drug-target binding kinetics and effects on antibiotic susceptibility
We next investigate how changes in physicochemical parameters will affect bacterial susceptibility to antibiotics. To this end, we consider the quinolone ciprofloxacin, which is one of the most commonly used antibiotics. The effect of varying target molecule content has been evaluated experimentally [18], which allows us to investigate resistance mechanisms that affect the target copy number, such as gene duplications. We can also study the effects of changes in drug-target affinity since the antibiotic class of quinolones contains drugs with a very wide range of drug target affinities (KD ~ 10 -4 -10 -7 ) with similar molecular size and similar modes of action. However, as (to our knowledge) for the majority of antibiotics, the exact rates of drug-target binding and unbinding as well as target occupancy at MIC are not known. We therefore first perform a uni-dimensional sensitivity analysis to explore how these unknown parameters will affect the predicted time-kill curves (supplementary material). Surprisingly, we find that the actual rates of drug-target binding, kf and kr (and not only the relative values, i.e. the affinity KD = kr/ kf) change antibiotic susceptibility profoundly. The faster the drug binds its target, the fewer bacteria can replicate before the antibiotic is fully active, and therefore the less antibiotic is needed to suppress bacteria ( Figure S2). While this finding is exciting and may guide future optimization of antibiotics, antibiotic binding rates are often not known and in these cases it may be difficult to fit our model quantitatively to time-kill curves.
However, our qualitative predictions regarding net growth at varying drug concentrations depending on target molecule content correspond well with experimental data where the target (gyrase) content of bacterial cells was altered by overexpression. We obtain a decreasing MIC for a higher number of targets when assuming that each bacterium can tolerate only a specified number of bound targets.

Estimating the MIC without data on binding rates
To explore qualitatively how changing drug-target binding kinetics affects the MIC, i.e. the minimal antibiotic concentration at which bacterial growth ceases, we simplify our model. As described above, the distribution of bacteria along the x-axis is due to the velocity v = vf + vr and to the replication rate. Here, we neglect the reproduction of free target molecules due to bacterial replication. Under this approximation, we use a lower drug concentration to inhibit the growth than we really need, because we neglect the reproduction of target molecules due to replication.
For this reason, we have an estimate of the lower boundary of the MIC.
The velocity v is an indicator of the action of the antibiotics. When this action is going to stop (equilibrium between association and dissociation), it means that the velocity is approached to zero.
We also know that the number of drug molecules A(t) has really small variations in time (almost always less than 0.1 %), so we can consider here A(t) = A(0), a constant value.
Under these assumptions, we can estimate the value x0 where the velocity becomes zero, by substituting v=0 in the equation above, obtaining At the MIC the death rate exceeds the growth rate. If we want to inhibit the bacterial growth, we need sufficient drug to exceed xMIC, the amount of bound target required to reduce bacterial growth to zero under idealized conditions (neglecting replication). This means I ≥ MNO . For this inequality to hold true, CD has to be: where # is the molar mass of the drug and Q = , / # is the affinity.
With the sign equal, we have UVBWX = Q . This is a lower boundary for the MIC.

Qualitatively, this means that the MIC is linear with KD when keeping all other parameters
constant. This is consistent with our previous finding when deriving the MIC from drug-target binding alone in a very simplified model [13].
We tested this hypothesis with pooled MICs from clinical enterobacteriaceae isolates collected before the wide-spread use of quinolones in the clinic. Here, we assume that the target molecule, gyrase, is sufficiently conserved in enterobacteriaceae that the relative drug-target affinities of quinolones are similar across these species. We correlate the KD of these quinolones (measured in E. coli) to the MIC and find that the MIC is indeed proportional to KD in 11 enterobacteriaceae species with R 2 between 0.88 and 0.98 ( Figure 4). retrieved from [19,20]. Quinolone-gyrase affinities were obtained from [19][20][21][22]. Linear regression was performed with the software package R. window. Here, we return to our fit of the experimental data on ampicillin to determine whether our model reflects the dynamics of a susceptible bacterial population. We compare predicted net growth rates depending on ampicillin concentration between a wild-type strain and an arbitrary resistant strain. We assume here that this resistant strain has a 50-fold slower drug-target binding rate, and the maximum growth rate is the 85% of the maximum growth rate of the wild type. We then predict the antibiotic concentrations at which resistance is selected. In this case, the "competitive resistance selection window", i.e. the concentration range below the MIC for both strains where the resistant is fitter than the wild type, would start at 0.66 mg/l and end at 3.32 mg/l. The "resistance selection window", i.e. the concentration range above MIC for the wildtype but below the MIC for the resistant strain, would start at 3.32 mg/l and end at 154.5 mg/l. line is the wild type behavior, the red line is the resistant type. We predict a "competitive resistance selection window" between 0.66 mg/l and 3.32 mg/l. A "resistance selection window" between 3.32 mg/l and 154.5 mg/l. For a concentration, higher than 154.5 mg/l we observe both strains suppressed.

Discussion
The importance of the MIC of both wild-type and resistant strains in identifying optimal dosing strategies to minimize the risk of resistance is widely appreciated [9]. It is advisable to avoid drug concentrations above the MIC for the wild type but below MIC for the resistant strain because these concentrations allow only the resistant strain to grow. In addition, selection for resistance may also occur below the MIC of the wild type in if antibacterial concentrations are in the so-called "competitive resistance selection window", where resistant strain may outcompete the wild type [8]. However, we do not know a priori how a given resistance mutation affects bacterial growth at a given antibiotic concentration. Here, we present a mathematical model that allows in silico predictions of ranges of antibiotic concentrations that should be avoided, and how this range depends on the molecular mechanism of resistance.
To validate our model, we fitted it to three published time-kill curves of E. coli exposed to ampicillin. All physicochemical characteristics needed to parametrize our model were determined experimentally for ampicillin. While it was shown experimentally that at MIC, 84-99% of the target molecules are bound, we do not know how bacterial death rates depend on target occupancy. Therefore, we parameterized our model with all other known characteristics and fitted the death rate as a function of occupied target, δ(x), to the data. Our model fits the experimental data very well (R 2 =0.93-0.98). Our model suggests that the death rate may have a sigmoidal-shaped dependence on target occupancy. This death rate function δ(x) predicts that at MIC, 89% of the target is occupied, which is in excellent agreement with the experimental data. We then went on to predict bacterial net growth rates across a wide range of concentrations and predict an MIC of 3.32 [mg/l]. The dataset we worked with was used previously to obtain more accurate MIC estimates than those obtained with classical methods. The previous estimate for the MIC based on 7 timekill curves was 3.4 [mg/l], which is very close to our estimate even though we used only three time-kill curves. We therefore conclude that our model is able to describe bacterial time-kill curves quantitatively based on chemical reaction kinetics.
Our next goal was to investigate how different physicochemical characteristics shape bacterial susceptibility depending on antibiotic concentrations. To this end, we performed a uni-dimensional sensitivity analysis. Surprisingly, we find that the binding (kf) and unbinding rates (kr), and not only the binding affinity (KD = kr / kf) has a profound effect on bacterial susceptibility. When keeping KD constant, drugs that bind more quickly are more effective because they prevent bacterial replication before the target molecules are reproduced. From a bacterial perspective, resistance mutations that slow down drug-target binding are more beneficial than those that speed drug-target dissociation. We then investigated whether our model can reproduce known effects of changes in target molecule content. We used strains overexpressing the target of the antibiotic ciprofloxacin, gyrase, and exposed them to a large range of ciprofloxacin concentrations. As previously shown, bacterial susceptibility increases with increasing target molecule content; our mathematical model is consistent with this finding [18]. This increases our confidence that our model is able to describe antibiotic efficacy when the intracellular reaction kinetics of drug-target binding changes. However, in our experiments we also observe an increased kill rate at very high concentrations of ciprofloxacin. This may be due to increasing DNA damage when ciprofloxacin binds gyrase, a mechanism not explicitly incorporated in our model. For many antibiotics only the drug-target affinity and not the actual binding rates have been published. While it is difficult to fit our model to time-kill curves in this case, we can simplify the model and make qualitative predictions about expected changes in MIC. Our model predicts that the MIC should be proportional to drug-target affinity. The antibiotic class of quinolones contains drugs with a wide range of drug target affinities (KD ~ 10 -4 -10 -7 ) with similar molecular size and similar modes of action, and we hypothesize that the MIC of different quinolones should be proportional to their gyrase affinity. We verified this hypothesis with data from a meta-analysis on the MICs of various quinolones in clinical enterobacteriaceae isolates collected before the widespread use of quinolones. This means that even in absence of measured binding rates, our model can be employed to understand qualitatively e.g. how much antibiotic to give to prevent resistance mutations that lead to a given change in drug-target affinity. It can also serve as a tool in drug development to predict how much drug-target binding has to be improved to achieve efficacy at a certain antibiotic concentration.
Finally, we demonstrate how our model could be used to predict resistance selection windows, thereby allowing designing drug dosing strategies that avoid exposures to antibiotic concentrations that select for resistance. Our model does not consider host immune responses to infection, which may further modify our expectations regarding treatment success [23][24][25]. Nevertheless, given the urgent need to preserve the efficacyof existing antibiotics and the need to develop new agents [26], we see a promising role for mechanistic models that can suggest the most promising dosing strategies based on the physicochemical and biochemical characteristics of drug-target interactions.
Taken together, our model offers a novel in silico platform that may aid in finding antibiotic treatment strategies that minimize the emergence of resistance.