Tunable Stochastic Pulsing in the Escherichia coli Multiple Antibiotic Resistance Network from Interlinked Positive and Negative Feedback Loops

Cells live in uncertain, dynamic environments and have many mechanisms for sensing and responding to changes in their surroundings. However, sudden fluctuations in the environment can be catastrophic to a population if it relies solely on sensory responses, which have a delay associated with them. Cells can reconcile these effects by using a tunable stochastic response, where in the absence of a stressor they create phenotypic diversity within an isogenic population, but use a deterministic response when stressors are sensed. Here, we develop a stochastic model of the multiple antibiotic resistance network of Escherichia coli and show that it can produce tunable stochastic pulses in the activator MarA. In particular, we show that a combination of interlinked positive and negative feedback loops plays an important role in setting the dynamics of the stochastic pulses. Negative feedback produces a pulsatile response that is tunable, while positive feedback serves to amplify the effect. Our simulations show that the uninduced native network is in a parameter regime that is of low cost to the cell (taxing resistance mechanisms are expressed infrequently) and also elevated noise strength (phenotypic variability is high). The stochastic pulsing can be tuned by MarA induction such that variability is decreased once stresses are sensed, avoiding the detrimental effects of noise when an optimal MarA concentration is needed. We further show that variability in the expression of MarA can act as a bet hedging mechanism, allowing for survival in time-varying stress environments, however this effect is tunable to allow for a fully induced, deterministic response in the presence of a stressor.

The system is robust to changes in transcription and translation rates, as changes in one rate can be compensated with changes in the opposite direction in the other rate to keep the maximum number of molecules fixed; consequently, the system dynamics are conserved for a wide range of balanced transcription and translation rates.
Dimerization and dimer disruption rates, k dr and k -dr , were assumed to be consistent with the cI dimerization rates in E. coli [8], nevertheless, the system is robust to variation in folding and dimerization parameters; 100fold changes in the values do not have a high impact in noise or protein levels.

Feedback variants and feedback strength
To obtain the variants lacking one or both feedback loops, the parameters were modified as detailed in Table S2. For the Only Negative network, the transcription and MarA degradation rates were changed to allow for a controlled comparison, where the same steady-state is achieved for all salicylate concentration (external equivalence, Fig. S2) and all but two parameters remained unaltered (internal equivalence, Table S2) [9].
In the Table S2, wt is the nominal value of the parameters given in Table S1. Table entries describe how this parameter was modified in each of the four feedback variants. For example, for the rate k a the Wildtype values are used for those variants that keep the positive feedback (i.e. Wildtype and Only Positive), and the parameter is set to zero for the other two variants. In the third row, all the transcription rates are multiplied by the constant specified in the table. The following parameters were set to obtain the same steady state MarA levels for the four variants when uninduced: c 1 = 1/5.416, c 2a = 44.282, c 2b = 1/1.520, c 3 = 3.705. All other model parameters remained the same.

