Tunable phenotypic variability through an autoregulatory alternative sigma factor circuit

Genetically identical individuals in bacterial populations can display significant phenotypic variability. This variability can be functional, for example by allowing a fraction of stress prepared cells to survive an otherwise lethal stress. The optimal fraction of stress prepared cells depends on environmental conditions. However, how bacterial populations modulate their level of phenotypic variability remains unclear. Here we show that the alternative sigma factor σV circuit in B. subtilis generates functional phenotypic variability that can be tuned by stress level, environmental history, and genetic perturbations. Using single-cell time-lapse microscopy and microfluidics, we find the fraction of cells that immediately activate σV under lysozyme stress depends on stress level and on a memory of previous stress. Iteration between model and experiment reveals that this tunability can be explained by the autoregulatory feedback structure of the sigV operon. As predicted by the model, genetic perturbations to the operon also modulate the response variability. The conserved sigma-anti-sigma autoregulation motif is thus a simple mechanism for bacterial populations to modulate their heterogeneity based on their environment.


INTRODUCTION
Cells live in a changeable environment and experience a wide range of environmental stresses. Bacterial populations have evolved strategies to survive these stresses. One strategy is for all cells to immediately respond to stress with the activation of the relevant stress response pathway (Hilker et al. 2016) . Alternatively, a bacterial population can generate a broad range of cellular states, which allows it to hedge its bets against the changeable environment (Veening et al. 2008) . Noise in gene expression has been proposed as a mechanism for generating phenotypic variability in genetically identical cells (Raj and van Oudenaarden 2008;Martins and Locke 2015) . However, how the bacterial population regulates individual cell decisions in order to modulate the fraction of cells that enter an alternative transcriptional state remains unclear ( Figure 1 A).
The σ V mediated lysozyme stress response pathway in B. subtilis is an ideal model system to examine how bacterial populations can tune their phenotypic variability. σ V is an extracytoplasmic function (ECF) alternative sigma factor. Alternative sigma factors replace the 'housekeeping' sigma factor, σ A , in the RNA polymerase holoenzyme and redirect it to regulons that control distinct regulatory programs. They have already been shown to display a high level of gene expression variability in B. subtilis Young et al. 2013;Cabeen et al. 2017;Park et al. 2018) and the σ V activation pathway is both well characterised and specific to one stress condition, which greatly simplifies analysis of its activation. σ V is the only pathway activated in response to C-type lysozyme (Ho et al. 2011;Guariglia-Oropeza and Helmann 2011;Ho and Ellermeier 2012) . Lysozyme is produced by animals as part of their innate immune system, and kills bacteria by cleaving the peptidoglycan of the cell wall between the N-acetylmuramic acid residue and the time-lapse microscopy of fluorescent reporters for σ V activity to characterise σ V induction dynamics in individual cells in response to lysozyme stress. We found that upon induction by sub-lethal levels of lysozyme, σ V is activated heterogeneously, with some cells activating σ V rapidly, whereas some cell lineages do not activate σ V for multiple generations. This heterogeneity is functional, as cells that respond to an initial sublethal stress are more likely to survive a subsequent lethal stress of lysozyme. Through experiment and modelling, we found that these dynamics can be explained solely by the autoregulatory feedback of σ V on its own operon. Our model predicted that the observed heterogeneity could be tuned by environmental history and by genetic perturbations, which we confirmed experimentally. The conserved sigma-anti-sigma autoregulation motif is thus a simple mechanism for bacterial populations to tune their levels of phenotypic variability.

