Chaos in a bacterial stress response

SUMMARY Cellular responses to environmental changes are often highly heterogeneous and exhibit seemingly random dynamics. The astonishing insight of chaos theory is that such unpredictable patterns can, in principle, arise without the need for any random processes, i.e., purely deterministically without noise. However, while chaos is well understood in mathematics and physics, its role in cell biology remains unclear because the complexity and noisiness of biological systems make testing difﬁcult. Here, we show that chaos explains the heterogeneous response of Escherichia coli cells to oxidative stress. We developed a theoretical model of the gene expression dynamics and demonstrate that chaotic behavior arises from rapid molecular feed-backs that are coupled with cell growth dynamics and cell-cell interactions. Based on theoretical predictions, we then designed single-cell experiments to show we can shift gene expression from periodic oscillations to chaos on demand. Our work suggests that chaotic gene regulation can be employed by cell populations to generate strong and variable responses to changing environments.


In brief
Chaotic processes are well understood in the physical sciences to generate seemingly random dynamics from purely deterministic mechanisms.The existence and role of chaos in biological systems are less certain.Here, Choudhary et al.  show that chaos explains the variability in the responses of individual bacterial cells to oxidative stress.

INTRODUCTION
The birth of chaos theory was highly significant because it made clear that unpredictable patterns in nature can arise without stochasticity, i.e., purely deterministically. 1 Many seemingly noisy systems were subject to reanalysis and reinterpreted as chaotic rather than stochastic, including examples from biology such as ecological dynamics, [2][3][4][5][6][7] gene expression, 8,9 immune system dynamics, 10,11 neural signal dynamics, [12][13][14][15] circadian rhythms, 16,17 and heart beats. 18However, the underlying causes of dynamics in biological systems are not as well understood as in the physical or chemical sciences. 19,20As a result, the inference of chaos in biology often rests upon mathematical models alone, 3,4,12,[21][22][23][24] which is not sufficient to demonstrate that chaos actually occurs in the biological system itself.
Empirically, chaotic dynamics have been inferred in biological data by statistical detection tools. 2,13,25,26However, it is challenging to distinguish chaotic from stochastic causes in this manner because the inference methods are highly sensitive to measurement noise and random fluctuations that are inherent to all biological processes. 27This problem was recently illustrated by Toker et al., who applied new analysis tools to both physical and biological systems. 28Although they found evidence for chaos in physical and simulated biological data, they concluded that measured heart rate data are stochastic, in spite of a large number of papers having previously concluded that they are chaotic. 18Chaos, then, has proved much more difficult to evidence in biological systems than in physical ones, leaving the importance of chaotic dynamics for biology in doubt.
In response to environmental threats, bacteria have evolved stress responses that drive rapid physiological adaptation. 291][32][33] We serendipitously discovered chaotic behavior in a theoretical model of one of the major bacterial stress responses, specifically the oxidative stress response of the model species E. coli.This discovery allowed us to leverage the detailed understanding and tractability of E. coli to overcome the typical challenges faced when studying chaos in biological systems.Under high hydrogen peroxide (H 2 O 2 ) stress, a cell will strongly induce the expression of proteins that remove H 2 O 2 within the cell and thereby lower the concentration of H 2 O 2 in its vicinity.Our work shows that this response, in combination with the responses of surrounding cells, perturbs the regular periodicity of cell growth, driving the stress response dynamics from periodic oscillations to chaotic fluctuations.By identifying the drivers of chaos, we are able to predict the conditions when it will be present and when it will be lost, and we validate these predictions empirically.In this way, we provide clear experimental evidence of chaos in a living system.Our results further suggest that chaotic gene regulation could be common in the bacterial responses to diverse types of stresses and offer functional benefits.

RESULTS
Chaos is predicted in a bacterial stress response Genetically identical bacteria often display considerable cell-tocell variability in their responses to the environment 34 (Figure 1A).Stress responses can be particularly variable, 35 and it is often assumed that this variability results from stochastic processes inside the cell.7][38][39] Such behavior could arise when the induction of detoxifying enzymes reduces the local concentration of stressor molecules, 40 leading to dynamic feedbacks.Indeed, recent work found that cellular heterogeneity in the response of E. coli to oxidative stress by hydrogen peroxide (H 2 O 2 ) is driven by cell-cell interactions rather than intracellular noise, 41 where the response of each cell is determined by the hydrogen peroxide scavenging activity of its neighbors.This project began with the goal of better understanding the impact of this response on cell-to-cell variability via a numerical model (Figures 1 and S1).The model captures three coupled processes (Figures 1B and S1).First is the dynamics of the intracellular stress response (named S in the model).In E. coli, these occur when H 2 O 2 oxidizes the transcription factor OxyR, which induces expression of H 2 O 2 -scavenging enzymes AhpCF and KatG, and the glutaredoxin-1 GrxA that converts oxidized OxyR back to its reduced form 33,42 (Figure 1C).Second is the inhibitory effect of H 2 O 2 on the growth rate of the cells (named G in the model) (Figure 1D) and, finally, there is the impact of cell-cell interactions (named I in the model) (Figure 1E).These interactions arise because the uptake of H 2 O 2 by one cell can lower the concentration of H 2 O 2 for surrounding cells.
We began by solving a one-dimensional (1D) version of the model, where cells are exposed to a constant external H 2 O 2 concentration from one direction (Figure S1; Video S1).This simple geometry allowed us to explore the behavior of the system at steady state and is highly amenable to empirical testing via microfluidic growth trenches (using a device called the ''mother machine''), which are commonly used in experiments with bacteria. 41,43We then followed the stress response (GrxA expression level) in so-called ''mother cells'' located at the base of the cell group farthest from the source of treatment; the other cells are termed ''barrier cells'' (Figure 2A; Video S1).Given that the model is purely deterministic, with no noise terms, we were surprised to observe seemingly random fluctuations in the stress response of individual mother cells (Figure 2A; Video S1).Moreover, repeated runs of the model yielded highly variable stress response trajectories, which we initially found confusing because model parameters were identical throughout.The only source of variability was that each run of the model began with unsynchronized cells, i.e., at random points in the cell cycles.To explore whether this was the source of the variation, we ran the model with synchronized barrier cells but shifted the initial cell cycle progression of only the mother cell by 2.5$10 À4 % and 5$10 À4 %.The resulting three trajectories were initially indistinguishable but began to diverge significantly after $6 h post treatment with H 2 O 2 (Figures 2B, S2A, and S2B; Video S1 [bottom]).Once diverged, the dynamics became completely different for the three model runs.This extreme sensitivity to small differences in initial conditions of a deterministic model suggests chaotic dynamics. 28,44e next visualized response trajectories as phase diagrams across a range of H 2 O 2 concentrations.These show the defined and closed orbits of periodic oscillations for lower H 2 O 2 concentrations, but at higher H 2 O 2 , we observed the dense aperiodic orbits that are indicative of chaotic fluctuations (Figures 2C and  2D).These fluctuations can also be seen in a bifurcation diagram, which shows the extrema values of the fluctuations for cells as a function of H 2 O 2 concentration.The form of the resulting diagram is characteristic of a chaotic system that shifts from periodic to chaotic regimes as H 2 O 2 concentration increases (Figure 2E).To test formally for chaotic dynamics, we computed the Lyapunov exponent l, which is positive for a chaotic system where a small perturbation in initial conditions leads to exponential divergence of the trajectories.As expected from the bifurcation diagram, cells in the model at higher H 2 O 2 concentrations showed chaos (l > 0), while lower concentrations caused predominantly periodic oscillations in responses (l < 0) (Figure 2F).Note that the frequent and abrupt transitions between periodic and chaotic behavior that we observed in our simulations are typical even for the simplest mathematical models of chaos. 1 An autocorrelation curve of the aperiodic traces also decreased quickly over time, again indicative of chaos, whereas periodic oscillations showed characteristic autocorrelation peaks (Figure 2G).
To explore the generality of our observations, we replaced the detailed model of the oxidative stress response with a simpler version for a generic stress response (Figure 3A).Here, cells take up toxic molecules from their surroundings and protect themselves by producing enzymes that reduce intracellular toxin concentrations (Figures 3B and 3C).Similar to the oxidative stress model, we observed that an increased toxin concentration leads to chaotic enzyme expression dynamics (Figures 3D and  3E).Furthermore, the system becomes more prone to chaotic behavior as the catalytic efficiency of the enzyme or its expression rate increase (Figures 3F-3I).