Algorithm
Step 1. Initialize cells by simulating 6-12 hours of growth, with or without antibiotics depending on the specific computational experiment (described below). The final state of the system and the levels of each reactant in the cell are stored to initialize the next simulation round.
Step 2: Simulate all individual cells for a fixed time (5, 75, or 540 minutes depending on the experiment, as described below), given an identical antibiotic time course, using as a start point the final state for that cell in the previous simulation round.
Step 3: Calculate the cost using the MarA levels for the cell. The cost of expressing genes in the mar regulon (e.g. efflux pumps, membrane proteins, stress response genes) is calculated as described in the Methods, using MarA mean levels from every 10th data point to calculate the weighted average of the cost over the entire period.
Cells with total costs above the arbitrary threshold of 0.975 are determined to be dead. However, the results presented here are not sensitive to the exact value of the threshold. Changing the threshold moves the maximum concentration of antibiotic that cells are able to tolerate, but does not affect the overall trends observed.
Step 4: Using the total cost, the number of daughter cells derived from each individual cell is calculated by: The cost calculated from a short simulation (75 minutes if antibiotic is present, or 540 minutes if it is not) is used to extrapolate the final number of cells after 12 hours if antibiotic is present, or 72 hours if it is not. These durations take into consideration the fact that a cell is more likely to grow in an antibiotic-free environment than when antibiotics are present. The extrapolation assumes that the distribution of costs for the population is constant over time.
The total number of daughter cells for each variant is obtained and the proportion of daughter cells for each variant is used to set the levels in the new population. That is, the dead and underperforming cells are replaced by growing cells such that the new ratio between populations is the same as the final ratio in the previous population. This allows us to maintain a constant number of cells, while representing the growth of the total population.
Three classes of growth simulations were conducted. The details for each are described below.
Growth in a fluctuating environment: At each round, the antibiotic was either "on" or "off", with probability of being present equal to 0, 0.25, 0.50, 0.75 or 1 for different runs. If antibiotic was present, the concentration was drawn from an exponential distribution with mean of 3mM and maximum concentration of 4.5mM. Setting a maximum concentration avoids the occurrence of complete population extinction. 200 cells of each variant were simulated, starting with an initialization period of 12 hours of growth with no antibiotic. Next, we simulated 50 rounds of antibiotic stress. If one variant reached 100% of the population prior to 50 rounds, simulations were stopped since this variant would persist for the remaining rounds. For each round, every cell was simulated for 540 minutes in the absence of antibiotic or 75 minutes in the presence of antibiotic.
If antibiotic is present, an initial 5 minute simulation was run in order to determine which cells survived (those with costs below the 0.975 threshold). Here we assumed that if a cell is not able to accumulate enough MarA to grow with a cost below the threshold during the initial 5 minute period it will die before being able to respond to the environment change. Similar results are obtained if the threshold or the time span are varied, as long as the cells are not allowed to reach steady state. Following this 5 minute simulation, a 30 minute simulation was run to allow the surviving cells to reach steady state, then another 75 minute simulation was run at steady state to calculate the average cost for the cell. The average cost was used to calculate the number of replications in 12 hours (without antibiotic) or 72 hours (with antibiotic). Eight independent replicates for each probability of antibiotic addition were run.
Growth in a constant environment: The antibiotic concentration was kept fixed. 200 cells of each variant were simulated, with a 12 hour (if no antibiotic is present) or 6 hours (if antibiotic is present) initialization period. Over the course of 10 rounds every cell was simulated for 540 minutes (in the absence of antibiotic) or 75 minutes (in the presence) and average costs were calculated as indicated for the fluctuating environment. This average cost was used to calculate the number of replications in 12 hours (without antibiotic) or 72 hours (with antibiotic). Three independent replicates were run for each antibiotic level.
Fraction of surviving cells after a pulse of antibiotic: 10,000 cells of each variant were simulated. A pulse of antibiotic was introduced after a 12 hour initialization period with no antibiotic present. The surviving cells were those with cost below 0.975 after the 5 minute simulation. Five independent replicates were run.

Deterministic model
Applying the law of mass action, a description of the system based on the chemical reactions given in Table S1 was developed using ordinary differential equations: For nullcline calculations, the derivatives of all state variables were set to zero. Using Mathematica, we reduced the system into equations to obtain lines showing dA dt = 0 and To evaluate stability of the equilibrium point we calculated the eigenvalues of the Jacobian for the full system of differential equations, evaluated at the equilibrium point.

Explicit Cell Growth Model
Unless otherwise noted, the models treat cell growth and division implicitly, including the dilution rate in the effective protein degradation rate and assuming that cells maintain a constant volume. To test whether this modeling assumption was valid, we conducted simulations that explicitly model cell growth for the Wildtype system with 0mM salicylate (Fig. S7). Model parameters and structure are the same as in the implicit growth model with the following changes: the protein dilution rate as set by cell division was subtracted from the degradation rates given in Table S1 [10], cell volume was modeled as growing exponentially with time [11], and the propensity values were updated as a function of volume (Table S1). We assumed a cell division time of 24 minutes, based on the E. coli doubling time in log phase with rich medium [12]. At each cell division event, the cell volume was divided in two and the molecules were partitioned following a binomial distribution with probability 0.5 [13]. Explicit modeling of cell growth and division has little effect on the results. Importantly, pulses in the MarA are observed in models that treat cell growth and division both explicit and implicitly.