Individual cells activate σ V heterogeneously under lysozyme stress
To characterise σ V activation dynamics, we first constructed a B. subtilis reporter strain containing a chromosomally integrated fluorescent reporter for σ V activity, P sigV -YFP. This strain also contained a reporter for the housekeeping sigma factor σ A (P trpE -RFP). P trpE -RFP expression was used as a constitutive control and to aid image analysis (Methods). We then used time-lapse microscopy to examine σ V activity in individual cells grown in the mother machine microfluidic device (Wang et al. 2010) . This device allows long-term tracking of hundreds of individual cells trapped at the ends of channels. It also allows fast switching between media conditions during imaging (Methods). We first measured σ V expression dynamics in response to a step change to a sub-lethal concentration of 1 µg/ml lysozyme  Figure S2 and Movie S1). The first cells begin to respond to the lysozyme stress with raised P sigV -YFP levels within 2 frames (20 minutes), suggesting rapid activation given the ~15 minute maturation time of the YFP reporter protein ( Figure S3 ). However, our time-lapse movies revealed that P sigV -YFP was heterogeneously activated at the single-cell level ( Figure 1 C). While 10% of cells reached the half maximum of σ V activity within 70 min of being exposed to lysozyme, it took 200 min, (approximately 4 generations) for 90% of all cells to activate σ V (  To test whether the observed heterogeneous activation of σ V in response to lysozyme was due to the growth conditions in the mother machine, we investigated σ V activation dynamics in liquid culture ( Figure S5 ) and using an alternative microfluidic device ( Figure S6 ). Both conditions showed similar heterogeneous σ V activation to that observed in the mother machine. Therefore the observed heterogeneity in P sigV -YFP is independent of the experimental setup. Previous work has shown that heterogeneous gene expression could be due to intrinsic noise in the expression of the P sigV -YFP reporter (Elowitz et al. 2002) , rather than due to heterogeneous σ v activity. To test whether the observed heterogeneity was due to the intrinsic variability of the P sigV -YFP reporter we constructed a strain containing chromosomally integrated P sigV -YFP and P sigV -CFP reporters and exposed it to 1 µg/ml lysozyme in liquid culture. In these snapshot experiments, the expression of P sigV -YFP and P sigV -CFP was highly correlated (R 2 = 0.85) ( Figure S7 ). Thus the observed heterogeneity in P sigV -YFP reflects changes in σ v activity and not intrinsic variability of the P sigV -YFP promoter.
Next, we asked whether the observed heterogeneity in P sigV -YFP is modulated by the level of lysozyme applied. We examined P sigV -YFP expression after application of 0.5, 1, 2, and 4 µg/ml lysozyme. These values were all below the previously reported minimal growth inhibitory concentration of 6.25 µg/ml (Ho et al. 2011) . To measure the distribution of σ V activation times, for each time point we calculated the fraction of cells that had crossed the half-maximum of their respective final σ V value (meaning the cell had activated σ V ). We found that when increasing the lysozyme concentration from 0.5 µg/ml to 4 µg/ml the heterogeneity in σ V activity was reduced ( Figure 1 E, Figure S2 and Figure S8 ). The time for 90% of cells to activate their σ V pathway decreased from 300 min (approximately 6 generations) for 0.5 µg/ml to 100 min (approximately 2 generations) for 4 µg/ml lysozyme ( Figure S4 ). At the same time, the fold change in induction between the unstressed σ V activity and the steady state σ V activity under lysozyme stress increased from ~190 for 0.5 μg/ml to ~520 for 4 μg/ml ( Figure 1 F and Figure S1 ). We also observed that under a 4 µg/ml concentration of lysozyme some cells (8% and 21% in two different repeat experiments) appeared sick and were wider than usual cells. These cells also overshot their new σ V activity steady state before relaxing to it. We removed these cells from our analysis ( Figure   S9 ), although including them did not affect our results ( Figure S10 ). Taken together, our results reveal that the level of phenotypic variability in sigV expression is tuned by stress levels, with P sigV -YFP heterogeneity reduced as lysozyme levels increase.
We next asked whether differences in individual cell states before applying a lysozyme stress could predict how a cell responds to the stress. There was no correlation between either growth rate or cell size before lysozyme application and response time after lysozyme application ( Figure S11 and Figure S12 ). Alternative sigma factors in B. subtilis display gene expression variability even in the absence of stress Park et al. 2018) . We found that, in the absence of stress, the coefficient of variation of P sigV -YFP was over 4 times higher than that of the constitutive control (0.65±0.29 versus 0.14±0.05) ( Figure S13 ). Cells which had elevated levels of P sigV -YFP before the addition of stress were more likely to activate σ V instantaneously ( Figure S14 and Figure S15 ). However, some cells with low P sigV -YFP levels could also activate σ V instantaneously ( Figure S14 and Figure S15 ). Therefore, while noise in P sigV -YFP expression before stress contributes to the observed heterogeneity in σ V activation after the addition of lysozyme, it is not the only cause. Next, we examined whether activating σ V had any effect on the survival against future lethal concentrations of lysozyme. Cells were first exposed to a priming stress of 1 μg/ml of lysozyme for 30 min, which induced σ V heterogeneously. Thereafter cells were exposed to 20 μg/ml of lysozyme for 20 min before removing lysozyme from the media ( Figure 2 A).
Surviving cells were defined as cells which had not lysed at the end of the movie, 280 min after the addition of 20 μg/ml of lysozyme. We found that the short exposure to sub-lethal concentrations of lysozyme (during the priming stress) improved survival to subsequent lethal stress levels ( Figure 2 B) from 0% ( Figure S16 A) to 12 ± 5% survivors ( Figure S16 B).
Surviving cells had on average 1.57 ± 0.33 fold higher P sigV -YFP levels than perishing cells immediately before the application of the lethal concentration of lysozyme ( Figure 2 C).