Observational data also suggest chaos
We next performed experiments on E. coli populations growing in mother machine chips under constant H 2 O 2 treatment, where  S3D) demonstrated an excellent quantitative agreement between our theory and the experiments (Figure 4A; Video S2).
Both showed similar spatio-temporal response dynamics and the same gradient in stress response level along the row of cells (Figure 4A; Video S2).Experiments also matched the theoretical prediction that sudden H 2 O 2 treatment triggers an induction of stress response, which coincides with a transient dip in the cell elongation rate followed by adaptation (Figures 4B, 4C, S3E, and S3F).Importantly, like the model, the mother cells in experiments displayed large fluctuations in stress response level during constant H 2 O 2 treatment (Figure 4D).The autocorrelation curves from these traces decreased quickly over time, indicating a lack of periodicity in the dynamics, consistent with chaotic behavior (Figure 4E).We further applied the ''chaos decision tree algorithm'' of Toker et al., 28 which is an analysis pipeline that uses the permutation entropy to categorize dynamics as stochastic or deterministic.The pipeline classified 92% of the measured stress response trajectories as deterministic (Figures 5A, 5B, and S3M-S3Q), supporting the prediction of the model that the fluctuations are predominantly a consequence of deterministic chaos and not caused by noise.Most of the 8% of trajectories that were categorized as stochastic corresponded to cells that had died at the onset of H 2 O 2 treatment (Figures 5A, 5B, and S3M-S3Q).Active growth dynamics are hence required for chaotic response behavior.
As another test to distinguish deterministic chaos from noise, we applied the Grassberger-Procaccia algorithm. 45This algorithm identifies the correlation dimension (or fractal dimension), which should be low when the response fluctuations are driven by a deterministic process with a small number of effective variables, but tends to infinity for a truly stochastic process. 46We found a finite correlation dimension of $2 for the GrxA dynamics in experiments and simulations, indicating determinism in the system (Figure 5C).The analysis, therefore, is consistent with our model prediction that deterministic response fluctuations are generated by a simple cyclic cell growth pattern that generates oscillations in the number of barrier cells.Indeed, GrxA dynamics were negatively correlated with changes in the number of barrier cells in experiments, consistent with the expectation that an increase in the number of barrier cells reduces the local H 2 O 2 concentration and thus the stress response level (Figure S4A).
The use of the mother machine allowed us to follow cell trajectories over long periods, which is important for our ability to test the model's predictions.However, the 1D structure introduced by the mother machine is also potentially unrepresentative of the way that bacteria normally grow and respond to stresses.To study cells in a more realistic setting, we imaged two-dimensional (2D) microcolonies with continuous H 2 O 2 treatment (Figures 5D, S5A, and S5B; Video S3).Although these experiments cannot monitor cells over a long time, we again observed substantial heterogeneity in responses between individual cells, both within and across microcolonies, which is consistent with chaotic behavior (Figures 5E, S5C, and S5D).Moreover, as the colony grows and expands, the stress response eventually decreases and becomes more uniform (Figure S5E).This behavior is predicted by the model: scavenging reduces the H 2 O 2 concentration in a larger colony and, with this, the potential for chaotic responses is predicted to decrease (Figure S5F).

