Stochastic Gene Expression Influences the Selection of Antibiotic Resistance Mutations

Abstract Bacteria can resist antibiotics by expressing enzymes that remove or deactivate drug molecules. Here, we study the effects of gene expression stochasticity on efflux and enzymatic resistance. We construct an agent-based model that stochastically simulates multiple biochemical processes in the cell and we observe the growth and survival dynamics of the cell population. Resistance-enhancing mutations are introduced by varying parameters that control the enzyme expression or efficacy. We find that stochastic gene expression can cause complex dynamics in terms of survival and extinction for these mutants. Regulatory mutations, which augment the frequency and duration of resistance gene transcription, can provide limited resistance by increasing mean expression. Structural mutations, which modify the enzyme or efflux efficacy, provide most resistance by improving the binding affinity of the resistance protein to the antibiotic; increasing the enzyme’s catalytic rate alone may contribute to resistance if drug binding is not rate limiting. Overall, we identify conditions where regulatory mutations are selected over structural mutations, and vice versa. Our findings show that stochastic gene expression is a key factor underlying efflux and enzymatic resistances and should be taken into consideration in future antibiotic research.


Introduction
Efflux pumps and drug-inactivating enzymes allow bacterial cells to evade the damage caused by antibiotic drugs. Bacterial pathogens with efflux and enzymatic resistances are ubiquitous and often are serious public health concerns (Li et al. 2015;World Health Organization 2017). For instance, several Gram-negative bacteria species with carbapenemase and ESBL activities are listed as critical priority pathogens for research and development of new antibiotics by the WHO (World Health Organization 2017).
Studies to quantify antibiotic resistance are routinely performed at the bacterial population level, that is, through protocols such as minimum inhibitory concentration (MIC) or time-kill curves. However, the action of the drug occurs within the cell on the molecular level. Bridging these different scales remains a challenge.
The action of antibiotic drugs usually involves small numbers of molecules and binding sites, such that intrinsic stochasticity could have a significant effect on the dynamics. This noise impacts how we interpret bacterial assays that extract population-averaged behaviors, which normally ignore cellto-cell heterogeneity. A recent study has suggested that many experimental phenomena, such as postantibiotic effects, can be explained by these noisy within-cell dynamics (Abel Zur Wiesch et al. 2015).
Another source of randomness is in the mechanism of gene expression (Paulsson 2005). Transcription typically occurs in bursts, a phenomenon often described by the socalled two-state or telegraph model of gene expression (Paulsson 2005;Lionnet and Singer 2012;Kumar et al. 2015). In this model, gene transcription switches stochastically between an active state with a constant rate of mRNA production, and an inactive state without transcription. Although certain genes may have complex transcription regulations, the two-state model is widely used when modeling both prokaryotes and eukaryotes due to its simplicity (Paulsson 2005). The duration of active and inactive states can range from minutes to multiple cell generations (So et al. 2011;Lionnet and Singer 2012;Hammar et al. 2014). The intrinsic stochasticity of gene expression is often the dominant source of randomness for gene-product numbers, and the separation of noise intensities has been used to perform mathematical analysis of gene-expression systems (Lin and Galla 2016). While the expression level of resistance enzymes or efflux pumps generally correlates with the level of observed phenotypic resistance (Zwart et al. 2018), the effects of stochasticity on resistance remain unexplored.
Enzymatic resistances can be enhanced by regulatory or structural gain-of-function mutations, affecting the gene expression or the enzyme efficacy, respectively (Yang et al. 2003;H€ andel et al. 2014;Berrazeg et al. 2015;Blair et al. 2015). Both regulatory and structural mutations have been described for a number of antibiotic resistances (Toprak et al. 2012;H€ andel et al. 2014;Berrazeg et al. 2015). However, the selective conditions that favor one type of mutation over the other have not been fully explored.
Here, we provide a computational model that investigates the effects of stochastic gene expression on resistance and resistance evolution. Our model describes the dynamics of within-cell processes of the drug-target-efflux system and accounts for stochasticity in transcription and translation, as well as drug diffusion, binding, and removal. Using this model, we consider regulatory and structural resistance mutations to study their behavior on the molecular, cellular, and population levels. Our goal is to explore how stochastic gene expression influences the survival and extinction of different types of resistance mutations under antibiotic treatments.