RsiV
We next attempted to understand how the single-cell activation dynamics of the σ V pathway are generated . First, to test whether the heterogeneity that we observed in σ V activation times was due to some cells activating a different stress response pathway, we analysed the genome-wide transcriptional response of cells to 1 μg/ml lysozyme in the wild type and in the Δ sigV knockout. As previously reported, only the sigV operon was induced by lysozyme (Guariglia-Oropeza and Helmann 2011) in the wild type ( Figure 3 B). The Δ sigV strain did not show any induction ( Figure 3 C). This suggests that no other pathways are induced by lysozyme. We therefore hypothesized that the observed P sigV -YFP heterogeneity is due to the σ V circuit itself (Ho et al. 2011;Hastie et al. 2013;Hastie et al. 2014) .
To test which components of the σ V circuit play a role in the heterogeneous activation of σ V , we constructed IPTG inducible overexpression constructs of all genes in the sigV operon ( sigV , rsiV, oatA, yrhK ) and in the σ V circuit ( sipS, rasP ) ( Figure 3 A). We then investigated σ V

Mathematical modelling reveals that P sigV -YFP expression dynamics can be explained by the 'mixed' σ V autoregulatory feedback loop
Based on our transcriptome and overexpression analysis we constructed a model of a simplified sigV operon consisting of sigV and its anti-sigma factor rsiV, with positive autoregulation of the operon by active σ V ( Figure 3 F). We did not include oatA , as it is not required for the heterogeneous activation of σ V ( Figure 3 F). RsiV binds σ V to form an inactive complex σ V -RsiV. On stress application, RsiV is degraded and σ V is free to activate the operon. The model was implemented using mass action kinetics, with hill equation terms for production rates. Chemical Langevin equations were used to model the intrinsic noise of the pathway, creating a system of stochastic differential equations (See methods for further details, (Gillespie 2000) ). Using a coarse-grained parameter search, we were able to find parameters that capture the heterogeneous σ V activation in response to lysozyme stress ( Figure 4 A and Figure S18 ). We were also able to recreate several aspects of the experimental data, including the dependency of induction time ( Figure 4 B) and steady state levels of sigV expression on increasing levels of lysozyme stress ( Figure 4 C). In the absence of lysozyme stress, the system is monostable and only a state with low σ V activity is accessible ( Figure S19 ).
Application of stress introduces a second steady state of high σ V activity, whilst also making the low sigV state only quasi-stable ( Figure S19 B). Noise in sigV expression causes a stochastic switch to the high σ V state. The stochastic nature of the switch generates the broad distribution of induction times observed ( Figure S18 ).

