A simple regulatory network coordinates a bacterial stress response in space and time

Bacteria employ diverse gene regulatory networks to protect themselves from stressful environments. While transcriptomics and proteomics show that the expression of different genes can shift strongly in response to stress, the underlying logic of large regulatory networks is difficult to understand from bulk measurements performed at discrete time points. As a result, it remains challenging to predict how these regulatory networks function at a system level. Here we use time-resolved single-cell imaging to explore the functioning of a key bacterial stress response: The Escherichia coli response to oxidative stress. Our work reveals a striking diversity in the expression dynamics of genes in the regulatory network, with differences in the timing, magnitude, and direction of expression changes. Nevertheless, we find that these patterns have a simple underlying logic. Firstly, all genes exhibit a transient increase in their protein levels simply due to the slowing down of cell growth under stress. Controlling for this effect reveals three classes of gene regulation driven by the transcription factor OxyR. Downregulated genes drop in expression level, while upregulated genes either show pulsatile expression that decays rapidly or gradual induction, dependent upon transcription factor binding dynamics. These classes appear to serve distinct functional roles in cell populations. Pulsatile genes are stress-sensitive and activate rapidly and transiently in a few cells, which provides an initial protection for cell groups. Gradually upregulated genes are less sensitive and induce more evenly generating a lasting protection that involves a larger number of cells. Our study shows how bacterial populations use simple regulatory principles to coordinate a stress response in space and time.


Introduction
Bacterial stress response regulons are highly complex in that they commonly involve a suite of genes performing diverse functions, which display a variety of expression patterns under stress 1 .However, these regulons are also typically under the control of a single master transcriptional regulator.This begs the question of how such simple regulation is able to orchestrate the timing and magnitude of so many genes.Bulk transcriptomic and proteomic studies have documented an association between stress tolerance and gene regulation for a wide range of stresses, including under starvation 2,3 , oxidative stress 4 , osmotic stress 5 , heat stress 5,6 , pH stress 7 , DNA damage 8,9 , phage infection 10 , and antibiotic exposure 11 .Moreover, these methods often identify tens to hundreds of genes that shift in expression with stress treatments, which suggests a complex regulatory network underlying each response 12 .
However, bulk population measurements are limited in their ability to link gene regulation with cellular phenotypes such as growth rate, because they lack single-cell resolution.In addition, due to cost and practicality, genome-wide gene expression data is typically recorded at only a few time points, meaning that temporal dynamics are not followed in any detail 13,14 .
An important alternative to bulk expression measurements is to use fluorescent reporter constructs that allow the gene expression in single cells to be followed.This approach has indicated the importance of temporal order in the activation of different genes within a regulon, e.g. for the SOS response, antibiotic stress, flagellar synthesis, and metabolic pathways 8,9,16,18,[21][22][23][24] .Spatial cell-cell variability can also emerge due to variability in the concentrations of stress agents and metabolites in cell populations 15,17,25 .The result can be complex spatio-temporal gene expression patterns under stress, which are shaped by cell-cell interactions, feedbacks between cell growth rate and gene expression 20,26,27 , and environmental fluctuations 19,20,28,29 .It is also becoming increasingly clear that stress survival is only partly determined by the protective responses of each individual cell, but strongly dependent on the collective protection provided by cell populations.However, to date, most single-cell studies have focussed on one or a few reporter genes for a given stress or disregarded the potential for cellular interactions, which leaves open the question of how whole stress response regulons are coordinated in cell populations.

Genes in the OxyR regulon show variable responses to stress
We imaged gene expression dynamics using time-lapse microscopy of E. coli cells, each carrying one of 31 plasmid-based transcriptional reporters that have an OxyR-regulated promoter followed by a fast-maturing GFP 61 or sCFP3 (for PkatG, PahpC and PgrxA 28,62 ).We collected data for ~10,000 individual cells growing inside microfluidic growth trenches ('mother machine type') during continuous H2O2 stress from ~2 hours before until 6 hours of treatment at 3-min time intervals, resulting in a rich dataset of ~200,000 expression measurements per reporter gene [Figure 1A].The microfluidic chip provides an opportunity to disentangle genuine gene regulatory dynamics from environmental changes because the continuous inflow of fresh media with H2O2 keeps the treatment concentration constant despite changes in the H2O2 scavenging capacities of cells [Figure 1B].Another feature of our approach is to study the AB1157 strain of E. coli that carries an amber mutation in RpoS 63 .This allows us to isolate the role of OxyR in gene regulation without the influence of stationary phase regulation of certain oxidative stress response genes by RpoS [64][65][66] .
With this setup, we were able to observe that all 31 genes showed an initial expression pulse within 10 min after the start of treatment [Figure 1C], but expression of different genes diverged after this initial pulse.This variability fits with observations on the expression dynamics of the oxidative stress response of Pseudomonas aeruginosa 67 .After 90 min of continuous H2O2 treatment, most genes remained upregulated while others became downregulated relative to the untreated expression level [Figure 1C].To further characterise differences in the regulation between genes, we computed the coefficient of variation (CV, variance normalised by the mean) across the expression levels of the 31 genes.The gene-togene variation was low in untreated cells, greatest during the transient expression peak shortly after treatment, and intermediate during steady-state with continuous H2O2 treatment [Figure 1D].These patterns suggest that the genes vary not only in expression magnitude but also show pronounced differences in induction timing and dynamics.