Overview
We model a population of elongating and dividing cells in a boundless environment with an antimicrobial drug. The drug molecules can diffuse across cell membranes, and they influence the cell's death rate by binding to drug-specific targets. A drug efflux system, which is subject to stochastic gene expression, removes drug molecules from the cell. We focus on mutations that affect this system, either through regulatory effects or changes to the enzymatic efficacy.
Each cell is modeled as an independent agent, such that there are no between-cell interactions. To uphold independence, we assume a constant concentration of drug outside the cell. Spatial effects are ignored. The within-cell model takes the following processes into account: cell elongation and division; intracellular production of drug-target protein; entry of drug molecules into the cell; the interaction of target proteins with drugs and subsequent cell death due to antibiotics; and the transcription, translation, and enzymatic activity of the resistance-mediating proteins. From here on, we refer to the resistance-mediating protein as efflux protein, as it can remove drug molecules from the cell either as an efflux pump or drug-inactivation enzymes.
Each cell is described by seven discrete variables (shown in table 1), which are updated stochastically using the adaptive tau-leaping Gillespie algorithm as implemented in the adaptive tau package in R (Cao et al. 2007). This is an approximation of the full stochastic simulation algorithm (Gillespie 1977). The interactions between the within-cell variables can be seen in the model schematic in figure 1 and are described in detail below. In addition, we model the elongation of the cell deterministically. All cell parameter values are reported in supplementary table S2, Supplementary Material online.

Cell Physiology and Division
The cells we model here are Escherichia coli. We assume that the cells are cylindrical with cross-sectional diameter d. A newly divided cell has length ' 0 . The cell elongates exponentially until it reaches twice its length at birth (Campos et al. 2014), which is achieved in time interval t G . At a time t after birth, the cell length ('), surface area (A), and volume (V) satisfy: For each cell, a small amount of zero-expectation Gaussian noise (SD ¼ 5%) is added to the mean generation time (t G ! N ½t G ; ð0:05t G Þ 2 ) to desynchronize cell divisions and to increase cell-to-cell heterogeneity. Division happens instantaneously once the cell reaches twice its birth length, and produces two equal-length daughter cells. All variables in table 1 are divided binomially between the two daughter cells upon cell division, with the exception of the DNA activity D E , which is directly inherited by both daughters.

Target Production
The target proteins usually fulfill essential functions in the cell, such as ribosomal subunits, RNA polymerase, DNA gyrase, and cell-wall components (Andersson and Hughes 2010). The numbers per cell of several targets have been well-characterized and used for computational studies (Abel Zur Wiesch et al. 2015). The intracellular concentrations of these essential proteins may impact the cell's growth, for instance as in the case with ribosomes (Greulich et al. 2015). Therefore, we assume that these target proteins are subject to various intracellular regulations that maintain a stable and optimal concentration within the cell. This assumption is consistent with previous findings showing essential genes have low protein-expression noise (Silander et al. 2012). In the model, this is achieved through constant production of targets over the duration of the cell's life. The production rate C is chosen such that, on an average, the number of targets doubles in a generation. That is, , such that: As the half-life of bacterial proteins is usually on the scale of multiple hours (Larrabee et al. 1980), we neglect protein

Efflux Production
The transcription and translation of efflux proteins are explicitly modeled: we expect the expression of these nonessential genes to be more stochastic than that of the drug targets. The gene activity switches on(off) with a fixed rate C On (C Off ), and therefore follows the telegraph process of gene expression (van Kampen 2007). Once the gene is activated (D E ¼ 1), mRNA is transcribed with a constant rate s. The mRNA molecules degrade with rate c. On an average, b proteins are translated during the lifetime of each mRNA molecule. As above, we assume that these proteins do not decay but are diluted only by cell division. Given that these processes are intrinsically stochastic, a brief period of gene activation may not always produce efflux mRNA and/or protein.

