Dynamical model of antibiotic responses linking expression of resistance genes to metabolism explains emergence of heterogeneity during drug exposures

Abstract Antibiotic responses in bacteria are highly dynamic and heterogeneous, with sudden exposure of bacterial colonies to high drug doses resulting in the coexistence of recovered and arrested cells. The dynamics of the response is determined by regulatory circuits controlling the expression of resistance genes, which are in turn modulated by the drug’s action on cell growth and metabolism. Despite advances in understanding gene regulation at the molecular level, we still lack a framework to describe how feedback mechanisms resulting from the interdependence between expression of resistance and cell metabolism can amplify naturally occurring noise and create heterogeneity at the population level. To understand how this interplay affects cell survival upon exposure, we constructed a mathematical model of the dynamics of antibiotic responses that links metabolism and regulation of gene expression, based on the tetracycline resistance tet operon in E. coli. We use this model to interpret measurements of growth and expression of resistance in microfluidic experiments, both in single cells and in biofilms. We also implemented a stochastic model of the drug response, to show that exposure to high drug levels results in large variations of recovery times and heterogeneity at the population level. We show that stochasticity is important to determine how nutrient quality affects cell survival during exposure to high drug concentrations. A quantitative description of how microbes respond to antibiotics in dynamical environments is crucial to understand population-level behaviors such as biofilms and pathogenesis.