Phenotypic variability is tuned by doubling copy numbers of sigV operon genes
To test the assumptions of our model, we predicted the effect on σ V activation dynamics of introducing a second copy of each component of the sigV operon ( sigV, rsiV , or both sigV and rsiV ), with each component driven by the sigV promoter. Our model predicted that these perturbations would modulate the heterogeneity of the system's response to lysozyme stress. A second copy of rsiV meant that larger fluctuations in sigV expression were required to kick the system into the high σ V expression state ( Figure 5 B and Figure S20 ). This increased the cell to cell variability in response times to lysozyme stress, as compared to the WT in the simulations ( Figure S20 ). Conversely, with a second copy of sigV , or a second copy of both sigV and rsiV, smaller fluctuations in sigV expression are required to kick the system into the high σ V expression state. Cell to cell variability in response times to lysozyme therefore decreases in the simulations ( Figure 5 C and Figure S20 ). In fact, a second copy of sigV caused an increase in the levels of σ V expression even before the addition of lysozyme in the simulations ( Figure S20 ). To test these predictions, we constructed strains containing a second copy of either sigV , rsiV or sigV-rsiV driven by the sigV promoter. As predicted by our model, snapshots of sigV expression after induction by 1 μg/ml of lysozyme revealed that a second copy of sigV or sigV-rsiV reduced the observed heterogeneity in σ V activation, whilst a second copy of r siV caused an increase in the heterogeneity in σ V activation ( Figure 5 C). Finally, adding a second copy of sigV alone caused the activation of σ V even in the absence of stress, as predicted. These results validate our model assumptions, and also show that the heterogeneity in sigV expression is easily tuned by simple genetic perturbations.