Drug Diffusion, Interaction, and Efflux
The antibiotic drug outside of the cells is kept at a constant concentration c out . The drug molecules diffuse through the cell envelope bidirectionally with a fixed diffusion rate r. The overall diffusion rate is proportional to the surface area of the cell, A(t). We stochastically model Fick's law of diffusion for the bidirectional movement, such that transport is driven by a concentration gradient. Once a drug molecule has entered the cell, it can either diffuse out, bind to its specific target, or be captured by an efflux protein. These binding events result in the formation of a drug-protein bound state. Binding occurs with the forward rate constant k f , reflecting the binding affinity of target and efflux proteins. The actual binding rate also depends on the concentration of the drug and the proteins that interact inside the cell. Drug binding to target and efflux are both considered reversible unless stated otherwise, and the bound states dissociate with backward rate k b . This event releases both the drug molecule and the protein into the cell. To compete efficiently with targets, efflux proteins should have at least the same or higher drug-binding rates than target proteins. Therefore, we set the basal efflux binding rate to be the same as that of targets, making them indistinguishable to the drug. The efflux protein has the capacity to remove the drug molecule once bound, through catalysis or transport, and this happens with rate k cat .
All of the above reactions are summarized in table 2.

Cell Death and MIC Fraction
Here, we consider bactericidal drugs that introduce a drugdependent death rate for each cell. Therefore, the presence of drug has no impact on the cell generation time. We combine two previously published models to implement a realistic death mechanism: the drug-induced death rate, logð10ÞdðqÞ, is an increasing function of the fraction of bound targets, q (Abel Zur Wiesch et al. 2015), but has a maximum value (Regoes et al. 2004). Under these conditions, the number of cells n follows nðtÞ ¼ 10 ½w max ÀdðqÞt nð0Þ, where w max is the population growth rate in the absence of antibiotic. The fraction of bound targets, q, also correlates with the concentration of bound targets. The death rate should satisfy the following criteria: (1) dð0Þ ¼ 0: In the absence of antibiotics the death rate is zero; (2) dðq MIC Þ ¼ w max : When the external drug concentration is 1ÂMIC-which is the minimum concentration at which no bacterial growth is observed, the fraction of bound targets is given as q MIC . The death rate at this MBE value should equal the maximum growth rate in the absence of antibiotic (w max ), such that the net growth rate is zero; (3) dð1Þ ¼ w max À w min : When all targets are bound, the death rate is maximal. The net growth rate is then w min , which is the minimum measured population growth rate.
In our model, the death rate takes the sigmoidal form as in Regoes et al. (2004), where the constants A and B are chosen such that the above criteria are satisfied, that is, The values of w max and w min have been determined for multiple antibiotic compounds (Regoes et al. 2004). Finally, the shape parameter j in equation (3) determines the steepness of the sigmoidal response. For j ! 1, we recover the step function where death only occurs if the fraction of bound targets is greater than a threshold (in this case, the threshold would be q MIC) .
Through equation (3), we can relate the within-cell parameter q MIC to the externally measured (population level) MIC of the drug. Empirically, MIC is determined as the lowest drug concentration to prevent a small bacterial inoculum from proliferating overnight into visible density, sometimes quantified specifically as > 90% reduction in growth rate compared with a drug-free control (Arthington-Skaggs et al. 2002;Tunney et al. 2004;Toprak et al. 2012). Alternatively, one could define the MIC as the concentration of drug at which the net growth rate is zero after a given time period; this is referred to as zMIC by Regoes et al. (2004). In the Supplementary Material online, we measure these quantities in simulations and compare them with different lineage survival probability thresholds following overnight experiments, such as IC50 (50% lineage survival), IC90 (10% lineage survival), and IC99 (1% lineage survival) (supplementary fig. S1, Supplementary Material online). The reason for using lineage survival probability is its computational efficiency compared with computing population growth rates. We find that IC90 is a suitable measure of MIC which correlates well with the zMIC drug concentration.
We determined values of q MIC by simulating lineages emerging from individual wild-type (WT) cells that grow for 20 h in a constant environment with c out ¼ 1ÂMIC. This time interval is comparable with standard overnight MIC experiments (Regoes et al. 2004). We then screened for q MIC, defined as the highest value of bound target fraction where >90% of the progenitor cells and their lineages become extinct. This screening process is highlighted in supplementary figure S2, Supplementary Material online, for the drug ciprofloxacin, from which we find an MIC fraction of 8.1%. Although this number may seem low, it could potentially be explained by the drug mechanism: ciprofloxacin binds to DNA-bound gyrases and fragments the bacterial chromosome via DNA double-strand breaks (Kampranis and Maxwell 1998;Tamayo et al. 2009). Given that each E. coli cell has $300 DNA-bound gyrases (Chong et al. 2014), this MIC fraction suggests that a cell is likely to die when its chromosome is fragmented by over 20 simultaneous doublestrand breaks, which seems plausible. The drug-specific parameters are listed in supplementary table S3