Transient growth inhibition causes a genome-wide expression pulse during stress
Our methodology enables us to analyse individual cell growth dynamics alongside gene expression.Importantly, this analysis revealed that the expression pulse at the start of H2O2 treatment precisely coincides with a period of transient growth inhibition [Figure 2A].During normal growth, protein concentrations stay constant in a cell because the expression rate is balanced with the dilution of molecules due to cell growth and division.H2O2 stress is known to rapidly stall cell growth which is triggered by an OxyS-mediated cell cycle arrest 68 .To remove the effect of growth dynamics on genes in the OxyR regulon, we normalised changes in expression level by the cell growth rate, to give a measure of the promoter activity 16 .
Correcting for growth rate effects revealed that 3 of the 31 genes (Phcp, PyfdL, PybjN) showed no significant change in promoter activity at any time [Figure 2B, C, Movie S1-S2].Moreover, another 9 of the genes (7 downregulated and 2 upregulated) showed only a transient change in regulation during the first ~100 minutes of treatment.These loci return to a baseline expression levels even though H2O2 is still present at a constant external concentration.
Focusing our characterisation on genes with lasting expression changes left 11 upregulated genes and 8 downregulated genes that showed a prolonged response to stress.
An interesting consequence of the effects of growth arrest on protein levels in the cell is that the genes that we identify as negatively regulated still show a positive pulse in protein levels after treatment [Figure 2B, C].In other words, the transcription rate from these promoters is downregulated in response to H2O2, but the reduction in growth rate outweighs this effect, resulting in an effective accumulation of proteins and the illusion of gene upregulation.These conclusions are supported by studying the constitutively-expressed Prna1 promoter that is not part of the OxyR regulon, which our analysis shows has a similar expression pulse but constant promoter activity [Figure 2C, Movie S3].That is, we find that the growth arrestdependent accumulation of proteins is a global genome-wide effect and not specific to oxidative stress response genes.

Upregulated genes display either pulsatile or gradual dynamics
We next focused on the 11 genes whose gene expression changes and remains upregulated with H2O2 treatment [Figure 2B, S1].Because we had found that gene-to-gene differences in expression were largest during the transient expression peak, we plotted the reporter expression level during the transient peak versus the level at steady-state [Figure 3D].This analysis reveals that the upregulated genes cluster on 2 distinct slopes.Seven of the upregulated genes (PkatG, PyaaA, PclpS, PhemH, PuxuA, PpoxB, PyaiA) all have a high ratio of peak/steady-state ~ 3, while the other upregulated genes (PgrxA, PtrxC, Pfur, PahpC) share a lower ratio of peak/steady-state ~ 1.5 [Figure 3E].The first class of genes showed a strong initial expression pulse that decays rapidly after its peak and stabilises at a low steady-state level (hence termed "pulsatile").In contrast, the second class of genes with lower peak/steady-state ratio showed a gradual induction to an elevated steady-state level (hence termed "gradually induced") [Figure 3B, C, Movie S1-S2].After classifying the up-regulated genes into two groups based on the peak/steady-state ratio, we also found that the pulsatile genes are generally activated earlier after the onset of H2O2 stress, achieving their maximal expression in 25.7±4.7 minutes compared to 37.5±7.8minutes post H2O2 treatment for the non-pulsatile genes [Figures 3F, S1, Movie S1-S2].