Experimental tests make or break chaos
Our model predicts the existence of chaotic behavior in a biological system-the oxidative stress response of E. coli-and our observational data are consistent with chaos.Together, these two approaches lend support for chaotic behavior and they reflect the typical standard of evidence in biological systems, where modeling predicts chaos and/or observation of seemingly chaotic dynamics are reported.However, there are problems with such evidence.Most obviously, a modeling prediction is just that; it does not demonstrate that chaos actually occurs in a biological system.Second, tests for chaos from observational data are challenging when the underlying causes of the dynamics are uncertain.
We, therefore, sought to leverage the tractability of our study system to provide strong evidence of chaos in a biological system via targeted perturbation.In particular, if the dynamics are indeed deterministic, then it should be possible to shift the responses away from the chaotic regime into the parameter space where periodic oscillations occur, whereas this should not be possible for stochastic fluctuations (e.g., caused by gene expression noise).
To evaluate this prediction, we returned to our model to identify changes that remove the chaos from the dynamics and shift them to periodic oscillations.This analysis revealed that feedback between each of the three components of the model-cell growth (G), cell interactions (I), and stress response (S) (Figure 1B)-is required for the emergence of chaotic dynamics (Figures S2C and S2F).Specifically, S is always required as it is the core of the stress response.Without the cell-cell interaction component of the model, no fluctuations were observed for the GrxA traces (Figure S2C).When the cell growth component was uncoupled from the model, i.e., the growth rate was unaffected by the H 2 O 2 treatment, then the response oscillations were no longer chaotic but periodic (Figure S2F).
From here, we identified parameter changes for each component of the model that are predicted to shift stress response dynamics from chaotic to periodic.These changes were as follows: (1) reduce cell growth rate to slow down the oscillations in the number of cells per trench (G) (Figures 6A and S2F-S2H), ( 2) reduce cell numbers to lower the effects of cell-cell interaction (I) (Figures 6B and S2C-S2E), and (3) reduce the stress response (S) (Figures 2 and 6C).For each case, we then devised a way to make this manipulation experimentally (Figures 6A-6C In each case, we followed the dynamics of the stress response in mother cells as before and used autocorrelation analysis to test for chaos.If the dynamics are periodic, one will see a characteristic autocorrelation that peaks at the frequency of the periodicity in the data.If the dynamics are chaotic, by contrast, no such peak in the autocorrelation is seen.However, if the fluctuations are stochastic, then we should not observe autocorrelation peaks under any condition.As expected, applying this test revealed no peak in the autocorrelation function for experimental conditions that yield chaos (cyan traces in Figures 6A-6C, S6A, S6C, and S6E).By contrast, in all three cases designed to remove chaos, we observe peaks in the autocorrelation function (black traces in Figures 6A-6C, S6B, S6D, and S6F).
In summary, we were able to identify three conditions that break the chaos in the model.We then demonstrate that making these manipulations in experiments also shifts cell responses away from chaotic to periodic oscillations.

A simple deterministic process explains chaotic oxidative stress response fluctuations
The above analyses all strongly support a deterministic origin to the response dynamics.However, this does not imply the responses are entirely devoid of any noise.The number of OxyR molecules and H 2 O 2 scavenging enzymes per cell is expected to fluctuate randomly due to gene expression noise. 47Might these fluctuations still be important for the response dynamics when combined with the deterministic causes?To investigate this, we added gene expression noise to our model by introducing a stochastic term in each molecular component of the stress response model (S*) and coupled it with the cell growth or cell interaction model components as before (Figure S4).Although the resulting response dynamics superficially matched aspects of the experimental data, they disagreed with key features.Specifically, none of the noisy response models (S*, S*+I, S*+G) showed both (1) loss of autocorrelation peaks for higher H 2 O 2 concentrations and (2) negative cross-correlation between response fluctuations and the number of cells per trench.These findings again support our conclusion that the (legend continued on next page) fluctuations are predominantly deterministic.Only when we coupled the noisy response with the full model of growth and cell interactions (S*+G+I) did we recover the dynamics seen in experiments, but this was already the case for the model without noise.Therefore, the model can tolerate the addition of noise but does not require it for the generation of unpredictable response fluctuations.
We find further evidence of the importance of deterministic processes over stochastic ones when we apply the Grassberger-Procaccia algorithm to the noisy response models.In all cases, the noisy models do not converge on a low correlation dimension, whereas the deterministic model had a correlation dimension of $2 (Figures S4G-S4I).This means that although the model includes multiple dynamic variables for each of the interacting cells in a trench, the response dynamics of one focal cell are actually driven by a deterministic process with only $2 effective variables.Strikingly, the experimental data had the same effective correlation dimension of $2 across all growth conditions used in our tests for chaos (varying growth rate, trench length, H 2 O 2 concentration), irrespective of whether the dynamics were periodic or chaotic (Figures 6D and S4G-S4I).This finding suggests that one simple deterministic process is responsible for generating all the observed response dynamics in experiments, from periodic to chaotic.Moreover, both our model and experiments suggest the cause of this determinism is strongly linked to the cell cycle (Figures S2F-S2H and S4A).Indeed, the timing of the peaks in the autocorrelation curves matches the cell cycle duration in each experimental condition, which is again predicted by the model (Figure 6E).
How does one process generate periodic oscillations under some conditions and chaotic fluctuations under others?Periodic oscillations occur at low stress levels when the cell growth rate is constant.Here, the number of cells in a given locality oscillates at regular intervals as cells grow and divide, leading to periodic changes in the H 2 O 2 concentration around each cell.At higher stress treatments, the growth rate of each cell becomes sensitive to the H 2 O 2 concentration.The cells then grow at variable rates and their local numbers change without a fixed period, leading to irregular oscillations in H 2 O 2 .In this way, the stress response dynamics transition to chaos.

DISCUSSION
Our work has identified chaos in a bacterial stress response.This finding shows that seemingly random phenotypic heterogeneity can be generated by deterministic rather than stochastic processes.That is, regulatory circuits are able to generate unpredictable outputs, even when the underlying mechanisms are entirely deterministic.Such cases are called chaotic because they have the property of amplifying infinitesimally small differences in the initial conditions to an extent that forecasting the long-term behavior is impossible-no matter how accurately the initial conditions can be defined.Although noise due to molecular fluctuations is certainly present in cells, our model and experiments show that phenotypic heterogeneity among cells can arise without the requirement for stochasticity.Moreover, our analyses suggest that the key driver of both the periodic oscillations and chaotic fluctuations is the determinism of the cell cycle (Figures 6, S2F-S2H, S6E, and S6F).
We have focused here on the behavior of cells as they grow in lines in channels of the mother machine microfluidic device.The great advantage of this study system is that one can follow stress response dynamics of individual cells for relatively long periods in a bacterial population where cell-cell interactions are still preserved.These relatively long time series from single cells were important for our ability to both identify and describe chaos empirically.However, our work suggests that the chaotic behavior seen in this system also occurs in more complex and realistic growth conditions.When we study cells in 2D colonies, we also see unpredictable stress response dynamics, which are consistent with chaos (Figures 5D and 5E).It will be interesting to understand how chaotic processes are affected in even more complex growth environments, such as submerged biofilms.On the one hand, additional variation in the local density and arrangement of cells should increase the potential for chaos, while on the other the increasing capacity to scavenge H 2 O 2 within larger communities may reduce stress levels and the duration of any chaotic dynamics.
Why have bacterial cells evolved to display such chaotic behavior?Our work shows that multiple processes interact to generate chaos, including the cell cycle and changes in local cell density.Another important factor is the strength of the stress response itself: a high expression rate and high catalytic efficiency of scavenging enzymes are critical to the transition from periodic oscillations to chaos (Figures 3 and 6).The evolution of chaotic behavior in the oxidative stress response may, therefore, lie in the benefits of a strong response for surviving stress, which is likely to provide a strong selective advantage to cells. 33,48This importance of a strong response for chaos mirrors the classic theoretical results on chaos in population biology.There, models predict that high population growth rate is needed to generate chaotic dynamics 44 because this causes the population to overshoot equilibria and over-compensate. 49haos also results in a variable response among cells.An intriguing possibility, therefore, is that chaotic responses could have the benefit of diversifying cell behavior as a population bet-hedging strategy against unpredictable stresses. 50Although microbes can harness intracellular molecular noise to generate phenotypic heterogeneity, deterministic chaos can achieve this these peaks are absent for chaotic dynamics (teal) in the case of (A) growth rate perturbation ( without noise disturbing the response accuracy of each individual cell.In fact, our study showed that the conditions that lead to chaos are exactly those where bet-hedging is valuable in principle, namely at high stress levels and in cell populations but not in isolated cells. Although our work is based upon the detailed study of one bacterial stress response, the existence of chaotic behavior may be widespread in cellular systems for related reasons. 3,51sing a generalized model, we find that chaotic cell responses are possible whenever the absorbance of a stressor reduces a cell's growth rate and lowers the stressor concentration of the surrounding cells (Figure 3).The feedback between these effects creates spatio-temporal dynamics that amplify small perturbations, leading to chaotic behavior.Hence, the survival strategies of bacteria and other cells exposed to stressors-such as antibiotics, antimicrobial peptides, or reactive chemicals-all have the potential for deterministic chaos.
The study of chaos in biology has received considerable attention, and there are many potential examples where unpredictable dynamics have been observed that appear chaotic. 28However, it remains challenging to identify chaos from either observational data or theoretical models alone.Here, we have presented a different strategy, which rests on the ability to manipulate chaos.In addition to the prediction of chaos and observational experiments, our model also correctly predicts the conditions where chaotic dynamics are lost.This close fit between model and measurements provides clear evidence for the existence of chaos in the single-cell dynamics of bacterial stress responses.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
We performed experiments with bacterial strains that were derived from E. coli K-12 AB1157.The description of genetic modification and growth conditions is described in sections below.REAGENT

Strains and plasmids
All experiments were performed with a strain derived from E. coli K12 AB1157 that was previously described in Choudhary et al. 41 The strain constitutively expressed P RNAI -mKate2 fluorescent marker for cell segmentation analysis and the flhD gene was deleted to inhibit flagellar motility allowing growth in mothermachine microfluidic chips.The OxyR response reporter plasmid carrying PgrxA-SCFP3 was derived from an E. coli promoter library of pSC101 plasmids. 54Each plasmid in the library contains the promoter region of a specific gene or operon in front of GFPmut2 fluorescent protein.We changed the GFPmut2 to the fast-maturing cyan fluorescent protein SCFP3 55 using Gibson Assembly (NEB).The promoter region was confirmed by sequencing and the plasmid was transformed yielding strain SU777 (AB1157, DflhD, P RNAI -mKate2, mutL-mYPet, carrying plasmid pUA139 PgrxA-SCFP3 with kanamycin resistance gene).Presence of the expected fluorescent protein signal was verified by taking microscopy snapshots.

Media and growth conditions
Cells were grown at 37 C for all experiments.Cells were streaked from glycerol stocks stored at -80 C on LB agarose plates with 25 mg/mL kanamycin.A single colony was picked and grown overnight shaking in 4 mL M9 minimal media.This media was prepared with M9 salts (15 g/L KH 2 PO 4 , 64 g/L Na 2 HPO 4 , 2.5 g/L NaCl, and 5.0 g/L NH 4 Cl), 2 mM MgSO 4 , 0.1 mM CaCl 2 , 0.5 mg/mL thiamine, MEM amino acids, 0.1 mg/mL L-proline, and 0.2% carbon source (glucose or glycerol).The next day, overnight culture was diluted 1:50 and grown shaking to OD600 $0.3 in 4 mL M9 minimal media.For loading cells in microfluidic chips, 0.85 mg/mL Pluronic F127 was added to the media to avoid cell aggregation.For experiments done under hydrogen peroxide treatment, the specific concentration of H 2 O 2 was added to the growth media immediately before the start of the experiment.LB (10% V/V) was added to M9 glucose media for certain experiments, as specified in the figures.

Microfluidics experiments
Mother machine chip preparation Single-cell imaging was performed using the 'mother machine' microfluidic device as described in 43,56 .The chip has a main channel for flow of media, branching into perpendicular growth channels (here called 'growth trenches') of dimension 1.2 mm width and 1.2 mm height and 25 mm length.The chips were made of polydimethysiloxane (PDMS, Dow Corning Sylgard 184 kit) polymer using a silicon wafer mold (Conscience).A 1:10 solution of polymerising agent and PDMS monomer were rigorously mixed and then poured onto the silicon wafer.This was placed in a vacuum chamber and pressurised to remove air bubbles.The device was then heated at 65 C in an oven for 2 hours to polymerise.For each experiment, one chip was cut out using a scalpel, and holes for inlet and outlet were inserted using a 0.75 mm biopsy puncher.The device was cleaned using 100% ethanol and dried with nitrogen gas.The cleaning was repeated 3 times.The PDMS chip was bonded on a glass coverslip (thickness No 1.5).These coverslips were first cleaned by sonication with acetone for 20 mins followed by isopropanol for 20 min, and then dried with nitrogen gas.The cleaned coverslip and PDMS chip were exposed to air plasma for 2 min and bonded at 95 C for 30 min.
Where indicated, a different silicon wafer was used to generate mother machine chips with shorter trenches of 10 mm length and 1.2 mm width and 1.2 mm height.Here a ''negative'' mould was prepared first using PDMS as intermediate from a silicon wafer (by mixing monomer and curing agent 1:5) that was then used to prepare the ''positive'' chips using the method explained above.Mother machine setup 1 mL of exponentially growing cells were spun down for 2 min at 6000 rpm.Cells in the pellet were resuspended in 100 mL of the supernatant and loaded in the microfluidic chips by pipetting through the inlet.The chip was then inserted into a custom-built centrifuge holder and spun at 5000 rpm for 10 min to aid the loading of cells into the growth trenches.50 mL syringes were filled with M9 minimal media containing Pluronic F127 and H 2 O 2 as indicated.The syringes were attached to silicon tubing (Tygon) and loaded onto syringe pumps (NewEra SyringePumpPro) to deliver media into chips at a constant flow rate of 2.5 mL per hour.Cells were initially grown without H 2 O 2 for $3 hours before switching the inlet media to a syringe containing H 2 O 2 .Microcolony chip preparation E. coli microcolonies were grown on 1% agarose pads made with M9 glucose + 10% LB media.The procedure of preparing these pads is shown in Figure S5A.Melted agarose solution was poured into the top of an empty 50 mL syringe plunger head wrapped with adhesive tape to act like a container.A glass cover slip was then placed on top of the taped cylinder.After the agarose was set, the cover slip was removed and 1 ml spots of overnight culture were dropped on the flat agarose surface and grown for 2 hours at 37 C.After 2 hours, the agarose was removed from the plunger head, thus leaving a conical dip in the agarose that acted as reservoir for adding H 2 O 2 treatment solution.The agarose pad was inverted to sandwich the cells between the agarose and a cover slip for imaging.A needle was inserted through the adhesive tape to continuously flow in growth medium with H 2 O 2 using syringe pumps.The medium dripped onto the conical dip and diffused through the agarose to reach the cells below.

Microscopy
Time-lapse microscopy of cells in mother machine chips Time-lapse imaging was performed using a Nikon Ti-E inverted fluorescence microscope equipped with 100x NA 1.40 immersion oil objective, motorized stage, sCMOS camera (Hamamatsu Flash 4), LED excitation source (Lumencor SpectraX), and operated with a perfect focus system.Exposure times were 100 ms for P RNAI -mKate2 (l = 555 nm) and 75 ms for sCFP3 reporter (l = 440 nm) using 50% of maximal LED excitation intensities.The excitation and emission lights were separated using a triband dichroic and individual emission filters.The microscope chamber (Okolabs) was maintained at 37 C throughout the experiments.Images were captured every 3 min for the 2 emission channels.Time-lapse microscopy of cells in microcolonies Time lapse imaging was performed on a Nikon Ti-E microscope equipped with a 100x NA 1.45 oil immersion objective, motorised stage, sCMOS camera (Photometrics Prime95B), LED excitation source (Lumencor SpectraX) and perfect focus system.Exposure times were 100 ms for P RNAI -mKate2 (l = 555 nm) and 75 ms for sCFP3 reporter (l = 440 nm) using 50% of maximal LED excitation intensities.The microscope chamber (Okolabs) was maintained at 37 C throughout the experiments.Images were captured every 3 minutes with phase contrast and the two fluorescence channels.

QUANTIFICATION AND STATISTICAL ANALYSIS
Mother machine data processing and analysis Time-lapse microscopy data were saved as .nd2files and visualized in Fiji. 52The data were processed using the BACMMAN plugin in Fiji as described in 53 and further analysed using custom Python and MATLAB scripts.Images were first pre-processed by BACMMAN using the P RNAI -mKate2 fluorescence channel to stack all individual growth trenches and correct for experimental drift in x-y coordinates and image rotation.The outlines of cells in the growth trenches were then jointly segmented and tracked over time based on the P RNAI -mKate2 fluorescence signal.The traces were visually inspected and manually corrected for errors in segmentation or lineage tracing using the BACMMAN software.The SCFP3 fluorescence was extracted by overlaying the cell masks from the P RNAI -mKate2 channel onto the SCFP3 channel and computing the mean intensity over the cell area.BACMMAN generated output in 3 excel files containing cell growth characteristic, P RNAI -mKate2 intensity data and SCFP3 intensity data.These files were then further analyzed using Python and MATLAB code as described below.

Microcolony image analysis
Segmentation of cells growing in microcolonies was performed based on the P RNAI -mKate2 fluorescence signal and using the MicrobeTracker tool in MATLAB 57 followed by manual correction of the segmentation masks.These outlines were then applied to the CFP channel and a MATLAB script was used to quantify the average intensity per cell area.Cell lineage tracing was performed manually and custom python code was used to plot CFP intensity traces.

Chaos decision tree algorithm
To categorise if the fluctuations in PgrxA-SCFP3 traces of individual mother cells are deterministic or stochastic, we applied the 'Chaos decision tree algorithm' as described by Toker et al. 28 The pipeline is available as MATLAB code.Briefly, the algorithm tests for stochasticity by computing the permutation entropy using the cyclic phase permutation algorithm.The permutation entropy quantifies the extent to which the values in a trace are ordered or random in time.The value of the permutation entropy for the original trace is compared to many randomly shuffled versions of the same trace in which the temporal order of the data points is removed while the mean and standard deviation are maintained.If the fluctuations are stochastic, then the permutation entropy is similar for the original and shuffled traces; otherwise the fluctuations are classified as deterministic.The chaos decision tree further tests for stationarity.i.e. whether statistical properties like the mean and standard deviation of a trace do not change over time; else the trace is classified as non-stationary.

Deterministic vs stochastic and live vs dead categorization of experimental traces
The MATLAB code described in Toker et al. 28 as explained in the section above was used to categorise PgrxA-SCFP3 traces of individual mother cells as showing deterministic or stochastic fluctuations.Traces were analysed at steady-state from 1-hour after the start of H 2 O 2 treatment.The traces were pre-processed with a moving-mean filter using a window of 3 time points (i.e. 9 min) twice.The rationale for this is that we are principally interested in the large-scale response fluctuations on the time-scale of the cell cycle ($50 to 100 min, depending on conditions).Only traces with at least 5.5 hours of data were included in this analysis (length of data required for entropy calculation in the pipeline).
PgrxA-SCFP3 traces were further categorised as originating from live or dead mother cells.A cell was considered dead if the length growth rate as computed by BACMMAN was below 0.0024 min -1 .

Autocorrelation analysis
Autocorrelation analysis was performed to distinguish between periodic and chaotic fluctuations of the PgrxA-SCFP3 traces of individual mother cells.Traces with at least 2-hours of data were analysed at steady-state from 1-hour after the start of H 2 O 2 treatment.The traces were pre-processed with a moving-mean filter using a window of 3 frames (i.e. 9 min) twice.Due to the slower growth of cells in M9 glycerol media, a window of 6 frames was used and traces were analysed at steady-state from 2-hours after the start of H 2 O 2 treatment.ACF curves were computed from the DPgrxA difference signal, which was calculated by subtracting the PgrxA-SCFP3 values of consecutive frames.The Python stattools.acffunction from the statsmodel library was used to output the autocorrelation value over a range of lag times.Mean ACF curves were computed by averaging the ACF values from the single-cell traces at each lag time.A peak finding algorithm was applied to the mean ACF curves to quantify the period of non-chaotic PgrxA-SCFP3 traces.This was done in python using the find_peaks function in the scipy.signallibrary with a prominence of 0.2.ACF curves for simulated GrxA traces were computed during steady-state from 700 to 5200 min after start of H 2 O 2 treatment at t = 50 min.The analysis was done as for experimental data but without moving-mean filtering.
Correlation Dimension to discriminate deterministic and stochastic processes We used the Grassberger -Procaccia algorithm 45 to estimate the correlation dimension (D corr ) of the GrxA response dynamics, which can be understood as the effective number of dynamic variables that generate the response fluctuations.A large number implies a stochastic process whereas a low number results from a deterministic process. 46o compute the correlation dimension, we followed the procedure described in Sandler et al. 46 For a given trace of {GrxA n } where n is any time point, the algorithm constructs a vector {GrxA n , GrxA n+1 , GrxA n+2 ,., GrxA n+E-2 , GrxA n+E-1 } considering them in E-dimensional space (E is also called embedding dimension).For a data set with N data points in an embedding dimension E, the correlation sum CðrÞ is quantified to then compute d corr : The Heaviside step function ƟðxÞ (ƟðxÞ = 1 if x > 0 and ƟðxÞ = 0 if x % 0) is used to compute the fraction of data points x i and x j that are within a distance r of each other.The values of r were evenly spaced in log scale from 0.01s to 3s, where s is the standard deviation of the fluctuation of the given GrxA trace.
The correlation dimension CðrÞ follows a power law for small r, such that CðrÞar Dcorr .Therefore, D corr can be estimated from the slope of a log-log plot of CðrÞ versus r.CðrÞ grows monotonically with r.
The D corr obtained for different values of E is plotted against E. For a deterministic process, the D corr vs E plot saturates and the saturating value of D corr corresponds to the effective dimension of the dynamic process, i.e. the number of dynamic variables that determine the response dynamics.In contrast, a truly stochastic process yields a straight line for D corr vs E, with E =D corr for any value of dimension E, implying that the process is determined by an infinite number of dynamic variables.
We applied a moving-average filter with a filtering window of 3 time points (i.e. 9 minutes) to smooth experimental PgrxA-SCFP3 traces before applying the Grassberger -Procaccia algorithm.The rationale for this is that we are principally interested in the largescale response fluctuations on the time-scale of the cell cycle ($50 to 100 min, depending on conditions).Nevertheless, applying the Grassberger -Procaccia algorithm to experimental data without smoothing also produced a low correlation dimension (Figures S4H  and S4I), consistent with the dynamics being overall deterministic.

Lyapunov exponent computation
Lyapunov exponent was computed for simulated GrxA traces of mother cells (positioned at the closed end of growth trenches) during steady-state from 700 to 5200 min after start of H 2 O 2 treatment.The traces were normalized in MATLAB using the ''normalize'' function that rescales the data with mean of 0 and standard deviation of 1.The Lyapunov exponent was computed in MATLAB using the ''phaseSpaceReconstruction'' function and ''lyapunovExponent'' function with fs = 10.

Cell growth model
The growth of rod-shaped E. coli cells was described by an adder model with equal length added to the cell over time between successive division events. 58The elongation rate g determines the increase in cell length lfrom time t to t + 1: The cell divides into 2 daughter cells of equal length when the total cell length added from the birth length exceeds 2 mm.We modelled the inhibition of the cell elongation rate by H 2 O 2 as a sigmoidal decay: where g 0 is the elongation rate without treatment, ½H 2 O 2 cell is the intracellular H 2 O 2 concentration, and the sigmoidal decay factors c 1 ; c 2 determine the sensitivity of the elongation rate to ½H 2 O 2 cell .

Cell-cell interaction model
Scavenging of H 2 O 2 results in concentration gradients that are determined by the spatial arrangement of cells and the diffusion of H 2 O 2 with diffusion coefficient D. Hence, ½H 2 O 2 external is a function of a cell's position x along the length of the growth trench.The spatial arrangement was modelled as a one-dimensional colony of rod-shaped cells with radius R C in cuboid growth trenches with square cross section W 2 .This is analogous to the growth trenches in a mother machine device 41 with H 2 O 2 treatment entering through the open end of the trench (where x = 0).The scavenging rate of H 2 O 2 is limited by the absorption of H 2 O 2 across the cell envelope 41 with absorption rate constant k abs .We modelled the spatial profile of ½H 2 O 2 external based on the reaction-diffusion equation, using a similar approach as described in Yang et al. 59 : We first consider the profile of ½H 2 O 2 external along the length of a single cylinder-shaped cell of length l t (the hemispherical caps at the cell poles will be considered below).At each position x, ðW 2 --pR 2 c Þcorresponds to the area around the cell where H 2 O 2 diffuses through the trench.The circular perimeter of the cell where H 2 O 2 is absorbed is 2pR c .At steady-state, v½H2O2 external vt = 0. Therefore, the reaction-diffusion equation becomes: 4) To solve this equation for ½H 2 O 2 external , we introduce the coefficient l: 5) We assume the non-adsorbing boundary condition at x =l t .The solution to this equation is: 6) Here, ½H 2 O 2 external;0 is the concentration at x = 0. We modify Equation 7to account for the shape of E. coli cells as cylinders capped with hemispheres.First, we multiply Equation 5 by the cell length l t : 7) Compared to the cylinder geometry, the free volume in the trench through which H 2 O 2 diffuses increases by 2pR 3 c --4 3 pR 3 c due to the two hemispherical caps.The surface area for H 2 O 2 absorption is unchanged.Hence: 8) The solution to the equation is½H 2 O 2 external = ½H 2 O 2 external;0 cosh x À l t l 0 cosh l t l 0 (Equation 9) Where, l 0 = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi [H 2 O 2 ] cell oxidises the transcription factor OxyR from its reduced to oxidised form, where K ox is the 2 nd -order oxidation rate constant.The total OxyR concentration was assumed constant, such that ½OxyR Ox = ½OxyR total À ½OxyR Red .Oxidised OxyR is reduced by GrxA with Michaelis-Menten kinetics where K red is the catalytic rate constant and h OxyR is the Michaelis constant.Hence: 10) The OxyR regulon includes numerous genes involved in various aspects of oxidative stress tolerance.We reduced this system to the two key H 2 O 2 scavenging enzyme genes katG and ahpC, and the glutaredoxin-1 grxA.Each gene has a basal expression rate (R grxA;basal , R katG;basal , R ahpC;basal ) and an inducible expression rate that depends on the concentration of oxidised OxyR (½OxyR total À ½OxyR Red ).The maximal induced-expression rates are given by K grxA;act , K katG;act ; K ahpC;act .The parameters h grxA;act , h katG;act ,h ahpC;act define the concentrations of ½OxyR Ox that give half-maximal induction of each gene.Gene expression is counteracted by dilution due to cell growth with rate g, given by the growth model (Equation 2).The gene expression dynamics follow:  14)

General stress response model
We formulated a general model to study the emergence of chaos in the response to a generic toxin.The purpose of this model was to abstract as many of the molecular details as possible and explore if a simple gene regulatory circuit is still capable of producing chaotic dynamics when coupled with the growth model and the cell-cell interaction model described above.The intracellular toxin concentration ½Toxin cell is determined by the influx of external toxin with rate constant R influx , and detoxification by an enzyme E where K cat is the catalytic rate constant and h E is the Michaelis constant: 15) The detoxifying enzyme is produced at an inducible expression rate that depends on the intracellular toxin concentration.The maximal induced-expression rate is K act and h act defines the toxin concentration that gives half-maximal gene induction.Enzyme expression is counteracted by dilution due to cell growth with rate g, given by the growth model (Equation 2).16) As above, the cell elongation rate depends on the intracellular toxin concentration: 17)