The heterogeneous activation dynamics of sigV are dependent on environmental history
Given that the σ V activation dynamics appear sensitive to the baseline levels of σ V , RsiV, and the σ V -RsiV complex, any perturbations to these levels should affect σ V activation times. We hypothesized that cells should have elevated levels of sigV operon components after a lysozyme stress is removed due to the recent high expression of the operon components.
The elevated levels of σ V stored in stabilized σ V -RsiV complexes should cause all cells to respond immediately to a subsequent reapplication of lysozyme. We investigated this hypothesis in our model. We simulated the application of lysozyme stress, which was then removed for the equivalent duration of several cell cycles The same level of stress was then reapplied and the heterogeneity in σ V activation disappeared ( Figure 6 A and Figure S22 ).
However, as the delay before the second application of lysozyme was increased (allowing system components to relax to pre-stress levels), the heterogeneity gradually returned ( Figure 6 B&C and Figure S22 ). Thus our simulations predict the existence of a temporary molecular memory of the environmental history. To test whether the σ V circuit exhibits such a memory, we grew cells in the mother machine device under 1 μg/ml lysozyme stress before removing the stress for 6h (approximately 7 cell cycles). All cells stopped activating σ V as soon as lysozyme was removed, as indicated by the decay in YFP fluorescence (

DISCUSSION
Here, we report a general mechanism for a bacterial population to tune its phenotypic variability based on stress levels, genetic architecture, and environmental history. Using quantitative single-cell microscopy and microfluidics we found that the activation of the alternative sigma factor σ V in response to lysozyme stress was heterogeneous. While some cells activated their σ V pathway immediately, others could take up to six generations to activate their σ V pathway. The observed phenotypic variability plays a functional role, as cells that respond to a sub-lethal stress were more likely to survive a subsequent higher stress application ( Figure 2 ). Through experiments and modelling we found that this heterogeneity could be understood by the 'mixed' positive and negative feedback of σ V activating both itself and its anti-sigma factor RsiV. Alternative sigma factors are a common regulatory system in prokaryotes, often controlling stress response and virulence pathways. The 'mixed' feedback loop is also a common motif, suggesting that this motif can be a general mechanism for bacterial populations to tune their phenotypic variability.
Our modelling and experiments found that recent exposure to lysozyme stress is remembered, even though the system turns off after removal of lysozyme stress. The key to this behaviour appears to be the 'mixed' feedback loop, as it allows amplified levels of inactive σ V -RsiV complexes in each cell after stress. These complexes can be cleaved by a subsequent addition of lysozyme, releasing σ V to activate its operon. Similar memories of previous stress have been observed in bacterial systems, although typically not to modulate phenotypic heterogeneity. For example other pathways such as the heat stress response in B. subtilis (Runde et al. 2014) or the oxidative stress response in yeast (Kelley and Ideker 2009) have been shown to have a memory of past conditions. Often this memory is facilitated through an autoregulatory positive feedback loop that can lock the cell into an ON state that is heritable through cell divisions (Lambert and Kussell 2014;Novick andWeiner 1957) (Zacharioudakis et al. 2007;Acar et al. 2005;Hashimoto et al. 2013;Biggar and Crabtree 2001) . However, it is difficult for a single positive feedback loop to allow the system to be OFF but primed for future stress, as we find to be the case for the 'mixed' feedback loop in the σ V pathway. In the future, it will be important to investigate whether the 'mixed' feedback loop also tunes the levels of phenotypic diversity by environmental history in other systems. Interestingly, computational studies have proposed that a 'mixed' feedback loop We have found that a combination of noise in gene expression and a mixed feedback loop can generate tunable phenotypic diversity in an alternative sigma factor circuit. Going forward, it will be important to observe whether alternative sigma factors more generally are used as a mechanism for bacterial populations to tune phenotypic diversity. This is particularly the case given that alternative sigma factors often control pathways critical to pathogenicity and resistance to antibiotics. Noise in sigma factor activity has already been observed in multiple alternative sigma factor circuits in B. subtilis Park et al. 2018) , as well as for the general stress response sigma factor in E. coli (Patange et al. 2018) . Additionally, an alternative sigma factor plays a role in generating phenotypic diversity in Mycobacteria (Ghosh et al. 2011;Sureka et al. 2008) . Two other obvious candidates for further study are the pathogens Enterococcus faecalis and Clostridioides difficile (Ho and Ellermeier 2019) , which also have a σ V pathway that is responsive to lysozyme.

Strains and Media
All strains are derivatives of the PY79 background strain (Table S1). Deletions were generated by replacing genes of interest with an antibiotic resistance cassette by recombination of a linear DNA fragment homologous to the region of interest. All strains had a Δ ytvA ::neo deletion insertion. YtvA is a blue light sensor and it was deleted to avoid any activation of the general stress response pathway σ B by the microscope illumination (Locke  Gaidenko et al. 2006) . For strains used in mother machine experiments, cells were made immotile by inserting a Δ haG :erm deletion cassette in order to improve the loading into the microfluidic device. Additionally, all strains contained a house keeping σ A promoter-driven mCherry for segmentation and as a constitutive control.  .
Cells were routinely grown in Spizizen's Minimal Media (SMM) (Spizizen 1958) . It contained 50 µg/ml Tryptophan as an amino acid source and 0.5% glucose as a carbon source.
Cultures were started from frozen stock in SMM and grown overnight at 30°C to an OD between 0.3 and 0.8. The overnight cultures were resuspended to an OD of 0.01 and regrown to an OD of 0.1 at 37°C.

Plasmids
E. coli strain DH5α was used to clone all plasmids. The cloning was done with a combination of non-ligase dependent cloning and standard molecular cloning techniques using Clonetech In-Fusion Advantage PCR Cloning kits. Plasmids were chromosomally integrated into the PY79 background via double crossover using standard techniques. The list below provides a description of the used plasmids, with details on selection marker and integration position/cassette given at the beginning. Note that all plasmids below replicate in E. coli but not in B. subtilis .