Introduction
Antibiotic responses in bacteria are remarkably dynamic and heterogeneous.Microbes typically carry regulated mechanisms of antibiotic resistance, which are repressed when the drug is absent to avoid the costs associated with expression of resistance genes [1,2].Therefore, when challenged by antibiotics, the cell has only a small window of time to activate its defenses before gene expression is halted by the drug action.Since expression of resistance genes is subjected to the strong stochasticity inherent to bacterial physiology, induction of the response is not always successful [3][4][5].Sudden exposures to antibiotics have been shown to result in phenotypic heterogeneity, diversity of response outcomes in single cells and complex growth patterns at the population level [6][7][8][9][10][11][12].To understand how microbial populations survive antibiotic treatments, we need models of antibiotic responses accounting for the dynamic and heterogeneous nature of antibiotic resistance.
Following antibiotic exposure, expression of resistance genes is controlled not only directly by regulatory mechanisms, but also indirectly by global effects of the drug on cell growth and gene expression, with a growing body of literature linking metabolism to the ability of bacteria to resist antibiotic treatments [13][14][15].Accumulation of resistance proteins in the cell interior depends on the resistance gene's expression rate and on the dilution of cell components due to cell growth, both of which are affected by the cell's metabolic environment.The presence of antibiotics itself alters the cell's growth dynamics and allocation of metabolic resources [16], which in turn affect expression and dilution of resistance, and consequently survival upon drug exposure.During drug responses, many antibiotics hamper the induction of resistance genes, which leads to further drug accumulation and decreased expression [17].These metabolism-mediated feedback mechanisms affect the course of antibiotic responses and can potentially amplify stochastic variations resulting in the coexistence of live and arrested cells [18].Stochastic generation of phenotypic diversity has also been studied in other systems, such as persistence, where a few cells from an isogenic population resist antibiotic exposures by growing slowly [19,20], or in heat stress responses in yeast [21].
Due to their small size, cellular components in bacteria are often present in small numbers and are subject to stochastic fluctuations.In particular, the transcription factors that regulate gene expression are often present on the order of tens of molecules per cell.Therefore, transcriptional regulation greatly increases stochasticity in gene expression, leading to strong variability in expression levels even among genetically identical cells [22,23].When this variability is high, it can be harnessed as a resistance strategy under stressful conditions [4,[24][25][26][27]. Bacterial colonies often harbor a subpopulation of cells with high levels of resistance proteins prior to drug exposure, increasing the chances of colony survival (heteroresistance) [28,29].However, even in the absence of strong pre-existing heterogeneity, subtle variations in the expression of resistance genes can be amplified during antibiotic responses by regulatory mechanisms controlling gene expression [30], diverging the course of responses among single cells, and ultimately leading to phenotypic heterogeneity at the population level.Therefore, isogenic populations growing under homogeneous environmental conditions still show remarkable diversity of outcomes at the singlecell level during antibiotic responses [24,31].
Here, we develop mathematical models of the dynamics of an antibiotic response, incorporating drug effects on cell growth and gene expression, to explain the emergence of heterogeneity during drug exposures.We start with a deterministic model that reproduces the progression of drug accumulation, expression of resistance and cell growth during a drug response.Then we develop a stochastic model to show how noisy dynamics can lead to phenotypic heterogeneity.We base our models on the classical E. coli tetracycline resistance tet operon, which displays many general characteristics of regulated antibiotic responses [26,32,33].The tet resistance mechanism is tightly regulated, controls a resistance gene that poses a significant cost for the cell, and is not directly regulated by other cellular processes [6,34,35].Cells carrying the native tet operon were shown to coexist in growing and non-growing states upon exposure to a large dose of tetracycline [6].The tet operon consists of two genes, an efflux pump TetA and its repressor TetR.In the absence of tetracycline, TetR represses both TetA expression and its own [36].TetR has a strong affinity for tetracycline, binding the drug as it enters the cell, which causes a conformational change resulting in loss of affinity for DNA and release of expression of both TetA and TetR.Efflux pump TetA is then rapidly produced, exporting tetracycline out of the cell.As the intracellular tetracycline concentration decreases, TetR resumes repression, avoiding a toxic overexpression of TetA [37].Since the tet mechanism directly senses the intracellular presence of the drug and elicits a fast and strong response, it is an ideal system to study antibiotic response dynamics and heterogeneity.
We describe the response dynamics of a widespread mechanism of regulated antibiotic response, which is broadly applicable to other responses where the cell needs to react quickly to rapidly changing environments [6,13,38], bringing the system out of equilibrium.Since the dynamics described in our model are not particular to antibiotic resistance mechanisms, many insights from this formulation can also be generalized to other transcriptionally repressed cellular responses that both sense and negatively act upon a chemical signal.

Dynamical model of the tetracycline response
To capture the dynamics of cell responses, we developed a mathematical model based on the main biochemical interactions involved in the E. coli tetracycline response (figure 1(a)).We integrate drug diffusion and accumulation into the cell, transcriptional regulation, and expression of resistance genes, as well as global effects of the drug on cell growth and gene expression.Tetracycline is a ribosome inhibitor [39], which reduces cell growth [40] and causes the cell to upregulate ribosome production in response [41].Altered ribosome levels then result in changes in the partition of the proteome [41,42], affecting expression of non-ribosomal genes including antibiotic resistance (described in detail below).Since in our model the cell growth rate is variable, a cell dies when the drug action causes ribosome function to cease, and the growth rate becomes zero.This system is largely governed by two feedback mechanisms: a stabilizing negative feedback provided by transcriptional repression [43][44][45][46][47], and a positive feedback mediated by the metabolism (growth rate-dependent) that can lead to bistability [18,48].
The model consists of a system of three differential equations that track changes in repressor TetR (r), efflux pump TetA (a), and intracellular drug (d) concentrations: Intracellular drug concentration changes over time according to three processes: the influx of drug from the environment into the cell, the export of drug by the efflux pump TetA, and the dilution of intracellular components due to cell growth.Since tetracycline is not actively degraded in the cytoplasm [49][50][51][52] and TetA and TetR are stable proteins [53,54], dilution due to cell growth is the main process driving down the concentration of these components [7] (except for efflux of tetracycline by TetA).Tetracycline enters the cell by diffusing through the cell membrane, with a rate proportional to the difference between the extracellular and intracellular drug concentrations (D (t) and d f ) with a diffusion constant K i .Although diffusion through the hydrophobic membrane is relatively slow, with a half-equilibration time around 45 min [54,55], intracellular tetracycline is already toxic for the cell at low concentrations in the nanomolar scale.Therefore, exposures to extracellular tetracycline in the micromolar scale, which tetequipped E. coli can resist, still results in intracellular drug rapidly reaching toxic levels.Efflux pump TetA exports tetracycline out of the cell efficiently, following standard Michaelis-Menten kinetics described by , where K a is the catalytic rate constant, and k a the Michaelis constant.In growing cells, the intracellular nanomolar concentrations of tetracycline do not significantly saturate TetA.As the cell grows, the drug is diluted in the cell interior, with the dilution rate equal to the growth rate λd.The cell growth rate is not fixed, and depends on drug action and metabolism, as detailed below.
In the cytoplasm, tetracycline strongly binds repressor TetR, which then undergoes a conformational change and loses capacity to bind its DNA binding site [6,21].Since biochemical binding and unbinding reactions happen at much faster timescales than the other relevant processes, we consider a chemical equilibrium d f + r f ⇌ r b between the unbound (free) forms of drug and TetR (d f , r f ) and the bound form r b , with an equilibrium constant K d .This equilibrium results in and , and the concentration of free TetR can be calculated by solving the resulting quadratic equation While the bound form of TetR is inactive, the free form transcriptionally regulates expression of both TetA and TetR.TetA concentration changes over time according to its synthesis f H a ( r f ) and its dilution due to cell growth λa.We model TetR regulation of TetA synthesis using a Hill function , which describes the equilibrium binding of free repressor TetR and its binding sites in the promoter region of TetA.A is the fully induced expression rate, in the absence of TetR repression.Free repressor r f decreases TetA expression, with r 0,a being the repressor level for half-maximal expression.The Hill coefficient n is a measure of how sharply expression rates transition between high and low levels around the threshold r 0,a , and is related to the cooperativity in TetR DNA binding.We use a Hill coefficient of n = 4 to reflect TetR binding two different DNA binding sites as a dimer [36,56] (SI).The factor f (λ, d) modulates the expression of resistance proteins according to global effects of drug action on cell metabolism (cell growth) and is detailed below.TetR concentration is similarly determined by its synthesis f H r ( r f ) and its dilution due to growth λr.Even in the absence of drug, gene expression depends on the cell growth rate, which is set by the quality of the nutritional composition of the immediate environment.The cell grows according to the total output of translation, and therefore the growth rate is proportional to the ribosomal content of the cell.Faster growing cells then harbor larger ribosomal content, and consequently less nonribosomal proteins.Additionally, since tetracycline is a ribosome inhibitor, it also causes the cell to upregulate ribosomal content, decreasing non-ribosomal content.A fraction ϕ Q ≈ 45% of the proteome, thought to consist of housekeeping genes, rescales gene expression proportionally to the growth rate, balancing synthesis and dilution in order to maintain constant expression levels [41].To accommodate for changes in ribosomal content, another variable fraction ϕ P of the proteome adjusts its expression level in the opposite direction.In our model, we account for the effects of intracellular drug concentration and nutrient quality on gene expression and cell growth by incorporating a proteome partitioning framework developed by Scott et al [41].According to this framework, all proteins in the cell fall within one of three categories: (1) Sector Q proteins whose concentrations are not affected by changes in ribosomal content, caused by changes in cell growth or translational inhibition by drug action (e.g.housekeeping genes), ( 2) Sector R of ribosome-affiliated proteins, whose concentration increases with the growth rate or translational inhibition, and (3) Sector P of all other proteins, whose concentration decreases with growth rate or translational inhibition to compensate for the increase in ribosomal proteins.
As the fraction ϕ Q of the cell's proteome that consists of proteins not affected by translational inhibition does not change, the allocation of resources towards ribosome-affiliated proteins and all others adds to a fixed portion of the proteome, with the two sectors varying according to the cell's translational (κ t ) and nutritional (κ n ) capacities.κ t relates to ribosomal function, measured by the global rate of translation elongation, and κ n relates to nutrient quality, or the capacity of the culture medium to support growth.Here, we refer to the base value of the cell's translational capacity in drug-free medium as κ 0 t , which has a universal value of 4.5 h −1 for E. coli [41] (henceforth, we use subscript or superscript 0 to indicate quantities under full-growth drug-free conditions).The translational capacity κ t is reduced by the presence of intracellular drug, which binds and inactivates ribosomes, and can thus be modeled by , where K ribo is the dissociation constant for drug-ribosome binding.Lower K ribo values correspond to stronger inhibition, and vice versa.The nutritional capacity κ n reflects the nutrient quality of the medium, where media with higher κ n allow faster growth rates.Here, we focus our analysis on situations where nutrient quality does not change, such as in exponential growth in liquid cultures or stable growth in microfluidic devices, where κ n can be determined as a fixed value that fits the growth rate allowed by the culture medium under drug-free conditions (see the SI).Otherwise, decreases in κ n can be calculated to reflect nutrient consumption, such as in saturating liquid cultures or as in the nutrient gradients generated by spatial structure in biofilms [7].
According to the framework of the proteome partition, the fraction of the proteome ϕ R dedicated to ribosomal proteins varies linearly with the cell growth rate, from ϕ R = ϕ min R when λ = 0 to a theoretical maximum ϕ R = ϕ max R ≈ 55%, since the remaining ϕ Q ≈ 45% of the proteome is fixed and not affected by changes in metabolism [41].Therefore, the difference between maximal and minimal ribosomal fraction is the variable part of the proteome that can be occupied by the sector P proteins that are affected by translation inhibition and nutrient quality.The fractions of the proteome dedicated to ribosomal sector R and variable sector P proteins are given by ϕ The cell growth rate is proportional to the product of ribosomal content and translation capacity, and is calculated as , where ρ = 0.76 is the conversion factor between RNA/protein ratio and ribosomal fraction calculated for E. coli.
Since expression of TetA and TetR has been shown to depend on cell growth, we assume these proteins belong to the variable sector P. Therefore, without regulation of the synthesis rates H a ( r f ) and H r ( r f ) , we would expect the concentrations of TetA and TetR to scale by ϕ P .We can then calculate the dependence of TetA and TetR synthesis on both translation inhibition and nutrient levels.This dependance is composed of two factors, one reflecting changes in growth rate and another reflecting changes in proteome partition.(λ/λ 0 ) scales down gene expression to match the decrease in growth rate, which does not change expression levels.For proteins in the P sector, there is further modulation by a factor , where ϕ 0 P is the P-sector fraction under full nutrients and no drug.Therefore, the global metabolic effects on the expression rates of P-sector genes can be modeled by f = (λ/λ 0 ) This leads to expression level steady states proportional to ∼ κ t / (κ t + κ n ), as expected for proteins in the P sector (see ϕ P dependence above) and simplifies to f = 1 at full nutrients and no drug.We note here that although reducing either the nutritional or translational capacities result in decreased growth rates, they have opposite effects on the proteome partition.Therefore, while reducing the growth rate by nutrient limitation results in ribosome downregulation and increased TetA and TetR expression levels, reducing the growth rate by tetracycline exposure results in ribosome upregulation and decreased TetA and TetR levels.Next, we use our model to simulate the tetracycline response under different nutrient conditions and drug concentrations.

Dynamical model captures the dynamics of cell responses
We begin by analyzing the time course of a typical response to a sudden exposure to tetracycline (figures 1(b)-(g)).Initially, as TetA levels are low, the drug quickly diffuses into the cell and accumulates in the cytoplasm, reducing the cell growth rate and slowing down gene expression.As the incoming drug quickly binds and inactivates repressor TetR, expression of both TetR and TetA is released shortly after exposure, although initial accumulation is slow.As TetA levels begin to increase in the cell membrane and export tetracycline back out of the cell, drug accumulation in the cytoplasm slows down and eventually reverses.As intracellular drug levels decline, cell growth and gene expression recover, accelerating TetA accumulation.When intracellular drug returns to low levels, TetR is released, resuming repression.TetA, TetR and tetracycline then equilibrate to steady-state levels, which depend on both TetR regulation and on the effects of the proteome partition on gene expression.This time course of the response dynamics qualitatively reproduces experimental measurements of expression of resistance genes and cell growth (figures 1(c)-(f)), obtained in a resistant population of E. coli cells suddenly exposed to a large dose of tetracycline in a microfluidic device [6,57] (data obtained from [6]).
To understand how tetracycline concentration affects cell growth and survival, we simulated the drug response to exposures of increasing drug doses.As drug concentration increases, the growth reduction experienced by the cell at the beginning of a response also increases in both duration and magnitude (figures 1(h)-(j)).At high drug concentrations, synthesis of TetR and TetA is also significantly reduced during this state of translational inhibition, bringing the cell to a quasi-arrested slow growth state.This state can be escaped if the cell maintains cell growth, however slow, to eventually accumulate TetA to sufficient levels to kickstart drug export (figure 1(i)).With enough TetA, the cell enters a positive feedback loop where reduced intracellular drug concentrations leads to higher growth rates and faster TetA production, resulting in further reduction of intracellular drug levels and faster cell growth.As the drug dose increases further, the cell is trapped in the slow growth state for increasingly longer times before eventually recovering.At very high drug doses, however, recovery is no longer possible, with the cell reaching sufficient translational inhibition such that TetA is not produced at high enough rates to initiate a recovery (figure 1(j)).The cell then cannot counteract the influx of drug, causing the growth rate to be further reduced by the same positive feedback mechanism.
Increases in drug concentration result in longer recovery times, defined as the time it takes for the cell to recover to the average of its minimum and final growth rates, up to a threshold drug concentration D thr where the recovery time tends to infinity (figure 1(k)).Past this threshold, the cell is permanently arrested following exposure and does not recover.We next determined the stable growth rates at the end of the antibiotic response for different doses of tetracycline, to examine the effect of the drug during steady state.Increasing drug doses cause a reduction in steady-state growth, while still keeping it at relatively high levels up to the D thr threshold drug dose (figure 1(l)).At higher drug doses, steady-state growth is low, corresponding to the arrested state of cells that do not recover.Therefore, the outcome of the tetracycline response is binary, with cells either permanently arrested or recovering to relatively high growth rates, as also observed experimentally.Slowgrowing cells are observed only in the transient following drug exposure, and eventually resolve into either arrest or recovery.
Cell survival to drug exposure also depends on the nutrient condition of the culture medium.To determine how growth conditions affect the outcome of antibiotic responses, we simulated the tetracycline response for decreasing values of the nutritional capacity, emulating the nutrient gradients commonly found in structured microbial communities (figures 2(a)-(d), reproduced from [7]).Decreases in nutritional capacity result in lower growth, both in the presence and in the absence of drug.However, cells growing under lower nutritional capacity can recover growth at higher drug doses, showing a significantly higher D thr and reduced recovery times (figures 2(e)-(j)).Therefore, slow-growing cells in poor media show higher resistance to drug exposures, suggesting a trade-off between cell growth and drug resistance also found in other studies [19,48,58].This result explains cell survival in biofilm microfluidic experiments (figures 2(a)-(d)), where spatial structure causes nutrient gradients decreasing from the surface to the interior of the colony.In such experiments, sudden exposure to tetracycline results in the arrest of fast-growing cells at the surface of the colony, while slow-growing cells in the interior survive exposure and grow to regenerate the surface.Taking the effects of the drug on the proteome partition into account is essential to explain this advantage of slowgrowing cells (compare figures 2(e) and S1 and the accompanying text in the SI).
Taken together, these results suggests that at high drug doses close to the threshold concentration D thr , small fluctuations could decide the cell's fate.Since drug responses resolve into either growth recovery or arrest, and both states exist within small margins of drug concentrations and nutrient conditions, we hypothesize that noisy dynamics can lead to the coexistence of live and arrested cells within the same population [18,21].Therefore, to study the emergence of heterogeneity during the tetracycline response, we next develop a stochastic model of the response dynamics, which incorporates noise as an integral part of its formulation.

Stochastic model of the tetracycline response
To simulate the generation of heterogeneity during antibiotic exposures, we developed a stochastic model of the response where the main biochemical reactions are considered as Poisson processes with known rates.The reactions, summarized in table 1, correspond to the same processes as described by the deterministic model, with drug efflux, transcriptional regulation, and drug action by ribosome binding implemented explicitly.TetA, TetR, intracellular tetracycline, ribosomes and TetR binding sites, both in free and bound forms, are considered as discrete quantities, resulting in a discrete system with known transition rates between states.Stochastic models incorporate noise into their formulation and can be simulated numerically to obtain probability distributions of possible outcomes of the response.
Tetracycline passively diffuses into and out of the cell through the membrane with the same rate K i , resulting in a net diffusion as given by the deterministic model.Inside the cell, tetracycline binds TetR reversibly into an inactive compound.Drug efflux follows Michaelis-Menten dynamics, starting with a reversible binding of tetracycline to TetA, followed by transport of the drug to the extracellular medium.Protein synthesis is combined into a single reaction, with a rate that depends on ribosome availability (detailed below).Repressor TetR regulates protein synthesis by binding to its DNA binding sites and blocking transcription initiation, and therefore TetA and TetR syntheses are modulated by the occupancy of the binding sites.TetR binds two sites O 1 and O 2 as a dimer in the promoter region of the tet operon (figure S2 and the accompanying text in the SI).We do not explicitly model the dimerization of TetR [36], but instead consider that the binding rate of TetR binding sites depend on the square of the concentration of TetR monomers, reflecting the binding of TetR dimers [24,25,59].TetA is synthesized from a single strong promoter P A , which is active only when both binding sites are free.TetR is synthesized from two independent promoters of similar strength: P R1 is active when both binding sites are free, and P R2 is active when O 1 is free.Dilution of cellular components is modeled as a death process, with a rate equal to the cell growth rate λ.The reaction rates K i , K a , R, A, and λ are the same as in the deterministic model (as summarized in table S1).The forward and backwards rates for the equilibrium constants K d , K ribo , K b and k a were chosen based on typical values for the timescale of each reaction: binding of small molecules to proteins is in the order of milliseconds in bacteria, while binding of transcription factors to DNA is in the order of seconds.
To incorporate the effects of drug action and nutrient quality, we explicitly consider a ribosome pool of variable size N R , which changes according to the nutritional and translational capacities of the cell.We set a theoretical maximum size of the pool N max R corresponding to the maximum ribosomal content ϕ max R .While it is impractical to set N max R to typical numbers of ribosomes per cell (tens of thousands), a reduced pool still captures the effects of ribosome binding by the drug on protein synthesis, cell growth and proteome partition.Then, the proteome partition is calculated from the size of the ribosome pool, with Correspondingly, the P sector fraction can be calculated from is the minimum ribosome content when growth is zero, corresponding to ϕ min R , and can be obtained from We can therefore calculate the cell growth rate as Tetracycline binds and inactivates ribosome, with an equilibrium constant calculated from the drug concentration necessary for half-repression of cell growth, determined experimentally.The reduction in translation capacity then reflects the inactivation of a fraction of the ribosome pool and is calculated as κ t = κ 0 t N Rf /N R , with N Rf being the number of free ribosomes.
The theory of proteome partition assumes that ribosome production is regulated to achieve optimal levels in each nutrient and drug condition, with the resulting proteome partition affecting the expression Dilution due to cell growth of proteins in the P sector.Therefore, we take the approach of calculating the ribosome production rate based on the translational and nutritional capacities κ t and κ n , and then using the resulting changes in the size and availability of the ribosome pool to calculate the effects of drug action in the growth rate and expression of TetA and TetR.In the absence of drug, the average size of the ribosome pool is determined by the nutritional capacity of the medium and the maximal translational capacity as . Therefore, we define a basal rate of ribosome production R 0 R = λ 0 N 0 R to generate a pool of the expected size.Finally, we calculate the dependencies g = (λ/λ 0 ) that modulate the synthesis of ribosomes and TetA/TetR under drug exposure, respectively.These dependencies have a common factor , which adjusts gene expression to the reduced growth rate without changing expression levels.Additionally, ribosome synthesis is modulated by ϕ , reflecting the upregulation of ribosomes when translation is compromised, and TetA/TetR synthesis is modulated by reflecting changes in the proteome partition.We note that the factor λ/λ 0 that affects ribosome production is calculated with the growth rate extracted from the pool itself, which is also used to calculate the dilution rate, and therefore does not affect the expected pool size.Since we calculate ϕ R /ϕ 0 R using κ t and κ n instead of fractions calculated from the pool, the expected pool size does not depend on itself.However, ϕ P /ϕ 0 P is calculated from the pool, so that gene expression responds to changes in the pool size.

Stochastic model shows large variations in response dynamics
We started by simulating the system during a drug response to moderate drug concentrations, first in the absence of drug for a period of 1000 min, and then introduced an extracellular drug concentration of 100 µM (figures 3(a)-(e)).The stochastic model closely follows the TetA, TetR and intracellular drug trajectories observed in the deterministic model following drug exposure, with expression rates depending on transient changes in partition fractions and ribosomes.Both in the absence and presence of drug, the size of the ribosome pool equilibrates to the expected steady-state levels.The high influx of tetracycline into the cell following exposure is counteracted by a temporary upregulation of ribosome production, which maintains protein expression, but reduces the resources available to P-sector proteins (TetA and TetR).
The full model, despite being exact, is computationally costly due to the high number of reactions being considered, with a single trajectory taking up to 24 h to simulate.To generate a large number of trajectories efficiently, we identify the fastest processes in the system (figure S4), which contribute less noise, and approximate their effect using the adiabatic limit where these processes can be considered as being in equilibrium.We find that most computational power is spent on fast binding/unbinding reactions of tetracycline with its ligands.Eliminating explicit tetracycline binding greatly improves computational speed while keeping the most significant sources of stochasticity.The system then shows stochastic noise coming only from the slower synthesis/degradation, drug import/export and DNA binding reactions.Free tetracycline molecules reversibly bind TetR, TetA and ribosome, with equilibrium constants K d , k a and K ribo , respectively.Therefore, the copy number of the free and bound forms can be estimated from the total number, such that is the number of bound TetA molecules and Ribo is the number of free ribosomes.Since the export of tetracycline out of the cell is proportional to the number of bound TetA complexes, we arrive at an export rate of K a a T d f / ( k a + d f ) as in the deterministic model.Finally, we calculate the number of free tetracycline molecules from where b denotes the bound forms.This equation is solved numerically to obtain d f at each time step in the simulation.
In this reduced stochastic model, only total copy numbers of the chemical species are modeled explicitly, while free and bound forms are calculated from the equilibrium constants and used to calculate the remaining reaction rates.Simulations of the reduced model had a 40-fold faster runtime (∼0.6 h per trajectory), and the overall response dynamics still closely followed the full model, displaying strong variability of recovery times.We tested ribosome pool sizes of N max R = 100 (figures 3(f)-(j)) and N max R = 1000 (figures 3(k)-(o)).Although the reduction of the ribosome pool introduces noise into the system, drug responses show qualitatively similar behavior in either case.The size of the ribosome pool is therefore a key parameter measuring the amount of stochasticity introduced into the system by the drug action on cell growth and gene expression.These results suggest that such adiabatic approximations are valid, and the reduced model can be used to generate probability distributions of response outcomes.

Exposure to high drug doses results in coexistence of growing and arrested cells
We now use the reduced model to simulate a large number of trajectories of the system to obtain a probabilistic distribution of the system's behavior.For each set of parameters, we simulated 1000 trajectories, which allows us to determine the shape of the distribution of recovery times and measure the probability of survival to drug exposures within a 5% tolerance.In each simulation, the system starts with the expected concentrations in the absence of drug, calculated from the deterministic model, and is simulated for 1000 min to reach equilibrium before drug is introduced.The system is then simulated in the presence of drug until steady-state dynamics are reached.In multistable systems, stochasticity can make the response reach different stable points.Therefore, trajectories starting from the same initial conditions can take increasingly different paths, resulting in either recovery or arrest and showing widely different recovery times.Within an isogenic microbial population in a homogeneous environment, different trajectories can be seen as the fates of individual cells, with the probability distribution of response outcomes reflecting the generation of heterogeneity in the population during drug responses.For each drug and nutrient condition used in the simulations, we calculate the probability that responses will result in growth recovery, and, from the subset of recovery trajectories, we calculate the distribution of recovery times.
Simulations of the stochastic model during exposure to high drug doses result in large variations in intracellular drug accumulation, eventually leading to the coexistence of recovered and arrested trajectories (cells).This regime, which includes short recoveries and cell arrest (infinitely long recoveries), corresponds to the very long recovery times observed in deterministic simulations when the drug concentration approaches D thr (figure 4(a)).While in some trajectories the cell was able to curb the influx of drug relatively quickly, in others intracellular drug reached much higher levels.Recovery times strongly depend on the maximum drug levels reached inside the cell.When intracellular drug reaches high levels, expression of TetA is compromised, and the cell becomes trapped in a slow-growth state.Taken together, these results suggest the existence of a semi-stable lowgrowth state, from where the system can escape to recovery if TetA concentrations are high enough.When intracellular drug levels are kept to lower levels, recovery is much quicker.Therefore, the initial dynamics of the response is subject to strong fluctuations, and it is crucial to determine cell fate.
At very high drug doses, in regimes where deterministic simulations no longer recover, stochastic simulations may still result in growth recovery, although they become increasingly rare and recovery times tend towards higher values.For exposures to drug concentrations past D thr , we no longer see an increase in recovery times, with distributions still centered at ∼7 h.However, larger drug levels result in decreased survival, and at very high drug concentrations no cells survive (figure 4(b)).These results agree with the experimental observation that cell fate is decided within the first few hours after exposure, after which recovery is rare.
Next, we investigated how nutrient conditions affect the dynamics of the antibiotic response during an exposure to a moderate drug concentration.For high nutritional capacities, average recovery times are faster than those predicted by the deterministic model (figure 4(c)).Since recovery depends on TetA levels reaching high enough levels to reverse the influx of drug, stochasticity in TetA expression increases the probability that this threshold is reached earlier.However, unlike in the deterministic model, lowering the nutritional capacity does not result in shorter recoveries (figure 4(d)).Although the probability of recovery decreases with lower nutrient quality, recovery times remain similar among cells that do recover.Overall, since recovery from drug exposure is essentially a process of threshold crossing, stochasticity in the expression of resistance results in faster recovery than what is predicted by the deterministic model under most growth conditions.

Discussion
Our model of the classical tetracycline response in E. coli combines biochemical reaction dynamics with gene regulation and global effects of the drug on metabolism to provide a framework to understand microbial antibiotic resistance in the context of dynamic and heterogeneous populations.Upon drug exposure, expression of resistance does not only depend on direct regulation by transcription factors but also on global effects on protein expression linked to the metabolic state of the cell [17,42,58].And since restoration of metabolic functions (i.e.cell growth) also depends on the expression of resistance, this metabolism-mediated link introduces an additional feedback mechanism in the control of drug responses: failure to quickly deploy resistance genes results in higher levels of intracellular drug, leading to further reduction in expression of resistance and further drug accumulation.This positive feedback mediated by metabolic effects can act as a switch and is known to generate bistability and the coexistence of growing and arrested cells in the presence of antibiotics [18,30,48,60].We note that this model does not include cell death explicitly, which is consistent with tetracycline being a bacteriostatic drug.Rather, when the intracellular drug concentration is too high, the cell is permanently arrested in a state with no cell growth or gene expression, which can be reversed if the drug challenge is relieved.Therefore, by considering the interplay between drug action, cell growth and expression of resistance, we are able to build a comprehensive model that explains the progression of antibiotic responses in single cells and, at the same time, captures the emergence of heterogeneity in the populations and colony-level collective behaviors.
Our models were able to recapitulate the complex dynamics observed experimentally in two different microfluidic experiments.A single-cell microfluidic experiment measuring phenotypic diversity during antibiotic responses found the emergence of remarkable diversity of growth and gene expression within isogenic populations of tet resistant E. coli strains [6].These experiments observed the coexistence of three phenotypes: fast-growing recovered cells, arrested cells with little TetA expression, and temporary moribund cells that grow slowly and eventually either recover or arrest.Our deterministic model was able to generate all three phenotypes observed in singlecell responses, including the existence of the semistable low-growth moribund state, from where the system can escape to recovery if TetA concentrations are high enough.Simulations of our model suggest the existence of distinct stable phenotypic states, suggesting future analysis to determine the exact nature of these phenotypes, and particularly what differentiates the intriguing moribund state from complete arrest.As is already done in our stochastic model, our deterministic model can also be extended to include the dynamics of the proteome partition, which are currently considered to be in equilibrium.
Our stochastic model recapitulates the emergence of phenotypic heterogeneity during drug responses, which is observed even in isogenic populations in homogeneous environments, explaining the coexistence of stable states corresponding to recovered and arrested cells.Gene expression in bacteria is known to be noisy, with large cell-to-cell variation, particularly in regulated genes.Stochastic models capture this variability, describing how naturally occurring noise in cellular processes is propagated through mechanisms of regulation to generate different outcomes during antibiotic responses.As observed experimentally, our simulations find that faltering expression of resistance upon exposure results in slow growth and delayed recovery, with large variation of recovery times.We find that at large drug concentrations, while the majority of the microbial population is arrested upon exposure, small subpopulations can still survive and regenerate the population.Although stochasticity during the induction of antibiotic responses also generates phenotypic diversity in the population, with the coexistence of live and arrested cells, this mechanism is distinct from persistence.While in the mechanism we describe phenotypic heterogeneity emerges after drug exposures, with surviving cells growing in the presence of the drug, persistence relies on pre-existing phenotypic variability, and dormant persistent cells that survive drug exposures recover growth only after the drug is removed from the medium [61].Moreover, the model can quantitatively predict how environmental factors such as the drug dose or the nutritional quality of the medium determine the distribution of cell states.Therefore, our models can be used to determine the patterns of population-level growth and expression of resistance, which result from the sum of all single-cell phenotypes.
The framework developed with this specific model can be used as the basis to develop large-scale stochastic simulations of whole drug-resistance pathways, which will be useful in identifying new drug targets.However, while our full stochastic model is comprehensive, its simulations are slow due to the necessity of serially executing reactions spanning a wide range of timescales.Therefore, the model is limited in the number of molecules that can be simulated in tractable time, especially with regards to the ribosome pool.Future implementations of this model could use more sophisticated versions of the Gillespie algorithm such as tau leaping to reintroduce stochasticity from fast reactions that is lost with the adiabatic approximations [62,63] or population dynamics algorithms [64].Increasing the efficiency of the simulations allows the consideration of larger systems with more realistic numbers of components, resulting in more accurate quantitative predictions.
A second microfluidic experiment [7] characterized how expression of resistance is coordinated during antibiotic responses across structured microbial populations (biofilms), where consumption of nutrients by cells closer to the surface generates chemical gradients towards the interior of the colony [7,[65][66][67].This experiment found that the quick arrest of fast-growing cells near the surface of the biofilm causes a redistribution of nutrients towards the interior and reactivates previously dormant cells at deeper layers (figure S5).Because lower metabolism increases resistance levels, reactivated dormant cells can survive exposure and repopulate the biofilm, maintaining population-level growth even at drug concentrations high enough to kill whole planktonic populations.Interestingly, the same three phenotypes from the single-cell studies were observed: fastgrowing cells with moderate TetA expression close to the surface of the biofilm, slow-growing cells with high TetA expression further into the biofilm, and arrested cells with little TetA expression in the interior.To study biofilms, our model can be extended with the addition of a spatial component, tracking nutrient concentration in a continuous field, as well as the introduction of nutrient consumption [7].In another study, a spatial representation of our model was able to recapitulate the dynamics of spatially heterogeneous growth patterns and expression of resistance-most notably the reactivation of dormant cells upon drug exposure and the persistent accumulation of resistance in the dormant subpopulation.
We find that slow-growing cells express higher levels of resistance genes and are better suited to survive sudden exposures.The dynamics of drug responses is characterized by markedly different gene expression levels between fast-and slowgrowing cells.Even in the absence of drug, TetA and TetR expression were experimentally determined to decrease linearly with the cell growth rate, which is predicted by the theory of proteome partition (figures 2(c) and (d)).In poorer media, slowgrowing cells are able to grow in the presence of the drug because of their higher expression of resistance.Therefore, the nutrient condition of the medium surrounding the cell affects its resistance profile.These results are in line with recent studies showing that competition for resources within polymicrobial communities affects both susceptibility to antibiotics and the genetic evolution of resistance [68].The effects of nutrient limitation are particularly relevant for biofilms, where nutrient gradients from the surface towards the interior of the colony generate an array of metabolic states that underlie collective resistance [69][70][71][72][73][74][75][76][77].
The framework developed here will be useful in investigating other important questions regarding antibiotic responses.While the simulations themselves only cover time scales up to one to two days, an analysis of the sensitivity of the system to parameter changes could also give insights into the types of mutations that can be expected from the evolution of antibiotic resistance mechanisms under regimes of repeated drug exposures on longer time scales [78][79][80][81].Apart from refining resistance proteins to be effective against specific antibiotics, such mutations could also target regulation to improve the responsiveness of the resistance mechanism.Our model can also be extended to analyze bactericidal drugs (tetracycline itself has bactericidal effects at high doses [82]).This could be modeled by the introduction of a probability of cell death or permanent cell damage, depending on intracellular drug levels.Considering bactericidal effects could lead to interesting new dynamics of the system, which will be important to explain the nature of slow-growing phenotypes and the recovery from drug exposures.
Microbial drug resistance in the real world is often at odds with lab measurements, with infections often returning from remission.Our model bridges the gap between the dynamics of drug responses at the singlecell level and the resulting collective behavior of the population, and helps to understand how subpopulations of microbial cells are able to resist exposure to high drug doses and regenerate the colony.A quantitative description of how cell responses are regulated in complex environments is crucial to understand community-level behaviors such as antibiotic resistance, pathogenesis, and biofilms, which often can be explained without invoking additional specialized mechanisms.

Simulation of the deterministic dynamical model
We numerically integrate the system of differential equations using the ode113 function in MATLAB [83], using values of external drug concentration and nutrient quality as input parameters.The parameter values used in the simulations are summarized in table S1 and were either estimated from our experimental data [6], or obtained from literature [26,84].The code for the simulations is available at our GitHub page: https://github.com/schultz-lab/Phys-Biol-2023.

Simulation of the stochastic model
We simulate this system using the classical Gillespie's stochastic simulation algorithm [85], which generates a large number of independent trajectories that are used to calculate probability distributions of the system states over time.Instead of using fixed time steps, the Gillespie algorithm calculates two probability distributions: one for the time that is needed for the next reaction event to occur and a second distribution characterizing which of the possible reactions will occur next.By choosing two random numbers from a random number generator, a value from the time distribution and a reaction are chosen.The time is then increased, and the system is updated accordingly.Each realization of this algorithm represents one trajectory.Transient time-dependent probability distributions can be obtained from a sufficiently large number of trajectories with different random seeds.

Figure 1 .
Figure 1.Drug concentration affects recovery from antibiotic exposures.(a) The tet resistance mechanism.Tetracycline (Tc) diffuses across the cell membrane, and binds repressor TetR, thereby releasing expression of TetA and TetR.TetA then exports the drug outside the cell.Tetracycline also binds ribosomes, inhibiting protein translation.Solid arrows represent protein synthesis, solid bi-directional arrows represent binding/unbinding, blunted arrows represent repression, dashed arrows represent import/export of tetracycline, and dotted arrows represent ribosome content promoting protein synthesis.(b) Schematic of the microfluidic device used for measuring growth and expression of resistance in single cells during drug exposures.The experimental data in this figure is reproduced from [6].(c) Average tetracycline, TetA, and TetR levels over the course of a response, measured in 40 tetracycline-resistant single E. coli cells during a microfluidic experiment.(d) Tetracycline, TetA and TetR levels over the course of a response as predicted by our deterministic model.(e) Three examples of time courses of single cells, resulting in different cell fates.Red and green colors represent expression of TetA and TetR, respectively, measured with fluorescent reporters.(f) Cell growth in single cells following a sudden tetracycline exposure (thin lines).Yellow, red, and black lines correspond to recovered, moribund, and arrested cells, respectively.Thick lines represent average growth within each group.Moribund cells that grow slowly following drug exposure eventually either recover or become arrested.(g) Varying the dissociation constant for drug-ribosome binding, the model reproduces the responses seen in single cells.(h)-(j) Growth rate, TetA, and TetR levels over the course of a tetracycline response, calculated at different extracellular drug concentrations.Higher drug concentrations lead to larger decreases in growth rate, eventually leading to arrest.(k) Recovery times increase with extracellular drug levels, approaching a vertical asymptote at a threshold drug concentration D thr .(l) Growth rate at the end of the tetracycline response for different extracellular drug concentrations.The final growth rate drops sharply around the threshold drug concentration D thr .

Figure 2 .
Figure 2. Slow-growing cells survive exposures to higher drug concentrations compared to fast-growing cells.(a) Schematic of the microfluidic device used for measuring growth and resistance expression in a bacterial colony.(b) Expression of resistance genes over time across the colony.The nutrient supply channel (not shown) runs along the top of the trap and flushes away any cells growing out of the trap.The experimental data in this figure is reproduced from [7].(c) Our tetracycline resistance model, which considers proteome partitioning, predicts a linear decrease in the expression of both resistance proteins with growth rate.(d) Our experimental data also measured a linear reduction in TetA and TetR levels with the growth rate in the absence of tetracycline.(e) Final growth rate at the end of the response for different combinations of extracellular drug concentration and nutritional capacity.Lower nutritional quality allows resistance to higher drug concentrations.(f) Threshold drug concentration D thr is higher for lower nutritional capacities.(g) Final growth rate varying with extracellular drug concentrations, for different nutritional capacities.(h)-(i) Growth rate over the course of a tetracycline response for different nutritional capacities, shown at two drug concentrations.(j) Recovery times varying with extracellular drug concentrations, for different nutritional capacities.Cells with lower nutritional capacities have shorter recovery times and maintain growth up to higher drug concentrations.

Figure 3 .
Figure 3. Stochasticity during antibiotic responses to moderate drug concentrations.Time courses of responses to a sudden drug exposure at time 0. Vertical dashed lines indicate drug exposure.(a)-(e) Full stochastic model, with tetracycline binding reactions considered explicitly, and a pool of 100 ribosomes.(a) Cell growth during the drug response.(b) Intracellular tetracycline, TetR, and TetA concentrations, normalized to maximum values.Thin lines show different trajectories obtained from the stochastic model.Dashed lines indicate the corresponding concentrations obtained from the deterministic model.(c) Ribosome pool, showing total and free ribosomes.(d) P-and R-sector fractions of the proteome.(e) Distribution of recovery times obtained from simulating ∼1000 trajectories.(f)-(j) Reduced stochastic model, considering tetracycline binding reactions in chemical equilibrium, and a pool of 100 ribosomes.This approximation maintains similar qualitative behavior and dynamics while being simulated in a 40-fold shorter runtime.(k)-(o) Reduced stochastic model and a pool of 1000 ribosomes.An increased ribosome pool results in decreased noise and a narrower distribution of recovery times.

Figure 4 .
Figure 4. Drug concentration and nutrient quality affect survival rates and recovery times upon exposure to high concentrations of tetracycline.(a) Distribution of recovery times among recovered trajectories, obtained from 1000 simulations for each drug concentration.The mean recovery time increases with extracellular drug concentration up to ∼7 h, obtained when drug concentration reaches the threshold D thr .Increasing drug further does not change recovery times, but results in decreased survival.Inset: Recovery times obtained from the deterministic model, with dashed lines representing the drug concentrations used in the stochastic simulations.(b) Probability of cell survival for different drug concentrations.Survival probability decreases with the drug dose and falls sharply for drug concentrations past D thr .(c) Growth rate in different trajectories during an exposure to 125 µM tetracycline (dashed vertical line).Stochastic fluctuations allow TetA to reach threshold levels necessary for recovery faster than in the deterministic model (dashed trajectory), resulting in faster recovery.Inset: Probability of survival for this set of simulations.(d) Recovery times for varying nutritional capacities, in the deterministic and stochastic models.The color gradient shows the probability of survival in simulations of the stochastic model.Inset: Recovery times obtained from the deterministic model for different nutritional capacities.The dashed line represents the 125 µM tetracycline concentration used in the stochastic simulations.

Table 1 .
Reactions of the full stochastic model.