Model simulations
Software Simulations of the model were performed using custom-written Python code.The following libraries were used: pandas, numpy, math, scipy, random, and matplotlib.

Growth model parametrisation
Under the experimental conditions of our study, the mean elongation rates of E. coli cells growing without H 2 O 2 in M9 glycerol, M9 glucose and M9 glucose + 10% LB media are g 0 = 0:025min -1 ,0:042min -1 and 0:065min h AhpC;act = 0:1 mM; h KatG;act = 0:18 mM; h GrxA;act = 0:1 mM; K AhpC;act = 0:2 mM min À 1 ; We estimated the calibration factor k GrxA;calibration which allowed to compare the PgrxA-SCFP3 intensity from experiments to the GrxA concentrations in the model simulations, such that: We varied k GrxA;calibration from 0 to 2500 with steps of 50 and minimized the MAE between GrxA model;calibrated and GrxA experimental .For a given value of k GrxA;calibration , the MAE was estimated from the PgrxA-SCFP3 intensity values of the cells at the open end of the growth trenches in M9 glucose media at steady-state treated with various concentrations of H 2 O 2 (12.5 mM, 25 mM, 37.5 mM, 50 mM, 62.5 mM, 75 mM, 82.5 mM, 100 mM and 500 mM).

MAE =
GrxA experimental;i À GrxA model;calibrated;i GrxA experimental;i Where i denotes the different H 2 O 2 concentrations.The MAE plotted against k GrxA;calibration showed a minimum at 1150 [Figures S3A and S3B].General stress response model parametrisation For the plots in Figure 3, we used the following parameter values to solve the general stress response model: h E = 1a:u:; K act = 0:1a:u:min À 1 ; h act = 1 a:u:; K cat was varied between 0 to 20min À 1 ; R influx = 1 min À 1 , ½Toxin external was varied between 0 to 200 a.u.The parameters used in the growth model and cell-cell interaction model were the same as for the H 2 O 2 stress response model.