Modelling explains how differences in OxyR binding generate pulsatile or gradual dynamics
What explains the different patterns of upregulation?The oxidative stress response involves feedbacks from various metabolic pathways 4,32,36,41,48,49,[53][54][55] and several genes in the OxyR regulon are co-regulated by more than one transcription factor 36,[64][65][66]69,70 . Diffeences in the expression of genes in the OxyR regulon, therefore, may arise in multiple ways.However, we hypothesised that the differences we are seeing in the two classes of upregulated genes has a simpler underlying cause in the properties of OxyR itself.To explore this, we turned to a published minimal mathematical model from our previous work, which uses ordinary differential equations to describe changes in gene expression over time.Here, H2O2 permeates through the cell membrane and oxidises OxyR, which leads to the induction of genes of the oxidative stress response regulon 71 [Figure 4A].For simplicity, we reduced the operon to include only grxA (glutaredoxin-1) that converts OxyR back to its reduced form, the genes ahpCF (alkyl hydroperoxidase) and katG (catalase) that scavenge intracellular H2O2, and a 'reporter gene' that is induced by OxyR but does not have a functional role (to mimic the GFP reporters used in experiments).Gene expression increases on induction by OxyR and decreases due to dilution by cell growth.The growth rate slows down with increasing H2O2 concentration.
We previously showed that this model recapitulates many aspects of the gene regulatory response and growth dynamics of E. coli with H2O2 treatment 71 .Here, we employ the model to understand the basis for the different expression patterns seen across the genes in the OxyR regulon.We first tested the model's ability to explain the passive gene induction caused by the growth arrest-dependent accumulation of proteins.As predicted, when the expression of the reporter gene was uncoupled from OxyR regulation in the model, genes still show the expression pulse during the transient growth arrest while the promoter activity stayed constant [Figure 4B].In addition, for a reporter gene that is downregulated by OxyR, we again observe a similar induction pulse followed by a decrease in expression [Figure 4C].
To further explain the data, we extended the model to explicitly describe the effects of OxyRdependent gene regulation based upon two parameters: Kind captures the maximal expression rate when the promoter is fully occupied by oxidised OxyR, and KD captures the dissociation constant of oxidised OxyR from the promoter.Varying only these two parameters allowed us to shift expression dynamics in upregulated genes from gradual to pulsatile as observed in the data.Specifically, changing KD affects the peak/steady-state ratio of the gene, with pulsatile genes being characterised by a high KD value [Figure 4D].Kind meanwhile determines the position of the gene along the slope in a peak vs steady-state plot.Our model predicts, therefore, that the two categories of gene regulation, pulsatile and gradual are the result of distinct promoter dissociation constants of OxyR [Figure 4D, E].In support of this, experiments showed that the dissociation constant of oxidised OxyR from the pulsatile gene promoter PkatG is an order of magnitude higher than for the gradually induced PahpC promoter 72 [Figure 4F].We conclude, therefore, that the divergent patterns in gene regulation in the regulon can be explained by coupling a simple OxyR redox switch with patterns of cell growth under stress.

Pulsatile genes protect against sudden stress
What is the benefit of pulsatile expression of some genes and gradual induction of others?
The model shows how the dissociation constant of OxyR affects sensitivity of genes to changes in H2O2 concentration.Specifically, promoters with a high KD for OxyR show an approximately proportional increase in activity with H2O2 concentration, whereas the activity of promoters with a low KD is less sensitive to H2O2 [Figure 5A].In support of this prediction, experiments show that rising H2O2 concentration causes a much steeper increase for the pulsatile PkatG promoter activity compared to the gradually induced PgrxA promoter [Figure 5A, S2].
Similarly, other pulsatile genes showed a higher sensitivity to H2O2 compared to gradually induced genes [Figure 5B, S2].In Mycobacterium tuberculosis, induction of katG showed higher sensitivity for H2O2 than ahpC 73 , matching the pulsatile behaviour of katG and gradual induction of ahpC in E. coli.The model further predicts that the pulsatile gene induction is a consequence of a response to a transient spike in intracellular H2O2 which accumulates until the increased level of scavenging enzymes tips the balance [Figure 5C].In support of this explanation, experiments show that the expression pulse of PkatG is much reduced when cells respond to a gradually increasing dose of H2O2 as opposed to a step treatment [Figure 5D].These analyses suggest that pulsatile genes are important in response to sudden stress, whereas gradually induced genes are more important during prolonged stress.Consistent with this hypothesis, deletion of the pulsatile katG gene makes cells extremely sensitive to sudden H2O2 treatment, but these cells were still able to survive a gradually increasing dose of H2O2 reaching the same final concentration [Figure 5E, F, Movie S4].That is, the expression pulse of katG appears to be critical for the rapid production of catalase enzymes to counteract the transient spike in intracellular H2O2 but katG expression becomes dispensable after the initial adaptation delay, as noted before 62 .H2O2 tolerance during prolonged exposure is then provided by the gradually upregulated AhpCF alkyl hydroperoxidase, which scavenges lower concentrations of H2O2 efficiently 30 .