Simulation Algorithm
To simulate a population of cells over multiple generations, we use the following algorithm: (1) Assign each initial cell a unique set of parameter values; (2) Identify the cell with the earliest birth time; (3) Run the adaptive tau-leaping simulation algorithm for this cell until the cell dies, or it reaches its predefined generation time; (4) If the cell survived, perform cell division to create two new daughters;   (2), until the experiment ends, no cells remain, or the number of cells is large enough that we can assume survival of the population for the duration of the experiment (200 unless otherwise stated).
We also describe the within-cell behavior by a system of ODEs, as shown in the Supplementary Material online. We refer to these equations as the mean-field solution, as it does not take into account any stochasticity, and instead reflects the average dynamics of an infinite number of cells.

Model Behavior
To test the validity of our model, we checked the distribution of target and efflux molecules across a large ensemble of WT cells growing in the absence of antibiotics (supplementary fig.  S4, Supplementary Material online). We found that the frequency of cells with active efflux DNA agrees with the prediction of the telegraph process ($2:5%), while the average number of efflux mRNA ($0:04) and efflux proteins ($4) per cell are within experimentally reported ranges (Cai et al. 2006;Taniguchi et al. 2010;Moran et al. 2013). The efflux protein copy number also follows a negative-binomial distribution, as predicted by the theory of their bursty dynamics (Paulsson and Ehrenberg 2000). Finally, the number of target proteins per cell is normally distributed, which is expected given its constant production rate.
In figure 2, we show an example trajectory of the intracellular model variables in the absence and presence of an antibiotic (ciprofloxacin). Starting from a single WT cell, we track only one daughter cell after each division. After six generations, we instate an external drug concentration, such that the drug can freely diffuse into the cell ( fig. 2A). The mean-field equations closely approximate the number of constantly produced intracellular drug targets ( fig. 2B). The number and fraction of drug-bound targets is also well approximated in the mean-field limit ( fig. 2C and D). However, the bursty dynamics of the efflux system are not well captured by the deterministic approximation ( fig. 2E-H). Here, a piecewise deterministic Markov process, which accounts for the gene-expression noise, would be a more appropriate approximation for this system (Lin and Galla 2016).
Although cells could accumulate multiple efflux proteins during the period of gene activity ( fig. 2G), we did not observe two or more drug-bound efflux proteins at the same time ( fig. 2H). The lack of multiple drug-bound efflux proteins could be caused by two factors: the binding rate of drug is slow enough to disallow near-simultaneous formation of multiple drug-efflux bound states, while the catalysis rate is high enough that any efflux-bound drug molecule is quickly catalyzed. Put together, it shows that here the rate-limiting factor for removing drugs by efflux is not the catalysis step, but the drug-binding process (at least for this set of parameters).