Simulation input
Simulation runs were initialised by specifying the user-defined parameters: number of time points T, time step duration Dt (1 min), [H 2 O 2 ] external , time of H 2 O 2 treatment (time point 50, unless mentioned otherwise), number of growth trenches (n g ), length of growth trenches (L trench ), initial cell lengths at t 0 .Unless otherwise specified, the initial cell lengths were drawn from a random distribution to match experimental conditions where cells are loaded into growth trenches from an unsynchronised culture.Cells are positioned in a straight row with their poles touching.The values used for simulations are given in the Methods S1.Simulation procedure Simulating the interlinked cell-cell interaction model (I), stress response model (S), and growth model (G) required discretisation in space and time.The simulation procedure is shown schematically in Figure S1.
First, the cell-cell interaction model predicts the external [H 2 O 2 ] external concentration that each cell is exposed to according to its position in the growth trench.The reaction-diffusion equation (Equation 9) is solved for the uptake of  2) from the growth model as inputs for each cell.
Thirdly, the growth model uses the intracellular [H 2 O 2 ] cell concentration from the stress response model as input to compute the cell elongation rate g for each cell (Equation 2), and changes the number, size, and positions of all the cells for the next time point (Equation 1).Cell elongation and division pushes cells towards the open end of the trench.If the summed cell lengths exceed the trench length then the outermost cell is removed at the next time point.
The outputs of the interdependent models at one time point are used as input conditions for the next time point.The simulation runs for a set number of time points.Parallel growth trenches are treated as independent simulation runs.