Pulsatile genes show higher cell-to-cell expression heterogeneity
We next asked to what extent the expression of pulsatile and gradually induced genes is coordinated within the same cell.Firstly, we found that the time of induction after the onset of stress was shorter and less variable across cells for the pulsatile genes [Figure 5G].Interestingly, for the control strain, the plasmid-expressed PgrxA displayed a slightly higher expression pulse compared to the PgrxA reporter on the chromosome, even after normalising for differences in the steady-state level [Figure 6C].This is a subtle but significant effect and explained by our model, which showed that the plasmid copy number increases transiently when the cell growth rate slows during the onset of H2O2 treatment [Figure 6C].This effective increase in gene dosage leads to an additional expression boost for the plasmid-based promoter.
Focusing on the dual reporter strain, we found that PkatG and PgrxA expression in the same cell showed relatively little correlation during the transient expression pulse [Figure 6D, Movie S5].However, both promoters exhibited substantial fluctuations in gene expression during prolonged treatment, and we found that these fluctuations were closely correlated [Figure 6E, Movie S5].This was evident from gene expression levels at discrete time points and from temporal cross-correlation analysis [Figure 6F].Together, our single-cell analysis revealed that pulsatile and gradually induced genes are tightly coordinated during prolonged stress, but appear to be uncoupled in regulation during the transient response to sudden stress.

Differentially-regulated genes show distinct spatial patterns in cell populations
Spatial patterns of gene expression are important in the oxidative stress response because of the ability of cells nearest the stress to protect cells that are further away.This phenomenon is remarkably effective, with each individual cell reducing the local H2O2 concentration by around 30% in its vicinity inside the one-dimensional microfluidic channel 26 .Within colonies, cells are protected in all three dimensions by surrounding cells leading to even steeper spatial H2O2 gradients from the edge to the interior of the population 28 .However, how this spatial structuring affects the whole oxidative stress regulon is unknown.We, therefore, asked whether the different temporal expression patterns we observe across the regulon were associated with different spatial expression patterns.
Quantifying expression levels in the channels of the microfluidic device revealed clear differences in the spatial patterning of regulation across genes.Downregulated genes exhibited an inverse gradient to that of the stress, with the lowest expression seen in the frontier cells that are closest to the H2O2 source [Figure 7A, C, Movie S2].Even passively induced genes that are not part of the OxyR regulon showed a spatial gradient in expression level that was driven by a more pronounced inhibition of growth for cells that are located closer to the H2O2 source [Figure 7C].Turning to upregulated genes, we see a steeper decay in expression in space as one moves away from the source of stress for pulsatile genes as compared to the gradually induced genes, in both our model and experiments [Figure 7, Movie S1].That is, pulsatile genes are strongly upregulated in relatively few cells that are close to the H2O2 source, which generates significant cross protection, such that only the most stressed cells show expression.The gradually upregulated genes then activate more evenly in larger numbers of cells to provide lasting protection.