ppsB::P trpE -mCherry ErmR
This plasmid was used to provide uniform expression of mCherry from a σ A -dependent promoter, enabling automatic image segmentation (cell identification) in time-lapse movie analysis. A minimal σ A promoter from the trpE gene was cloned into a vector with ppsB homology regions ) . The original integration vector was a gift from A. Eldar (Eldar et al. 2009) . The target promoter of sigV was cloned into the EcoRI/BamHI sites of AEC12 (Eldar et al. 2009) (gift from M. Elowitz, CalTech).

amyE::P sigV -mturq SpectR
The target promoter of sigV was cloned into the EcoRI/Nhe1 sites of plasmid amyE ::3XCFP Spect R . Target promoter sequences are described below.

A Nikon inverted Ti-E microscope (Nikon, Amsterdam, Netherlands) with a Nikon Perfect
Focus System (PFS) hardware auto-focus was used. All images were acquired using either a Photometrics CoolSnap HQ2 CCD camera (Photometrics, Tucson, AZ, USA) or a Photometrics Prime sCMOS camera (Photometrics, Tucson, AZ, USA). The microscope stage was enclosed in an incubator (Solent Scientific, Segensworth, UK) which was set to 37°C for all experiments. Illumination was provided by a LED lamp (CoolLED, Andover, UK).
Epifluoresence was provided by a Lumencore Solar II light engine (Lumencore, Beaverton, OR, USA). Chroma filters (Chroma, Bellows Falls, USA) #41027 for the RFP channel, #49003 for the YFP channel and #49001 for the CFP channel were used. All experiments were done with a phase 100x Plan Apo (NA 1.4) objective (Nikon, Amsterdam, Netherlands).

CellAsic Experiments
Cultures were grown overnight to an OD of 0.1, and then resuspended to an OD of 0.001 for loading into CellAsic B04 microfluidic chips (Merck, Darmstadt, Germany). Cells were loaded into the microfluidic chips with a pressure of 4 to 6.5 psi for 2 to 4 s. The following lysozyme concentrations were investigated with the CellAsic setup: 0 µg/ml, 0.5 µg/ml, 1 µg/ml and 4 µg/ml. Fresh media was perfused into the chip with a pressure of 1 psi. After 60 min of growth in standard SMM, the media was switched to media containing lysozyme. Cells were imaged at regular intervals (every 10 min) and the acquired movies were analysed with the standard Schnitzcells package for MATLAB ) .