Simulation visualization
We generated videos to visualise the simulation results using custom Python code and the following libraries: Image, ImageDraw from PIL, tifffile.We drew rows of rod-shaped cells according to the outputs of the growth model and the intensity of each cell was given by a linear grayscale conversion of the GrxA value from the stress response model.Simulations with multiple growth trenches were combined into the same image and converted into .tifffiles.The tiff files were concatenated over time into videos using Fiji. 52Further data visualisation and extraction (e.g.lineage tracing) was performed on the simulated videos using BACMMAN software 53 as for experimental data, described above.

Model of the oxidative stress response with noise Stochastic differential equations
The stress response model was modified using the Langevin approach [63][64][65] to account for stochasticity in gene regulation.This approach involves adding a stochastic noise term to the set of differential equations ( Equations 10, 11, 12, 13, and 14).Specifically, for the deterministic differential equation of the form: = fðXÞ (Equation 18) where X = ½OxyR Red ; ½GrxA; ½KatG; ½AhpC, we add a noise term h: Where W is the variable of a Weiner process satisfying the condition that W t À W 0 $ Ɲ. 66 Ɲ is a normal distribution with mean of m and variance of s 2 : Rearranging Equation 17gives the stochastic differential equation (SDE) as follows: dX = fðXÞdt + sdW (Equation 20)

Euler-Maruyama method
To solve Equation 20 numerically, we use the Euler-Maruyama method. 67,68Let Dt = t i+1 À t i where i represents the iteration step from 0,1, .to N (number of iterations).Therefore, Equation 20 can be rewritten as: 21) Where Z i represents the normal distribution (Weiner process derivatives) with mean = 0 and variance = 1.