Discussion
The rapid and coordinated regulation of stress response genes is critical for survival in a changing environment 74 .However, bulk measurements provide limited insights into how gene expression dynamics are coordinated in space and time within cell populations.To address this, we leveraged the power of single-cell transcriptional reporters to follow the bacterial oxidative stress response under continuous H2O2 treatment.Our work revealed a diverse set of dynamics across 31 oxidative stress response genes.Despite this diversity, our modelling and experiments show that OxyR binding dynamics together with the effects of stress on the cell growth rate can explain the different gene regulation patterns.Promoters with a negative Kind exhibit reduced expression rate, despite the fact that protein levels transiently increase due to a slow-down of cell growth.For upregulated genes, promoters with a high KD drive a strong transient expression pulse after H2O2 treatment leading to pulsatile expression, whilst promoters with a lower KD activate genes more slowly and sustain elevated expression rates throughout prolonged treatment.Moreover, by using a dual reporter, we show that the expression of different genes in the stress response diverges during the initial rapid response to treatment but then becomes closely correlated at steadystate.
The dynamics we observe are broadly consistent with previous bulk measurements 4,36,49 for time points shortly after the start of treatment [Table 1].We observe differences at later time points [Figure S4].However, these differences are to be expected as bulk methods typically study gene expression in dense cultures after a single high dose of H2O2.H2O2 levels decay quickly in bulk cultures and the stress response abates, which is not the case with continuous supply of H2O2 in microfluidic chips 4,37,50 .
What is the evolutionary function of the different patterns in gene regulation that we observe?Our model and experiments suggest that the higher H2O2 sensitivity of pulsatile genes provides a rapid defence against sudden bursts of intracellular H2O2.Interestingly, several of the pulsatile genes including yaaA 38 , clpS 37 and hemH 41 are involved in iron regulation.Their pulsatile induction is likely beneficial to counter the rapid lethality of H2O2 experienced due to the Fenton reaction 62,[75][76][77][78] .Being more sensitive to local H2O2 fluctuations, pulsatile genes therefore display high cell-cell variability in expression magnitude and steeper spatial gradients compared to gradually induced genes.Of the two H2O2 scavenging enzymes, pulsatile catalase (katG) is important to protect against immediate stress experienced by cells on the periphery of a population while gradually induced alkyl hydroperoxidase (ahpCF) efficiently scavenges low H2O2 levels that reach the rest of the population.Gradually induced thioredoxins (grxA 78 , trxC 46 ) regulate the redox status of proteins, including OxyR.These genes, therefore, show a more consistent response both in space and time, with more cells contributing to protection at the steady-state.In this way, the action of a single transcription factor is able to orchestrate the stress response to generate both individual and collective protection.
Data and code availability: All the raw data collected for analysis in this study along with the custom-built python codes for analysis and simulations will be freely and openly available on the Oxford Research Archive.Any further information about data and code is available upon request by the corresponding authors.

Strains and plasmids:
We used strains derived from E. coli K12 AB1157 for our experiments.All the strains had a constitutively expressed Prna1-mKate2 marker aiding the segmentation of cells during microscopy and flhD deletion to stop the cells from escaping out of the microfluidic growth channels during time-lapse microscopy.The reporter plasmids for genes used in this study were obtained from a transcriptional reporter library of PSC101 plasmids 83 .The reporter plasmids comprised of GFPmut2 fluorescence protein after the promoter region of a given gene or operon and had a kanamycin resistance marker.GFPmut2 was replaced with SCFP3A fluorescent protein for PkatG, PahpC, and PgrxA using Gibson Assembly (NEB) 28 .The reporter plasmids were mini-prepped and transformed in our AB1157 background strain.Strains were selected on 25 μg/mL kanamycin resistance LB agarose plates and checked for fluorescence signal by microscopy snapshots.
The chromosomal reporter for PgrxA-SCFP3A was constructed by inserting PgrxA-SCFP3Akanamycin using λ-red recombination in endogenous loci (aslA) on the chromosome of DH5α strain expressing pKD46-ampicillin.The successful colonies were selected by growing them overnight at 37˚C on kanamycin plates and subsequently re-streaking again on kanamycin selective plates to remove the temperature sensitive pKD46 plasmid.Then, the insert was moved into our strain of interest that had Prna1-mKate2 marker and flhD deletion, using P1phage transduction.The successful colonies were selected on kanamycin plates and tested for the insert using colony PCR, after which, the colonies were restreaked and grown overnight thrice to remove the phage.Colony PCR was again performed before storing the strain at -80˚C.PgrxA-mYPet and PkatG-mYPet plasmids with ampicillin resistance marker were obtained by genetically editing PgrxA-SCFP3A-kanamycin and PkatG-SCFP3A-kanamycin plasmids using Gibson assembly.We first replaced CFP with mYPet fluorescent protein and the resultant plasmid was edited again using Gibson assembly to replace kanamycin with an ampicillin resistance marker.The final strain was checked by colony PCR, microscopy snapshot, and re-streaking on 100 μg/mL ampicillin plates.Finally, the PgrxA-mYPetampicillin and PkatG-mYPet-ampicillin plasmids were transformed in the PgrxA-SCFP3Akanamycin chromosomal fluorescent reporter strain to construct the dual reporter strain.The strain was imaged to test for dual fluorescence using microscopy snapshots and selected on plates supplemented with 25 μg/mL kanamycin and 100 μg/mL ampicillin.