Wafer Fabrication
We used two different microfluidic designs of the mother machine : 1. All the mother machine data in this paper were acquired with a microfluidic design based on the original mother machine (Wang et al. 2010 Figure S26 a different mother machine design was used. It was based on a mother machine design first described by Norman et al. (Norman et al. 2013) . The main difference to the Jun lab design is that the channels in which cells grew were longer and that they have a shallow side channel for better perfusion of media to the cells. As a substrate for the PDMS chip fabrication we used an epoxy in the fume hood and were ready to use the next day (Gruenberger et al. 2013) .

Experimental Preparation
The chips were plasma bonded to a glass bottom dish (HBSt-5040, Willco Wells, Amsterdam Netherlands) by treating the chip and the glass dish surfaces with a plasma (Femto Plasma System, Diener, Ebhausen, Germany) for 12 s at power of 30 W and a pressure of 0.35 mbar. The chips were then put into the oven for 10 min at 65°C, in order to strengthen the bonding. In the meantime, a bovine serum albumin (BSA) passivation buffer was prepared at a concentration of 20 mg/ml in SMM. The chips were then passivated with this buffer for 1 h at 37°C.

Memory Experiments
For all shown memory experiments, the JLB221 strain (ΔespH) was used instead of JLB130, unless stated otherwise. This was done in order to prevent clogging of the microfluidic chip during these long term experiments due to the production of extracellular matrix . The two strains had the same P sigV -YFP dynamics in response to lysozyme ( Figure S23 & Figure   S24 ). Cells were inoculated in SMM and grown overnight to an OD of 0.3 to 0.8. In the morning, the cells were resuspended to an OD of 0.01 in SMM + 1 μg/ml lysozyme and grown to an OD of 0.1. Cells were then loaded as described in the "experimental preparation" section.
In the mother machine the cells were grown for 4h with 1 μg/ml lysozyme before removing the stress. The stress was then reapplied after 6h or 12h.

Cumulative Activation Time
The activation time was calculated as the time between switching from SMM to lysozyme and the time point when then P sigV -YFP passed its half maximum ( Figure S1 ). The cumulative fraction was then calculated as the cumulative fraction of cells with P sigV -YFP values higher than the half maximum of their final values (representing cells that have activated) ( Figure S1 and Figure S8 ).

Removal of Overshooting Cells
For 4 μg/ml lysozyme the P sigV -YFP activity overshot its steady state activity. Overshooting cells were sick and also wider. Thus overshooting cells could be removed based on their width. First the single-cell width traces were smoothed with a gaussian filter and their maximum value was determined. All single cell traces with a maximum width larger than the mean width + 6 sigma of cells before the stress were removed ( Figure S9 & Figure S10 ).

Growth Rate
The instantaneous growth rate as shown in Figure S11 was calculated as in (Martins et al. 2018) : where len i is the length of the cell at frame i and Δt is the time between subsequent frames taken from the image metadata. Growth rates were only calculated within one cell cycle and not over a division.

Calculation of surviving cells in priming experiment
The fraction of primed cells which survived the 20 µg/ml lysozyme stress treatment in Figure   2 were defined as followed: Where Survivors is the fraction of cells that survive the 20 µg/ml lysozyme stress treatment.
min(#cells survivor_channel ) is the smallest number of cells in a channel during the movie for channels that still have at least one cell at the end of the movie (i.e to estimate the number of surviving cells we calculated the smallest number of cells in the channel between first frame after addition of 20 µg/ml lysozyme and the end of the movie). #cells before_stress is the number of cells in a channel in the last frame before addition of 20 µg/ml lysozyme stress.
Statistical analysis of the effects of higher P sigV -YFP levels before lysozyme stress application on response time As the distributions in Figure S14 were not normal we use the non-parametric Kolmogorov-Smirnov test (ks-test) to verify the null hypothesis that cells with higher P sigV -YFP values before lysozyme stress application would have the same activation time distribution as cells with lower P sigV -YFP expression. Cells with high P sigV -YFP were defined as cells with a YFP fluorescence higher than the mean YFP fluorescence before the addition of lysozyme plus one standard deviation. All other cells were defined as low P sigV -YFP. A p-value below 0.05 was used to reject the null hypothesis. The ks-test was performed using the built in MATLAB function kstest2 .
We found that for 0.5, 1 and 2 μg/ml lysozyme the null hypothesis was rejected. Thus cells with higher P sigV -YFP values before lysozyme stress application had a shorter activation time. For 4 μg/ml lysozyme the null hypothesis was accepted. Therefore cells with higher P sigV -YFP values before lysozyme stress application had the same activation time as cells with lower P sigV -YFP values. was transferred to a new microcentrifuge tube, to which 500 µl of absolute ethanol was added. After vortexing, the lysate was transferred to a RNeasy spin column and centrifuged 15 seconds at > 10000 rpm. RNA purification was then carried out following the instructions of the Quiagen RNeasy kit. RNA quality and integrity were assessed on the Agilent 2200 TapeStation, and RNA concentration was assessed using Qubit RNA HS assay kit. Library preparation was performed using ScriptSeq™ Complete Kit (Illumina, BB1224), for 2 µg of high integrity total RNA (RIN>8). The libraries were sequenced on a NextSeq500 using paired-end sequencing of 75 bp in length.

RNA-seq analysis
The raw reads were analysed using a combination of publicly available software and in-house scripts. We first assessed the quality of reads using FastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were aligned to the Bacillus Salmon (Patro et al. 2017) . Read counts for each gene were imported using the tximport R package (Soneson et al. 2015) . Genes differentially expressed between the WT and ΔsigV mutant or in response to lysozyme treatment were identified using the DESeq2 R package (Love et al., 2014).

Mathematical Model
We constructed a minimal mathematical model of the σ V circuit that was based on the core components that we found experimentally to modulate sigV heterogeneity ( Figure 3 ). The model consists of σ V , RsiV, and the σ V -RsiV complex. Free σ V that is not bound in a complex activates the production of both σ V and RsiV. Lysozyme stress was modelled as activating the cleavage of the σ V -RsiV complex, releasing σ V while degrading RsiV. We note that although heterogeneous activation can be simulated through just a positive feedback loop, our model is the simplest model that can reproduce the perturbations of rsiV and sigV that we observe experimentally.

Model reactions
We used the following set of biochemical reactions to model the dynamics of σ V , RsiV, and σ V -RsiV. Production: The production rate of σ V and RsiV was assumed to follow a hill function, where the production rate of both σ V and RsiV were identical and activated by σ V .
We assumed a constant and identical dilution rate for all three components of the system. Binding/Dissociation: We assumed that σ V and RsiV would bind to each other to form a complex and that this complex could dissociate. We also assumed the binding rate to be much higher than the dissociation rate. Cleavage: We assumed that the σ V /RsiV complex would be cleaved by lysozyme, with σ V getting released and RsiV degraded.

Model equations
The model was implemented using mass action kinetics, with a hill function representing the production rates. This generated the following system of ordinary differential equations: (Eq3) A stochastic version of the model was implemented using the chemical Langevin equations.
This generated the following system of stochastic differential equations:  Figure S19B) only exhibits a bistable type behaviour in a finite region, we note that for these parameter values, in the deterministic bifurcation diagram bistability remains as the stress levels approach infinity ( Figure S19A).
Thus, a deterministic version of the model with these parameters is unable to turn on, no matter the level of lysozyme stress, as stochastic fluctuations in the expression of the sigV operon are required for activation. Using a modified parameter set, we verified that heterogeneous activation dynamics are still observed in the stochastic model, even when the corresponding deterministic bifurcation diagram has a finite bistability region ( Figure S19 ).

Code availability
Code used for simulations can be found here: https://gitlab.com/slcu/teamJL/schwall_etal_2020. Code for analysis of data reported in this study is available upon request from the corresponding author.   shuts off σ V activation. Increasing stress levels from 1 µg/ml to 10 µg/ml recaptures the WT 1 2 3 4 5 σ V activation dynamics. n>=3 for all data shown. Only three randomly selected boxplots are shown for circuit components where n>3. (F) Schematic of simplified σ V circuit with only σ V and RsiV, where σ V (green) activates its own expression and that of its anti-sigma RsiV (orange, R).  simulations predict that a second copy of P sigV -sigV (2x sigV ) or P sigV -sigVrsiV (2x sigV-rsiV ) makes the σ V response to lysozyme homogenous while a second copy of P sigV -rsiV (2x rsiV ) should increase the heterogeneity. Bars represent the fraction of cells which have activated σ V 30 min after treatment with lysozyme. Cells which had activated σ V were defined as cells which had passed a threshold of the WT mean before lysozyme application plus six standard deviations. Error bars correspond to mean (N=999) ± SD (1000 bootstraps). (C) As predicted, adding a second copy of P sigV -sigV (2x sigV ) or P sigV -sigVrsiV (2x sigV-rsiV ) makes the σ V response to 1 µg/ml lysozyme homogenous in experiments, while a second copy of P sigV -rsiV (2x rsiV ) increased the heterogeneity. Each bar plot is the average of three biological repeats. The error bars correspond to the mean ± SD.