Noisy equations for oxidative stress response
To solve the equations numerically, the equations were non-dimensionalised in time t nd = t$g, where g is the elongation rate.Then, noise term ɣ gene was added to each equation, giving the following:  29) The noise term is added in the form ɣ gene = s gene Z i ffiffiffiffiffiffiffiffiffi ffi Dt nd p For our simulations Dt nd = 5$10 À 6 min À 1 and {s OxyR ; s GrxA ; s KatG ; s AhpC } = {0:0002; 0:02; 0:02; 0:02}.We chose values of s such that the magnitude of the fluctuations of the stochastic response model (S*) was similar to the magnitude of the chaotic fluctuations of the full deterministic model (S+G+I).MAE was computed between the CV (Coefficient of variation) of GrxA from the deterministic model and the stochastic stress response model solved for a range of concentrations of H 2 O 2 for s of [0.00005, 0.0001,0.0005,0.001,0.005,0.01,0.05,0.1].

MAE = CV deterministic À CV stochastic CV deterministic
This model was then coupled with the deterministic G and I models individually (S*+G, S*+I) or with both G and I models (S*+G+I), as required.

Figure 1 .
Figure 1.Modeling the oxidative stress response in bacterial populations (A) Environmental stress, such as H 2 O 2 exposure, induces heterogeneous responses in bacterial populations, which could be caused by stochastic or deterministic mechanisms.(B) H 2 O 2 affects cell growth rates (G) and triggers an intracellular stress response (S) that creates stressor gradients by cell-cell interactions (I).S-G-I feedback modulates H 2 O 2 concentration in space and time, both outside and inside the bacteria ([H 2 O 2 ] external , [H 2 O 2 ] cell , respectively).(C) Schematic of the core OxyR gene regulatory circuit corresponding to the stress response component of the model.It predicts the expression dynamics of proteins that scavenge intracellular H 2 O 2 (KatG, AhpCF) and control OxyR oxidation status (GrxA).Model output illustrated for constant [H 2 O 2 ] external exposure from 0 min.(D) The growth model describes the inhibition of cell elongation by H 2 O 2 .In turn, the growth dynamics feed into the stress response model by determining the dilution rate of enzymes.Top: cell elongation rate as a function of intracellular H 2 O 2 .Bottom: exponential growth and division cycles of a single cell without H 2 O 2 treatment.(E) Cell-cell interactions are described by a reaction-diffusion model where intracellular scavenging of H 2 O 2 creates a stress gradient from the edge to the interior of a cell population.Changes in the number and arrangement of cells in the population are determined by the growth model.See also Figure S1.

Figure 2 .
Figure 2. The model predicts chaos in the stress response (A) Oxidative stress response fluctuations in individual ''mother cells'' at the base of a one-dimensional population with ''barrier cells'' positioned closer to the H 2 O 2 source.The model produces seemingly random dynamics of GrxA protein expression level during continuous H 2 O 2 treatment from t = 0 min.The curves represent three independent simulation runs, starting with unsynchronized cells at random points in the cell cycles.(B) Stress response fluctuations diverge greatly over time, even if the differences in initial conditions are very small: here shown by the GrxA dynamics for 3 mother cells that differ very slightly in their initial stage of the cell cycle (2.5$10À4 % and 5$10 À4 % length differences).(C) (Left) Representative GrxA dynamics for a mother cell with continuous treatment at different H 2 O 2 concentrations (80, 140, and 260 mM).(Right) Histogram of counts of extrema detected for 3 mother cells for different H 2 O 2 concentrations.(D) Phase diagrams for the GrxA dynamics of the mother cells presented in (C), displaying bistable (80 mM) and multistable periodic oscillations (140 mM) and chaotic fluctuations (260 mM).(E) Bifurcation plot of the GrxA extrema values over a range of H 2 O 2 concentrations (n = 3 simulations per concentration).Vertical lines represent example traces in (C) and (D).(F) A positive Lyapunov exponent (l) shows chaotic divergence from initial conditions, computed for GrxA dynamics at different H 2 O 2 concentrations.Individual points represent single mother cells, with red dots for chaos (l>0) and black dots for periodicity (l % 0).Blue line and shaded region show mean ± SD of n = 3 cells simulated per H 2 O 2 concentration.(G) The autocorrelation function (ACF) distinguishes periodic and chaotic response fluctuations.Mean of ACF for GrxA of mother cells decreases steeply for chaotic traces under high H 2 O 2 treatment (1.7 mM, purple) and shows regular peaks for periodic traces under low H 2 O 2 treatment (140 mM, orange) (n = 3 simulations).See also Figure S2 and Video S1.

Figure 3 .
Figure 3. Chaos emerges in a general model of stress responses in cell populations (A) Uptake of toxins reduces cell growth rates (G) and triggers an intracellular stress response (S) that creates toxin gradients by cell-cell interactions (I).S-G-I feedback modulates toxin concentration in space and time, both outside and inside the bacteria ([Toxin] external , [Toxin] cell respectively).(B) Schematic of a generic stress response in which exposure to a toxin induces the expression of a detoxifying enzyme with rate K activation that removes toxin with rate K cat .(C) Model output illustrates the expression dynamics of the enzyme (maroon) and the intracellular toxin concentration (yellow) for constant external toxin exposure from t = 0 min without S-G-I feedback.(D) A positive Lyapunov exponent (l) shows chaotic divergence from initial conditions, computed for enzyme expression dynamics over a range of toxin concentrations (n = 3 simulations per toxin concentration).Higher external toxin concentrations lead to chaos.(E) Phase diagrams for the enzyme expression dynamics of a mother cell at toxin concentrations marked by vertical lines in (D), displaying periodic oscillations and chaotic fluctuations.(F) Higher K cat of the enzyme increases chaotic behavior.Lyapunov exponent for enzyme expression dynamics over a range of K cat values.(G) Phase diagrams for the enzyme expression dynamics of a mother cell at K cat values marked in (F), displaying periodic oscillations and chaotic fluctuations.(H) Higher expression rate K activation of the enzyme increases chaotic behavior.Lyapunov exponent for enzyme expression dynamics over a range of K activation values.(I) Phase diagrams for the enzyme expression dynamics of a mother cell at K activation values marked in (H), displaying periodic oscillations and chaotic fluctuations.
and S3G-S3L): (1) growth rates were reduced by switching to a lessfavored carbon source for E. coli (from glucose to glycerol) (Figure 6A; Video S4), (2) the protective effect of cell-cell interactions was reduced by manufacturing a modified mother machine with fewer cells in each growth channel (Figure 6B; Video S5), and (3) the strength of the stress response was reduced by lowering the concentration of H 2 O 2 (Figure 6C; Video S6).

Figure 4 .
Figure 4. Experiments on the oxidative stress response in E. coli reveal a good fit with the modeling predictions (A) Top: model simulation snapshot of GrxA expression after 90 min of 100 mM H 2 O 2 treatment.Bottom: snapshot of experiment with E. coli cells growing in a ''mother machine'' expressing PgrxA-SCFP3 after 90 min of 100 mM H 2 O 2 treatment.Scale bar, 10 mm.(B and C) Model simulation predictions and experimental data for mean GrxA expression (left) and mean cell elongation rates (right) under constant 100 mM H 2 O 2 treatment from t = 0 min for cells at different positions in growth trench (n = 100 simulated trenches and 3 experimental repeats).(D) PgrxA-SCFP3 dynamics of individual cells diverge greatly over time under constant 100 mM H 2 O 2 treatment from t = 0 min in experiments (5 representative mother cells shown).(E) The steep decay of the autocorrelation function (ACF) of the response fluctuations is consistent with chaos.Mean of ACF for PgrxA-SCFP3 of mother cells with 100 mM H 2 O 2 treatment (blue, 3 experimental repeats).ACF for individual cell traces shown in black (n = 100).See also Figure S3 and Video S2.