Mutations
We consider seven classes of cells throughout this study: knockout, WT, three types of regulatory mutants, and two types of structural mutants. Knockout (KO) cells lack efflux protein production. WT cells have the basal level of efflux production. The effects of a mutation are characterized by a multiplicative parameter l > 1. The three types of regulatory mutants we consider are: REG-ON, which has an increased rate of gene activation (C On ! lC On ); REG-OFF, which has longer periods of active transcription by reducing the inactivation rate (C Off ! C Off =l); and REG-BURST, which increases the translation rate or protein burst size from each translational event (b ! lb). The two types of structural mutants include STRUCT-BIND with improved binding affinity (k f ! lk f ) to drug molecules and STRUCT-CAT with improved catalytic rates to break down drug molecules (k cat ! lk cat ). The mutant cells also carry a cost, 0 < ( 1, which is associated with their modified function. We assume this cost affects all metabolic processes in the cell, leading to a longer generation time, as well as slower translation.
We therefore implement C ! ð1 À ÞC; b ! ð1 À Þb, and ht G i ! ht G i=ð1 À Þ. The KO cell type carries a negative cost, as we assume it grows slightly faster for lacking efflux protein production completely.

Lineage Survival
We first simulate the survival probability of cell lineages when faced with a constant concentration of ciprofloxacin for 20 h. As well as considering different mutant classes, we also vary the mutant effect parameter l. KO and WT classes are included as controls.  For the cell lineages that became extinct in figure 3, we looked at their extinction times (supplementary fig. S6, Supplementary Material online). Like the WT and KO cell types, the structural mutants and REG-BURST have their peak mean extinction time near the WT MIC. If extinction occurs at sub-MIC concentrations (which is a rare event, as shown in fig. 3), it is most likely to happen at the beginning of the simulation when cell numbers are small. At 1ÂMIC, the extinction time diverges as death rate equals growth rate at this drug concentration (so that the expected lineage lifetime becomes infinite). As the drug concentration increases beyond the MIC, more drug diffuses into the cells and the death rate increases, leading to faster extinction of the cell lineages. Compared with WT and KO, both REG-ON and REG-OFF mutants show a shift of the peak extinction time toward higher concentrations with larger mutational effects. This corresponds to the shift in their respective IC90 values. Finally, we note that the STRUCT-BIND lineages of cells show a broad distribution (large variance) of extinction times at high drug doses.

MBE
Combining the data on survival probability and extinction times, our results suggest that REG-ON mutations provide intermediate but homogeneous resistance. This is characterized by the high survival probability at low drug concentrations, which declines quickly as the drug concentrations increase until all cells eventually die. Resistance resulting from a REG-OFF mutation, however, seems heterogeneous. This is characterized by prominent deaths even at low concentrations but with some cell lines surviving through intermediate drug doses. The heterogeneity is even more pronounced in the STRUCT-BIND mutants: the low extinction times show that some cells are dying early, but the continued survival through intermediate and high doses suggest that a fraction of the cells is highly resistant. For STRUCT-BIND, the extreme resistance heterogeneity should be the combined result of low but highly stochastic gene expression and elevated binding rates of the resistance enzyme, such that the few cells with copies of the improved efflux protein are highly protected.

Growth Rates
An alternative measure of mutant performance is through their population-level growth rates across different drug concentrations. As well as measuring the level of resistance conferred, these results can be directly compared with empirical observations. Following the protocol of Regoes et al. (2004), we simulated populations of cells for a short duration and counted the number of live cells at 10-min intervals, before extracting the net population growth rate ( fig. 4).
Both regulatory mutants show similar dose-response profiles, which gives no indication of the differences in their lineage survival probability. STRUCT-BIND mutants, however, generally maintain a higher net growth rate under high drug pressure when compared with the other cells types. This is in agreement with figure 3, where a few cell lineages are able to proliferate at these high treatment intensities. STRUCT-CAT, WT, and KO cells show the classical sigmoidal dose response, and the WT dose-response curve agrees well with the data reported by Regoes et al. (2004) (fig. 4E).

Heterogeneity and Resistance
The number of efflux proteins per cell in REG-ON and REG-OFF mutants are nearly identical in mean for a given mutation effect l ( fig. 5). However, these mutant classes have very different distributions of efflux protein number across the population. For REG-ON mutants, the protein number distribution is unimodal such that cells are more homogeneous; with frequent transcription bursts (e.g., at l ¼ 50 and 200) all cells produce some efflux proteins. REG-OFF mutants, in contrast, display a bimodal distribution of protein number: A large fraction of cells within the population do not produce any efflux protein, while some produce more than the corresponding REG-ON mutant. This is the cause of REG-OFF cell death even at low drug concentrations, but with some survival at intermediate doses.
Overall, regulatory mutations only provide limited resistance. For REG-ON and REG-OFF with 200Â mutational effects, transcription is active for $80% of the time. Therefore, resistance cannot be improved significantly by augmenting the transcription dynamics beyond this level. Further improvements in resistance then rely on increasing the translation rate (resulting in more proteins per mRNA) or 10 −1 10 +0 10 +1 10 +2 10 −1 10 +0 10 +1 10 +2 Stochastic Gene Expression . doi:10.1093/molbev/msz199 MBE by improving the efficacy of the efflux protein itself. This is consistent with previous findings showing that resistance provided by overexpression is limited (Yang et al. 2003;Toprak et al. 2012).
For STRUCT-BIND mutants, the efflux protein number distribution is very similar to that of the WT as there are no changes in expression (apart from the small cost of carrying the mutation). Therefore, the majority of STRUCT-BIND cells in a population contain no efflux protein. However, those cells that do contain efflux proteins have a high level of resistance, due to the high efficacy of these proteins. This is the source of heterogeneity in the STRUCT-BIND mutant population.
To investigate this further, we consider the survival for a STRUCT-BIND mutant, but with different levels of gene expression stochasticity ( fig. 6). Concretely, we increase the frequency of gene activation (C On ! aC On for a > 1), and correspondingly decrease the efflux protein translation rate cb ! cb=a. In this way, we maintain the mean number of efflux proteins per cell ( fig. 6A) but alter the expression noise.
By reducing the gene expression noise (larger a), we note an increased IC50 but decreased IC90 ( fig. 6B), as the survival profile now resembles that of REG-ON mutants (see fig. 3A). With reduced expression noise, the efflux protein number in cells becomes more uniform and cells with extreme protein abundances disappear, making resistance homogeneous across the cell populations. Overall, this shows that while noise control reduces the amount of drugs needed to inhibit the cells, the antibacterial effect also diminishes when drug concentrations fall below the MIC.
Another source of heterogeneity stems from the partitioning of proteins at cell division. Our model assumes efflux proteins are binomially distributed between two daughter cells, and so far we have only considered the symmetrical case with binomial parameter P ¼ 0.5. It has been previously shown, however, that biased partitioning for efflux pumps exists and could cause long-lasting phenotypic heterogeneity (Bergmiller et al. 2017). To check the effects of partitioning bias, we varied the probability parameter P of the binomial distribution to favor one daughter cell with more efflux proteins at division (supplementary fig. S8, Supplementary Material online). Among WT lineages, IC90 is unaffected since most cells have no efflux proteins, while a few cells have some with low efficacy. In REG-ON mutants, however, we observe a significant increase in the lineage survival probability (IC90), as one cell can now accumulate a significant number of efflux proteins. In the STRUCT-BIND mutants, one cell in the lineage will accumulate highly effective efflux proteins, making it into a superresistant cell and hence increasing the IC90. These results initially carry across to zMIC as determined from growth rate measurements-increasing the efflux distribution bias results in a higher zMIC. However, in the most extreme case (P ¼ 1) where only one daughter receives all the efflux protein, the STRUCT-BIND mutant shows a decreased zMIC measurement. Here, at drug concentrations above the WT-MIC, only one cell maintains the resistance from generation to generation, and hence the maximum growth rate is zero. Our results therefore suggest that there is an optimum strategy to invest more in the fitness of one daughter which maximizes the resistance of the population. This is not seen in the REG-ON mutant as almost all cells are producing efflux mRNA and proteins frequently. REG-ON daughter cells that inherited no efflux protein could replenish their resistance rapidly and survive the antibiotic, thereby contributing to population growth.
Although it is commonly accepted that increased mean expression of resistance enzymes can enhance resistance (Zwart et al. 2018), our results suggest that expression noise is also important for the population-level dynamics.

Survival after Antibiotic Pulse
The survival and extinction time properties of the different mutants are expected to affect the outcome of pulsed antibiotic treatment, which reflects the situation in a patient where drugs are not administered in a manner that maintains a constant concentration. Therefore, we simulated pulsed treatments where cells are exposed to a one-time dose of antibiotic for a limited time period. For simplicity, the drug dose is expressed as a step-function without more-detailed pharmacokinetics. We then measured the survival probability for each mutant class across different treatment durations and concentrations ( fig. 7). This survival probability can also be interpreted as the genetic diversity, as it measures the fraction of unique lineages which survive the pulsed treatment.
REG-ON mutants have the highest survival probability at low drug concentrations (c out 2ÂMIC) or with brief drug pulses (e.g., t ¼ 30 min), as shown in figure 7A. This is due to their long extinction times that allow regulatory mutants to outlast drug concentrations beyond their IC90. STRUCT-BIND cells survive at higher concentrations with long drug pulses (

Systematic Analysis
So far, we find that STRUCT-BIND mutations provide most resistance, followed by REG-ON. STRUCT-CAT, however, do not increase the resistance against ciprofloxacin. For betalactams, it has been shown in experiments that increasing k cat provides a high level of resistance (Knies et al. 2017;Palzkill 2018). To check whether our results hold true for other conditions, we performed a systematic grid sampling of the model parameter space. Specifically, we vary six model parameters: binding rate, catalysis rate, diffusion rate, number of targets, number of efflux proteins, and number of drug molecules per cell at MIC, constructing 729 unique parameter combinations. For each of these, we consider the WT cell, as well as the REG-ON, STRUCT-BIND, and STRUCT-CAT mutants with an effect of l ¼ 200. The performance of each mutant class is quantified as their ability to increase the IC90 of the cell population relative to the WT cell. Full details of this procedure can be found in the Supplementary Material online.
Through this analysis, we find that the largest effect mutations are STRUCT-BIND, but on an average REG-ON outperform the other mutants ( fig. 8). REG-ON and STRUCT-BIND mutants frequently provide increases of >10-fold in IC90, MBE while STRUCT-BIND mutants could even reach >100-fold in some circumstances. STRUCT-CAT mutants have little-to-no effect most of the time, but it is possible for them to have up to 10-fold advantage over the WT. Concretely, STRUCT-CAT mutants perform better when: 1) the binding rate is high; 2) the average number of efflux pumps per cell is high; 3) the number of intracellular drug molecules are high; and 4) the WT catalysis rate is low (supplementary fig. S10, Supplementary Material online). In all four scenarios, catalysis becomes the rate-limiting step in the efflux of drug from the cell. In summary, while STRUCT-CAT mutations could enhance resistance, their contributions are fairly limited compared with STRUCT-BIND and REG-ON.