Media and growth conditions
Strains were stored in 20% glycerol stocks at -80˚C.Strains were streaked on freshly prepared LB agar plates supplemented with appropriate antibiotics for selection (25 μg/mL kanamycin and/or 100 μg/mL ampicillin) and incubated at 37˚C.One colony was picked from the overnight plates and added in minimal M9 media with appropriate antibiotics for overnight growth at 37˚C in a shaking incubator at 200rpm.The minimal media was prepared by mixing M9 salts (15 g/L KH2PO4, 64 g/L Na2HPO4, 2.5 g/L NaCl, and 5.0 g/L NH4Cl), 2 mM MgSO4, 0.1 mM CaCl2, 0.5 mg/mL thiamine, MEM amino acids, 0.1 mg/mL L-proline, and 0.2% glucose.
The next day, cells were diluted 1:50 in minimal M9 media and grown for ~3 hours for the microfluidic experiments.0.85 mg/mL Pluronic F127 was also added to this culture and minimal media to avoid cell aggregation when loading the cells on a microfluidic chip.Minimal media with or without H2O2 was flown continuously through microfluidic chips using syringe pumps.The concentration of H2O2 used was as specified in the figure legends or text, and was added to the minimal M9 media just before setting up the experiments.

Microfluidic chip preparation and setup
Microfluidic setup and preparation of 'mother machine' devices were performed as described in Choudhary et al 28 .