Figure 5 .Figure 6 .
Figure 5. Measurements of stress response dynamics in E. coli are consistent with chaos (A) Decision tree algorithm by Toker et al. 28 suggests that most mother cells in experiments display deterministic and chaotic response dynamics under 100 mM H 2 O 2 treatment.Pie-chart indicates the fraction of dead (red) and alive (green/black) mother cells detected as having stochastic (black/red) or deterministic (green) dynamics under 100 mM H 2 O 2 (n = 3,581 cells, 3 experimental repeats).(B) PgrxA traces (top) and their phase diagrams (bottom) for representative mother cell traces treated with 100 mM H 2 O 2 treatment from t = 0 min, which are classified as deterministic (green) or stochastic (red).(C) Bar plots show mean and standard deviation of maximal correlation dimension for experimental (black) and model (orange) GrxA traces of mother cell with (dark) or without shuffling (light) under 100 mM H 2 O 2 treatment, as computed by the Grassberger-Procaccia method.Random shuffling was performed as a control to remove temporal relation between data points.The low correlation dimension is consistent with determinism in experiments and simulations.(D) Stress response dynamics of cells growing in a colony are consistent with chaos.Snapshots of PgrxA-SCFP3 expression with 1 mM H 2 O 2 treatment from t = 0 min show cell-cell variability (scale bar, 10 mm).(E) Single-cell trajectories from the colony experiment in (D) are consistent with chaotic divergence of stress response dynamics.See also Figures S3, S4, S5, and S6 and Video S3.
1,806 and 1,991 cells, respectively; n R 3 repeats), (B) population size perturbation (1,003 and 1,440 cells, respectively; n R 3 repeats), and (C) H 2 O 2 concentration perturbation (1,003 and 1,361 cells, respectively; n R 3 repeats).(D) Bar plots show mean and standard deviation of maximal correlation dimension for GrxA traces of mother cells in experimental conditions shown in (A), (B), and (C) resulting in chaos (teal), or periodicity (black), as obtained from the Grassberger-Procaccia method.(E) The periods of non-chaotic oscillations correlate with cell cycle duration (interdivision time) over a range of growth rates in simulations and experiments.Mean and standard deviation of interdivision time and the time of the first ACF peak for simulated GrxA dynamics (left, n = 3 simulations per condition) and experiments (right) with 100 mM H 2 O 2 in M9 glycerol in 25-mm trenches (black), 50 mM H 2 O 2 in M9 glucose in 10-mm trenches (dark blue), and 25 mM H 2 O 2 in M9 glucose + 10% LB in 10-mm trenches (light blue) (966, 397, and 661 cells, respectively; n R 3 repeats).See also Figures S2 and S6 and Videos S4-S6.
Oxidative stress response modelThe oxidative stress response was modelled based on mass action kinetics using a set of 5 coupled ordinary differential equations (ODEs) to predict the dynamics of gene expression and intracellular H 2 O 2 concentration ([H 2 O 2 ] cell ) for a given external H 2 O 2 concentration ([H 2 O 2 ] external ).This is illustrated in the following schematic: [H 2 O 2 ] external by the outermost cell first which is exposed to the fixed concentration of [H 2 O 2 ] external in the growth media.This leads to a reduced external [H 2 O 2 ] external concentration for the cell located immediately beneath the outermost cell.The procedure is repeated to predict the external [H 2 O 2 ] external concentration from one barrier cell to the next until the mother cell is reached at the bottom of the population.Secondly, the stress response model predicts for each cell the changes in the intracellular [H 2 O 2 ] cell concentration, concentration of reduced regulator OxyR red , and concentrations of the stress response enzymes (GrxA, KatG, AhpCF) for the next time point (Equations 10, 11, 12, 13, and 14).It uses the external [H 2 O 2 ] external concentration from the cell-cell interaction model and the cell elongation rate g (Equation R grxA;basal + K grxA;act $ ½OxyR total À ½OxyR Red À ½OxyR total À ½OxyR Red Þ+h grxA;act The intracellular H 2 O 2 concentration is determined by the influx of external H 2 O 2 with rate R influx $½H 2 O 2 external , a basal endogenous production rate R H2O2;basal , and scavenging by catalase and peroxidase enzymes with Michaelis-Menten kinetics where K AhpC , K KatG are the catalytic rate constants and h AhpC , h KatG are the Michaelis constants: 60respectively.For the sigmoidal function describing the sensitivity of elongation rate to H 2 O 2 , the growth rate reduces to half when½H 2 O 2 cell = c 2 = 10 À 4 mM, with growth stalling at ½H 2 O 2 cell $ 2c 2 .The value c 2 was chosen in the sub-nanomolar range as E. coli cells can tolerate up to nanomolar ½H 2 O 2 cell in the absence of exogenous H 2 O 2 .60Thesteepness of the sigmoidal decay is given by parameter c 1 and was refined to match the variation of elongation rates for cells at different positions in a growth trench under H 2 O 2 treatment observed in the mother machine experiments.To this end, we varied c 1 over five order of magnitude from 200mM À 1 and minimized the value of mean absolute error (MAE) between g model and g experimental .For a given value of c 1 , the MAE was estimated for steady-state elongation rate values obtained for cells at different positions in trenches growing in M9 glucose and treated with various concentrations of H 2 O 2 (12.5 mM, 25 mM, 37.5 mM, 50 mM, 62.5 mM, 75 mM, 82.5 mM, 100 mM and 500 mM).Where i denotes the different H 2 O 2 concentrations and cell position denotes the position of cells in the trench from the open end.The MAE plotted against different values of c 1 showed a minimum at 2$10 4 mM À 1 [FiguresS3C and S3D].Cell-cell interaction model parametrisationThe reaction-diffusion Equation9was solved with a cell radius of R c = 0.575 mm and growth trench widthW = 1.2 mm.The diffusion coefficient of H 2 O 2 in water is D = 10 À 9 m 2 =s.The H 2 O 2 absorption rate60is k abs = 1:6 $10 À 5 m=s.Oxidative stress response model parametrisationThe model was parametrised using the following literature values:K Ahpc = 660s À 1 (Ref 61 ), h AhpC = 1:2 mM (Ref 61 ), K KatG = 490000s À 1 (Ref61 ), h KatG = 5900 mM (Ref 61 ), R H2O2;basal = 0:02 mM min À 1 (Ref 60,61 ), K ox = 0:1mM À 1 s À 1 (Ref 62 ), K red = 8mMs À 1 (Ref 62 ), h OxyR = 2583 mM (Ref 62 ), ½OxyR total =1 mM (Ref 62 ).The gene regulatory parameters were matched to the experimental data presented in this paper and in Choudhary et al. 41 : Red ½H 2 O 2 cell + K red ½GrxA ½OxyR total À ½OxyR Red À ½OxyR total À ½OxyR Red Þ+h OxyR + ɣ OxyR (Equation 22) K katG;act ½OxyR total À ½OxyR Red À ½OxyR total À ½OxyR Red Þ+h katG;act K ahpC;act ½OxyR total À ½OxyR Red À ½OxyR total À ½OxyR Red Þ+h ahpC;actTo solve the equations numerically using the Euler-Maruyama method, Equations 22, 23, 24, and 25 were converted to the form shown in Equation18, giving the following:½OxyR Red i+1 = ½OxyR Red i + Dt nd g À K ox ½OxyR Red i ½H 2 O 2 cell i + K red ½GrxA i ½OxyR total À ½OxyR Red i À ½OxyR total À ½OxyR Red i