Discussion
Efflux pumps and drug-inactivating enzymes cause some of the most important clinical resistances (Li et al. 2015;World Health Organization 2017). Overexpression of efflux pumps is also associated with virulence and enhanced mutation rates (Wang-Kan et al. 2017;El Meouche and Dunlop 2018). Although stochasticity is inherent to gene expression and underlies vital cell functions, its effects on antibiotic resistance have not been systematically addressed.
To study this, we constructed an efficient agent-based model that bridges the scales between within-cell molecular processes and population-level dynamics. We considered different classes of regulatory and structural mutations that enhance resistance, and evaluated their performances under different drug pressures. Previous models have utilized drug-target binding dynamics to help explain populationlevel phenomena such as persistence (Abel Zur Wiesch et al. 2015). Our model extends this approach by adding enzymatic reactions. Furthermore, unlike previous deterministic models, our stochastic approach explicitly captures the noise at each step of the relevant biological processes, allowing us to dissect the effects of each source of noise on bacterial behavior.
Regulatory mutations that increase either transcription frequency (REG-ON) or duration (REG-OFF) by the same factor have the same mean expression. Thus, a mean-field deterministic model would predict equal levels of resistance. However, increasing the duration of transcription does not reduce the noise of gene activation, leading to a subpopulation of cells with little or no resistance enzymes which remain drug-susceptible. This is not observed in the mutants with frequent gene activation, which instead show a homogeneous population of cells with equal resistance.
Structural mutations improve the biochemical properties of the resistance enzyme by increasing the binding rate to the drug (STRUCT-BIND) or the catalysis rate (STRUCT-CAT). While structural mutations may influence both binding and catalysis simultaneously, here, we disentangle them to elucidate their relative importance in resistance. We find that improving the binding rate has a significant effect on resistance, while improving the catalysis rate may have no effect (at least for our parametrization of E. coli exposed to ciprofloxacin).
To understand the robustness of these conclusions, we systematically investigated a large range of parameter space taking into account different target, efflux, and drug dynamics. Increasing resistance through faster catalysis (STRUCT-CAT) is effective when the number of efflux enzymes is greater than or equal to the number of targets, which may be the case when considering beta-lactamase resistance where increased catalytic rates have been detected (Knies et al. 2017;Palzkill 2018) in mutants with large increases in MIC (Hall 2002;Schenk et al. 2012). Furthermore, STRUCT-CAT mutants perform better when the binding rate is high or when the WT catalysis rate is low, such that catalysis is rate limiting. Our findings, however, suggest that drug binding is predominantly the rate-limiting step in the efflux of antibiotics from the cell.
In general, high regulatory noise results in heterogeneity in phenotypic resistance whereas low noise results in homogeneous resistance. By controlling the gene-expression noise, we show that noise reduction may facilitate bacterial inhibition of mutants by reducing their IC90. At the same time, mutants also have increased IC50, suggesting that drugs could be even less efficient when falling into sub-MIC concentrations. This is relevant for exploring the clinical potential of treatments which modulate gene-expression noise. Such noisemodulating chemicals have, for example, recently shown promising effects on reactivating HIV from latency, a process that relies on high gene expression noise (Dar et al. 2014). The success of HIV-latency modulators has provided a new concept in drug discovery and we envision that a similar approach may be tested to regulate resistance in bacteria.
In clinical antibiotic treatments, the drug concentration is not maintained at a consistently high level, as opposed to laboratory experiments. Here, we assume a constant drug scenario as a minimalist approach for two reasons: 1) computational efficiency as cells can be independently simulated; and 2) if the cells can modify the drug concentration, then we would also have to account for spatial variation of this concentration, which in turn would require biomechanical models of cell division and movement. The implementation of such a complex model would cloud the inferences we make from cellular stochasticity. On the other hand, the constantconcentration assumption is only valid when molecule numbers are high and when local fluctuations in the drug concentration are negligible. The latter may not apply if pharmacokinetics exist or if drugs are deactivated by the bacteria. One specific scenario that makes our assumption valid is if we assume the media is well-mixed: that is, the diffusion speed of the drug in the media is much faster than the diffusion across the cell membranes. If the volume of media is much greater than the total volume of all cells in the media, then the external drug concentration is approximately constant across the experiment. Further complications would be the release of drug upon cell lysis, which again would not be a problem in the aforementioned largevolume, high-density scenario.
While working with constant drug concentrations, we varied the duration of the drug dose to study how pulsed Sun et al. . doi:10.1093/molbev/msz199 MBE treatments affect bacterial survival. This is relevant as the mutant classes show different extinction times (supplementary fig. S6, Supplementary Material online), which could affect treatment success. Regulatory mutants have long survival times even at high concentrations, which leads to high survival probabilities as long as the drug is applied in a brief pulse. Thus, regulatory mutants could be more clinically problematic than structural mutants despite having lower IC90. This outcome will not be reflected in standard experimental protocols based on constant concentrations such as MIC measurements. Therefore, a comprehensive understanding of mutant dynamics requires more detailed assessment methods for cell death and growth across different antibiotic concentrations.
Our drug-target model is based on ciprofloxacin, gyrase, and the AcrAB-TolC efflux pump. This system has the advantage of relatively few drug targets per cell and drug that acts at low concentrations, which reduces computation time. We parametrized our model based on previous experimental findings. It can accurately reproduce cell population dynamics under various drug concentrations as well as within-cell protein distributions under stochastic expression ( fig. 4E and supplementary fig. S4B, Supplementary Material online). We therefore expect our model and its findings to be qualitatively robust.
Given that stochastic gene expression is a general mechanism in cell biology, our findings may also offer insights into general bacterial adaptation to changing environments. One example is the evolutionary reproducibility of adaptation via gain-of-function mutations: in certain biological systems, regulatory mutations seemingly take precedence over structural mutations and vice versa (Toprak et al. 2012;Blank et al. 2014;Lind et al. 2015). Here, we show that stochastic gene expression has temporal and heterogeneous effects on regulatory versus structural mutants, which helps explain why certain mutations could be advantageous.
In conclusion, the dynamics of gene expression could have a strong impact on the efficacy, and therefore the evolution, of resistance mutations. In particular, the intrinsic stochasticity of gene expression is a crucial determinant of evolutionary success. The framework we developed in this study now opens further opportunities to assess the impact of stochastic gene expression on bacterial adaptation, and the impact of molecular biology on evolution by extension.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online. All code and generated data have been archived at https://doi.org/10.5281/zenodo.3375506.