Time-lapse microscopy
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 Prna1-mKate2 (λ = 555 nm) , 300 ms for MutL-mYPet (λ = 508 nm), 75ms for GFPmut2 reporter (λ =470 nm) and 75 ms for sCFP3 reporter (λ = 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 and images were captured every 3 min.

Mother machine data processing and analysis
The microscopy time-lapse .nd2files for the experiments were processed using BACMMAN plugin 84 in Fiji 85 , which were subsequently analysed using custom made Python scripts.
BACMMAN first performs pre-processing which corrects for drifts during imaging and aligns growth channels spatially over time.To obtain cellular parameters, we used the Prna1-mKate2 marker to segment and mark the cell outlines against the background and then used this as a mask for other fluorescence channels.In all the fluorescence channels (CFP, GFP and mYPet), the mean fluorescence intensity inside the mask of each cell mask were computed.The cells were tracked over time using the segmentation masks of mKate2 channel to provide cell lineage information.BACMMAN generates output in separate excel files containing cell growth characteristic, Prna1-mKate2 intensity data and one for each fluorescence data.These files were then further analyzed using a custom python pipeline to compute cell parameters as described in the following section.

Cell parameter calculations
Elongation Rate: The instantaneous elongation rate at time  was calculated based on the logdifference in cell length   at time between consecutive frames as

Coefficient of variation:
The CV values were calculated as the standard deviation divided by the mean.
Relative expression in spatial dimension: Expression intensity of cells with atleast 1 barrier cells were divided by the intensity of outermost cell (i.e.cell with no barrier cells) for each growth channel at a given time point.

Cross-correlation analysis:
The temporal cross-correlation between the different reporter intensity traces of mother cells was computed using the statsmodel library in Python.
Correlation values from individual growth trenches were then averaged over all observed growth trenches.

Linear regression analysis
The linear regression analysis was performed in Python using the stats.linregressfunction in scipy library which computes the least square regression for a linear fit between two sets of data points.

Oxidative stress response model
We modelled the oxidative stress response model as shown in Figure 4A  This leads to the following equations:
A constitutive gene, that is not part of the OxyR regulon, was modelled such that it was expressed at rate   and diluted with growth: The However, the amplitude of the expression peak reached by individual cells was more variable for pulsatile genes compared to gradually induced genes [Figure5H, Figure S1, S3].Regulation of pulsatile genes, therefore, appears to prioritise precise control of the timing over the magnitude of the response [Figure 5G-H, S3].That is, this form of regulation function to activate a response very quickly, if not accurately, in the face of rapid onset stress.Next, we constructed a dual-reporter strain with two fluorescently labelled transcriptional reporters: the pulsatile PkatG-YFP reporter expressed from a plasmid and the gradually induced PgrxA-CFP reporter inserted on the chromosome [Figure 6A, Movie S5].We also constructed a control strain with two PgrxA reporters, one marked with CFP on the chromosome and one marked with YFP on the plasmid, to account for variability due to differences in fluorescent proteins (YFP versus CFP) and the difference in gene copy numbers of the reporters (chromosomal versus multi-copy plasmid) [Figure 6A, Movie S6].The mean expression dynamics of the dual reporter strains matched our previous analysis for PkatG and PgrxA in single reporter strains [Figure 6B, Movie S5].

Figure 2 :
Figure 2: Oxidative stress response genes show a range of different expression dynamics during constant H2O2 stress (A) Mean elongation rate (teal) and mean expression rate for all genes (black; dashed) of frontier cells treated with 100 μM H2O2 from t=0 min (shaded region).(B) Schematic illustrates the expression dynamics of the 5 gene regulation categories under H2O2 treatment (from t=0, shaded region) together with the number of genes in each category: (left) downregulation throughout the H2O2 treatment (green) or transient downregulation post treatment (gold); (middle) upregulation throughout the H2O2 treatment (purple) or transient upregulation post treatment (grey); (right) no significant change in promoter activity observed (black).Dashed lines represent the basal level of expression before treatment.(C) Mean promoter activity of frontier cells for 31 transcriptional reporters (ordered

Figure 3 .
Figure 3. Upregulated promoters show either pulsatile or gradual induction dynamics (A) Plot of the peak expression level (during the transient phase) versus the expression level at steady-state of frontier cells relative to basal expression for the indicated transcriptional reporters with 25 µM, 50 µM and 100 µM H2O2 treatment.The pulsatile (pink) and gradually induced (blue) genes cluster on two distinct slopes, shown with linear fits (dashed lines).(B) Peak/steady ratio obtained by linear regression of values for individual genes shown in panel A cluster into two categories (pink: pulsatile genes, blue: gradual induced genes, n≥3 repeats per gene).(C) Plot shows mean reporter expression levels of frontier cells for pulsatile (pink) and gradual (blue) genes under 100 µM H2O2 treatment provided at t = 0 min (shaded region) (n ≥ 1500 and n≥3 repeats per gene).(D) Kymograph of PhemH and PtrxC transcription reporter expression representative of pulsatile and gradual induction with 100 µM H2O2 from t=0 minutes, respectively (Scale bar = 10 µm).(E) Plot shows mean (bold) and single cell (thin) reporter expression levels of mother cells for PhemH (pink, pulsatile) and PtrxC (blue, gradual) under 100 µM H2O2 treatment provided at t = 0 min (shaded region).

Figure 4 .
Figure 4.A model of the oxidative stress response predicts the molecular basis of the different categories of gene regulation (A) Schematic represents the oxidative stress response model.The cells experience influx of H2O2 at Rinflux rate, causing OxyR oxidation and induction of stress response genes.GrxA (glutaredoxin-1) converts oxidised OxyR back to its reduced form and scavenging enzymes AhpC and KatG lower the intracellular H2O2cell concentration.The expression of the proteins is balanced by dilution due to cell growth, where the cell elongation rate g is a function of H2O2cell.OxyR also regulates a reporter gene which has no function in the response itself.(B) Effect of growth inhibition by H2O2 on gene expression dynamics.(left) Cell elongation rate affected by intracellular H2O2 concentration (orange, dashed) compared to constant growth

Figure 5 .
Figure 5. Pulsatile genes respond quickly to bridge the adaptation lag after H2O2 treatment (A) The response of pulsatile genes is more sensitive to changes in H2O2 concentration than gradually induced genes.Model outputs (left) and experimental data (right) for expression level at peak (bottom) and steady state (bottom) of PkatG (pulsatile gene, pink) and PgrxA (gradual, blue) across a range of external H2O2 concentrations.Dashed lines show linear fits.(B) Relative changes in gene expression at peak (top) and steady-state (bottom) of frontier cells with H2O2 treatment ranging from 25 µM to 100 µM H2O2.Pulsatile genes (pink) show higher dose-sensitivity compared to gradual induced genes (blue) (n≥2 repeats per gene per H2O2 concentration).(C) Model results of intracellular H2O2cell concentration with step

Figure 7 :
Figure 7: Spatio-temporal expression patterns across the OxyR regulon (A) Snapshots of cells with representative reporters for the three categories of gene regulation: PfhuF (downregulated), PyaaA (pulsatile) and PahpC (gradually induced) during peak expression (top) and at steady-state (bottom) with 100 µM H2O2 treatment.(B) Scavenging of H2O2 by bacteria creates H2O2external gradients in the growth trench from the source of treatment at the open end to the mother cell at the closed end.Model output for the expression of pulsatile (pink), gradually induced (blue) and downregulated (green) genes for cells across the growth trench (increasing number of barrier cells) relative to the expression of frontier cells.(C) Mean expression of 31 transcriptional reporters for cells with different number of barrier cells with 100 µM H2O2 treatment relative to the expression of frontier cells during peak expression (left) and steady-state (right) (pulsatile: pink, gradual: blue, down: green, constitutive Prna1: black, putative: grey; n≥3 repeats per gene).(D) Schematic depicting the spatial variation across a bacterial population in the expression of OxyR-controlled genes that are downregulated (left), pulsatile induced (center) and gradually induced (right) under H2O2 treatment.
log(  )− log( −t ) t .For calculating the elongation rates of cells at different positions in the growth trench, cells were tracked according to their initial position until the number of barrier cells decreased by 2. Gene expression: Reporter fluorescence intensity values were averaged over the area of each cell.Number of barrier cells: The number of cells that are located between the open end of a trench and the cell being analysed.Promoter activity: The instantaneous promoter activity at time  was calculated as the cell averaged rate of change of total fluorescence intensity of a cell.It was computed as 16 : is the mean fluorescent intensity of cell and A is the cell area.the instantaneous elongation rate of the cell.For calculating the promoter activity values of cells at different positions in the growth trench, cells were tracked according to their initial position until the number of barrier cells decreased by 2. Expression Rate: The instantaneous expression rate at time  was calculated based on the difference in cell intensity   at time between consecutive frames as   −  −t t .For calculating the expression rates of cells at different positions in the growth trench, cells were tracked according to their initial position until the number of barrier cells decreased by 2. Peaky expression: 90 th percentile of the mean fluorescence intensities from 12 to 90 minutes post H2O2 treatment.Steady expression: Average of the mean fluorescence intensities from 150 minutes post H2O2 treatment.Time to peak expression for single cells: Time in minutes post treatment until 100 minutes when individual cells show maximum fluorescence intensity.Peak expression for single cells: Mean fluorescence intensity for individual cells from 18 minutes until 90 minutes post treatment.This value was subtracted by the mean fluorescence intensity without H2O2 treatment i.e. basal fluorescence intensity.

Figure S3 :Figure S4 :
Figure S3: Cell-cell variability in gene expression magnitude for oxidative stress response genes in different regulation categories Coefficient of variation for peak (top) and steady state (bottom) gene expression values for frontier cells under 100 µM H2O2 treatment (pulsatile: red, gradual: blue, negative: green, putative: grey, constitutive Prna1: black).Error bars represent standard deviation (n ≥ 3 repeats per gene)

Table 1 : Genes chosen for this study along with their characteristics.
P: upregulation, N: downregulation; p-values using Mann-Whitney tests where n.s. for p>0.05 359 71d previously described in Choudhary et al71.Here, external H2O2 is provided at the rate   .Intracellular concentration of H2O2 ([ 2  2 ]  ) oxidises OxyR with rate   and converts it into oxidised form that induces multiple genes in its regulon.For simplicity, we modelled the induction of important genes of stress response regulon by OxyR: glutaredoxin-1 (grxA) that converts OxyR oxidised back to its reduced form with rate   , scavengers: alkyl hydroperoxidase (ahpC) and catalase (katG) that reduce [ 2  2 ]  and a reporter gene that does not directly or indirectly affect the OxyR regulation.All the genes are produced at a basal expression rate ( , ,  , ,  ℎ, ,  , ) and their expression is modulated upon induction by OxyR.The induction by OxyR is modelled as Michaelis Menten , ,  ℎ, ,  , ,  , ) and half-maximal induction Km is given by dissociation constant of OxyR from the promoter regions of the genes (KD) (  , ,  , ,  ,ℎ ,  , ).Finally, the gene expression reduces due to growth dependent dilution effect where growth rate () is a function of [ 2  2 ]  .
intracellular [ 2  2 ]  concentration is determined by the influx of external H2O2 with rate   • [ 2  2 ]  , a basal endogenous production rate   2  2 , , and scavenging by catalase and peroxidase enzymes with Michaelis-Menten kinetics where  ℎ ,   are the catalytic rate constants and ℎ ℎ , ℎ  are the Michaelis constants: Lastly, the effect of plasmid copy number variation is considered by modifying the equation for reporter gene expression (Eq 5) such that the gene induction is proportional to the plasmid copy number .The change in plasmid copy number is computed as a constitutive expression with rate   and diluted by cellular growth as given below.