Exploiting fluctuations in gene expression to detect causal interactions between genes

Consider two genes of interest, geneX with transcript X and geneZ with transcript $Z_k$. Our approach relies on engineering a passive reporter gene geneY with identical transcriptional control to geneX and mRNA lifetime, but with transcript that does not affect other cellular components, see Figure 2A. This could be achieved with a transcriptional reporter that is put under control of the same promoter as geneX and placed at a similar gene locus. One way to ensure that geneY interacts minimally with other cellular components would be to remove its start codon so that it will not be translated into a protein.

Co-regulated mRNA reporters have been engineered to satisfy the assumptions underlying Equation 1, with transcripts counted with smFISH (Baudrimont et al., 2019; Raj et al., 2006; Skinner et al., 2013). Alternatively, future improvements in RNAseq (Saliba et al., 2014; Lubeck et al., 2014; Svensson et al., 2018) accuracy might allow for reliable fluctuation measurements of transcript levels in single-cells through sequencing approaches.

The dual reporter assay has been frequently implemented as co-regulated fluorescent proteins (Schmiedel et al., 2015; Pedraza and Paulsson, 2008; Elowitz et al., 2002; Bar-Even et al., 2006). The invariant of Equation 2 directly applies when X and Y are co-regulated fluorescent proteins with first-order translation rates and maturation times (see Appendix 1 Section 4B for details of the proof). Fluorescent proteins can then be used to detect causal interactions in gene regulation as follows. The gene of interest geneX is fused to a fluorescent protein to make a functional fusion geneX-FP. A passive reporter protein $Y$ is made by introducing a spectrally distinguishable fluorescent protein Y under the control of the same (but distinct) promoter as geneX-FP. The expression level $z_k$ of geneZ can be measured either through a transcriptional reporter or a functional fusion protein with a third fluorescence reporter. Under the assumption that the fluorescent protein Y does not directly affect other cellular components, Equation 2 directly applies to the covariances of fluorescence levels as long as X does not causally affect $Z_k$ . As the normalized covariances in Equation 2 are independent of scaling factors, standard fluorescence microscopy methods can be used without the need to determine absolute numbers. Two stable fluorescent proteins with similar maturation times (Balleza et al., 2018) should be chosen to ensure that the assumptions underlying the class of systems are satisfied. Additionally, the translation rates of X and Y need not be identical but can be proportional with a fluctuating proportionality factor as defined in the transitions of Equation 1. As a result, different fluorescent proteins with different ribosome binding sites, mRNA secondary structures, and gene lengths can be used as long as the translation rates remain proportional.

The class of stochastic processes presented in Equation 1 constrains the dynamics of cellular abundances in stationary processes. However, experimental data typically report measurements of growing and dividing cells. Under the assumption that during cell division, molecular abundances are divided on average proportional to cell size, we can show (see Appendix 1) that Equation 2 must be satisfied by the dual reporter abundances when $Z_k$ is a component not affected by geneX. In this analysis of cellular growth and division, we allow for partitioning noise, division time fluctuations, asymmetric divisions, and arbitrary growth-rate dynamics (see Figure 2B).

Chemical reaction rates typically depend on molecular concentrations and not absolute abundance numbers. Under general assumptions (see Appendix 1 Section 12), we can show that the covariance constraint of Equation 2 describes molecular concentrations in growing and dividing cells, see Figure 2. This result assumes that the abundance of X does not affect cell volume or growth rate. This requirement can be intuitively understood, because if X affects cell volume, then it causally affects the concentration of $Z_k$ which depends on volume (see Appendix 1 Section 12).

The class of systems defined in Equation 1 assumes that the dual reporters are degraded in a first-order process with a constant rate parameter $\beta$, see Equation 1. This assumption can be relaxed: the invariant of Equation 2 generalizes to the class of systems in which the degradation rate constant is an arbitrary function of all the cellular components that are not affected by the dual reporters, as long as it does not decay to zero (see Appendix 1 Section 4). As a result, the degradation constant can vary in time with arbitrary extrinsic fluctuations. Additionally, for the class of co-regulated fluorescent proteins, the translation rate parameters can also vary with arbitrary extrinsic fluctuations.

Under the null hypothesis that there is no causal interaction from X to Z, Equation 2 is equivalent to $r=1$, where $r={\eta }_{\mathit{xz}}∕{\eta }_{\mathit{yz}}$ is the covariability ratio. Given 95% confidence intervals for an experimentally measured covariability ratio, we can conclude at the 2.5% significance level that there is a causal interaction from X to Z if $r=1$ falls outside of the 95% confidence interval of the data (see Appendix 1—figure 2 for an example).

In order to illustrate the generality of the above results and to perform a numerical sanity check, we simulated several stochastic birth-death systems in growing and dividing cells, see Figure 2. We modeled 10 systems made up of 6 network topologies in which a component Z is correlated with X and Y but is not affected by $X$ , see Figure 2A. These systems include cellular components with non-linear reaction rates and feedback loops, fluctuating degradation rates, and confounding variables that affect all observed components. We analyzed each system subject to three different cellular growth dynamics: periodic exponential growth, periodic linear growth, and linear growth with fluctuating cell-division times and division sizes, see Figure 2B.

To numerically test Equation 2, we generated random sample paths for abundances and concentrations using the Gillespie algorithm (Gillespie, 1977; Voliotis et al., 2016). System parameters were varied over several orders of magnitude. The numerical data verify that the normalized covariances of cellular abundances as well as cellular concentrations satisfy Equation 2 in growing and dividing cells, see Figure 2C. The same does not hold for the corresponding Pearson correlation coefficients, see Figure 2C. Note, due to finite sampling, any numerical simulation will necessarily show small deviations from the exact equality of Equation 2. Through statistical analyses and re-running systems for increased sampling, we confirmed that the (minuscule) deviations observed in numerical simulations were consistent with finite sampling; see Materials and methods.

In Appendix 1 Section 14 we simulate a larger example network made of 10 components in which X regulates a cascade of components. The results suggest that the relative distance down a regulatory cascade can be inferred from the degree of violation of Equation 2.

Experimental techniques to measure mRNA abundances or fluorescence levels potentially introduce significant measurement noise. For example, when counting mRNA abundances, smFISH can lead to probabilistic undercounting noise and fluorescence microscopy can report fluorescence levels that include photon noise, read noise, and segmentation errors, while flow cytometry introduces Poisson noise when used with bacteria (Galbusera et al., 2020).

The invariant of Equation 2 holds in the face of measurement noise as long as the noise is symmetric in X and Y measurements. We can show (Appendix 1 Section 16) that Equation 2 holds in systems with arbitrary multiplicative noise that is independent of the X and Y signals, and additive noise that is independent of $Z_k$. Additionally, the invariant holds in the face of systematic undercounting that introduces a binomial readout of the signal of interest, along with Poisson-Gaussian noise. If, however, the experiment introduces a different type of noise in X as compared to Y, then Equation 2 is no longer valid. It is thus important to choose similar measurement techniques while measuring the abundances of X and its passive reporter Y.

Inherent to our fluctuation approach is the direction of inference, that is violations of Equation 2 imply the presence of a causal connection, whereas agreement with Equation 2 does not imply the absence of causal interactions. The question thus becomes whether a biologically relevant genetic circuit with known causal connections ever breaks our covariance relation.

We engineered four synthetic circuits in which a TetR-YFP fusion protein (X) represses an RFP reporter (Z), and a passive CFP reporter (Y) is under the same transcriptional control as X, see Figure 3A. While the circuits differ in dynamics and connections, $X$ affects $Z$ in all cases. Appendix 1—figure 17 for details of the synthetic circuits #1–4. Using time-lapse fluorescent microscopy with cells growing in a microfluidic device (see Figure 3C), we measured the normalized covariances from populations of genetically identical cells with high accuracy while discarding all temporal information.

Figure 3

Download asset Open asset

We found that circuits #1 and #2 clearly broke the constraint, circuit #4 did not violate it, and circuit #3 showed significant deviations right at the limit of what can be reliably detected, see Figure 3D. Using the estimated error bars and a null hypothesis test at the 2.5% significance level, we find that the data from circuit #3 does not reject the null hypothesis of Equation 2, see Appendix 1—figure 2.

The above repressilator test circuits were chosen for their well-characterized genetic interactions and because oscillating abundances make for a convenient test of the constructs. However, large oscillations in Z in also make it much more challenging to detect violations of Equation 2 because, in those circuits, the intrinsic noise we exploit is small compared to the periodically varying dynamics of the circuit. Indeed, the two clear violators correspond to altered versions of the repressilator circuit that do not oscillate (see Appendix 1—figure 17). These non-oscillating test cases may in fact be the most relevant assessment of using Equation 2 in natural gene regulatory circuits. All the above circuits were made of genes encoded on a low copy number plasmid. We have not shown that our method can successfully detect causal interactions in a chromosomally encoded gene regulatory circuit. However, plasmid copy number fluctuations will reduce the degree of violation of Equation 2 by reducing the fraction of reporter variability that is due to the intrinsic noise we exploit (see Appendix 1 Section 15). A priori, causal connections between chromosomally encoded genes are thus expected to be easier rather than harder to detect with our method.

Overall, the experimental data confirm that Equation 2 has the fundamental power to experimentally detect causal interactions. However, the precision of current experimental techniques lies at the edge of what is necessary to detect physiologically relevant interactions when genes show large variability relative to the intrinsic noise we exploit. This highlights the importance of correctly estimating measurement uncertainty in our approach. Our error bars estimate sampling error and errors from non-even sample illumination, temporal drift, background fluorescence, and autofluorescence. These were estimated using a bootstrapping approach where the corrections and the normalized covariances were computed recursively to samples of the data, see Materials and methods, Appendix 1—figure 6.

Our crucial theoretical result is that Equation 2 must be satisfied by any circuit in which X does not affect Z. To experimentally test this prediction, we engineered test cases in which X and Y are co-regulated, while Z is a component that is not affected by X. This took one of two forms: First, we constructed eight different circuits in which RFP (Z) was expressed constitutively using the pRNA1 promoter as previously used as a segmentation marker in E. coli (Potvin-Trottier et al., 2016; Lord et al., 2019). We used various synthetic circuits, and the RFP reporter was either integrated chromosomally or included on the circuit plasmid (the latter leading to correlations between RFP and our synthetic components without introducing a causal interaction). Second, as an additional control, for each circuit, we also considered the cellular growth rate as the unaffected cellular property Z. While the growth rate affects reporter concentrations through dilution, we did not expect cellular growth to be affected by our synthetic circuits, which made growth rates a convenient additional negative control.

We found the data from all but one synthetic circuit consistent with our predicted invariant; see Figure 3E. Note, the false positives correspond to the same synthetic circuit: one data point is from taking Z as the RFP concentration, while the other data point is from taking Z as the growth rate. To rule out experimental error, we repeated the experiment twice while swapping the YFP and CFP reporters, which confirmed the result (see Appendix 1—figure 7). Additionally, we confirmed the result using a different cell segmentation pipeline (see Appendix 1—figure 8).

Taken at face-value, our method thus produced two false positives (out of 16 tests), suggesting that our approach could be valuable but somewhat unreliable. However, next, we present experimental evidence that indicates the violating data points were not false positives but may have correctly identified an unexpected causal effect from X on Z in the outlier circuit.

Three cellular growth dynamics were simulated. The first and second are linear and exponential growth, respectively, with constant division times and symmetric cell divisions. Here, $V(t)$ trajectories were generated analytically as periodic functions. The volume is reduced by a factor $1∕2$ at evenly spaced division times $\{{\tau }_i\}$ . Between division times ${\tau }_i$ , ${\tau }_{i+1}$ the volume is given by $V(t)=V_0\cdot (1+\frac{t-{\tau }_i}{t_d})$ for the linear case and $V(t)=V_0\cdot 2^{(t-{\tau }_i)/t_d}$ for the exponential case, where $t_d$ is the time between divisions, and $V_0$ is the volume right after a division. The third simulated volume dynamics is linear growth with stochastic division times and asymmetric divisions. Here, a constant linear growth rate $\frac{dV}{dt}=u$ is used. Division times $\{{\tau }_i\}$ and division factors $\{a_i\}$ are picked recursively: The ${\tau }_{i+1}$ is taken from a normal distribution with mean set to the doubling time:

where $\mathcal{N}(\mu ,\sigma )$ is the normal distribution with mean μ and standard deviation $\sigma$ .Until then the volume grows linearly with constant rate

At ${\tau }_{i+1}$ the volume is reduced by factor $a_{i+1}$ taken from a normal distribution with mean set to the ratio of cell volumes at the beginning and end of the cycle

As a result, cells that grow more (less) than double in size during a cycle tend to divide with larger (smaller) division factors, ensuring the volume trajectories do not eventually expand (decay) to infinity (zero).

At division, molecular abundances are reduced according to a binomial splitting with probability given by the division factor. For example, if the volume is reduced by a factor of 0.4 at a cell division, that is $V\to 0.4V$, then each molecule has probability 0.4 to remain in the followed daughter cell.

Normalized covariances were computed by integrating over the trajectories to obtain time averages for first and second moments. This is equivalent to using the distribution given by calculating the fraction of the total system time spent in each sampled state. In the ergodic regime, this distribution converges to the stationary distribution of the ensemble.

We simulated ten groups of systems defined by their rate functions (see Appendix 1—table 1 for details). For each group, each model parameter was picked randomly multiple times from the set {0.1, 1, 10}. This was done for each of the three volume dynamics we considered. If a simulated system gave an average abundance in one of the components less than 0.01 molecules, then the system was omitted to avoid numerical errors that arise from divisions of small numbers when computing the normalized covariances. In total, there are 6522 resulting simulations which are plotted in Figure 2C.

For each system, simulated trajectories were generated forty times, and then repeated until the percentage errors of the normalized covariances, taken as the standard error of the mean divided by the mean, reached less than 1% or the number of simulations reached 1000. Each trajectory ran for 20,000 cell divisions and started with a unique random number generator seed. Each component abundance was set to 1 at the beginning of each trajectory, and the cell cycle time was set to 0. We let the simulated cells reach 200 cell divisions before we start to compute the time average integrals for the moments, in order for the effect of the initial condition to dissipate, and we analyze the systems’ cyclo-stationarity state. The final normalized covariance for each system corresponds to the average taken over the ensemble of simulated trajectories with 95% confidence intervals given by twice the standard error. Distributions for the uncertainties for each class of systems are plotted in Appendix 1—figures 3 and 4.

For a given system, Equation 2 was tested by verifying that the ratio ${\eta }_{xz}/{\eta }_{yz}$ produced a 95% confidence interval that encompassed the predicted value of 1. The test was satisfied by 6369 out of 6522 systems (97.7%), consistent with the definition of confidence intervals for finite sampling. The percentage of outliers for each simulated process is plotted in Appendix 1—figure 5. Re-simulating the outliers for twice the number of simulations led to a reduction in the standard error by a factor of $1/\sqrt{2}$, with Equation 2 being satisfied by 143 of the 153 outliers (93.5%). Re-simulating the remaining 10 violators with four times the number of simulations reduced the errors by a factor of 1/2, and all of the 10 outliers satisfied Equation 2 when quantified at such high numerical accuracy.

Polydimethylsiloxane (PDMS; Sylgard 184 Silicon Elastomer, Thermo Fisher Scientific) was mixed at a 10:1 (monomer $:$ curing agent) ratio, poured on top of a 1.0 μm tall wafer and degassed for 1 hr at room temperature before baking for an additional 1.5 hr at 65°C. After careful removal of the PDMS from the wafer, individual PDMS chips were cut out with a razor blade. The inlet and outlet holes were punched with a 0.75 mm biopsy puncher (World Precision Instruments). The PDMS chips were sonicated in isopropyl alcohol (Thermo Fisher Scientific) for 30 min and dried at 65°C for 15 min. Glass coverslips (Thermo Fisher Scientific: 22x40 mm #1.5) were cleaned with 1M potassium hydroxide (KOH, Sigma Aldrich) for 20 min. The PDMS chips were bonded to the glass coverslips using a plasma cleaner (Oxygen flow rate at 45 sccm, power at 30W for 15 s, Tergeo Plasma Cleaner, PIE Scientific). The completed microfluidic device was heated to reinforce the plasma bonding at 100°C for 10 min, then 65°C for 30 min.

E. coli strains were grown overnight in LB with appropriate antibiotics to select for cells containing the constructed plasmids. At ∼ 3-4 hr prior to the experiment start time, the overnight cultures were diluted 1:100 in imaging media consisting of M9 salts, 10% (v/v) LB, 0.2% (w/v) glucose, 2 mM MgSO4, 0.1 mM CaCl2, 1.5 μM thiamine hydrochloride and 0.85 g/L Pluronic F-108 (Sigma Aldrich, surfactant to prevent cells sticking to the surface of a microfluidic device). At an OD600 of 0.2-0.4, cells were loaded into the main feeding channel of the microfluidic chip and centrifuged at 5000 × g for 10 min to push cells into the cell trenches. The feeding channels were then connected to syringes filled with imaging media using Tygon tubing, and the microfluidic chip was placed in a temperature-controlled incubation chamber set at 37°C. Media was pumped through the feeding channel using syringe pumps (New Era Pump System), first at a rate of 5 μL/min for ∼ 0.5-1 hr to get the cells comfortable in the trenches, then at a high rate of 100 μL/min for 1 hr to clear the inlets and outlets. The rate of media flow was then set back to 5 μL/min for the duration of the experiment.

Images were acquired using a Zeiss Axio Observer inverted microscope equipped with a 63x Plan-Apochromat M27 oil objective (NA 1.40), an Orca Flash 4.0 LT camera (Hamamatsu), and an LED epifluorescence illuminator (Zeiss Colibri 7). The experiments were performed inside a temperature-controlled incubation chamber set at 37°>C. To reduce photobleaching, the exposure time (100 ms) and light intensity (10-20%) were set low, with 16-bit CZI images taken every 5-8 min. Focal drift was corrected automatically with the Definite Focus 2 (Zeiss) monitoring and compensation system using an infrared laser (850 nm). In all experiments in which CFP, YFP, and RFP are measured simultaneously (i.e. all experiments except for the one shown in Figure 4D), the following filter sets were used for acquisition: CFP (Semrock FF02-475/20-25), YFP (Chroma CT560/39bp), RFP (Zeiss Filter Set 91 HE LED), with respective Zeiss Colibri 7 illumination wavelengths set to 430 nm (CFP), 511 nm (YFP), and 590 nm (RFP), along with the Zeiss TBS 450/538/610 beam splitter. For the experiment in Fig. 4D in which a gfp reporter is used to measure gadX expression, the following filter sets were used for acquisition: GFP (Chroma ET525/50m) and cy5 (Zeiss BP 690/50) with Zeiss Colibri 7 illumination wavelengths set to 475 nm (GFP) and 630 nm (cy5), along with Zeiss QBS 405/493/575/653 beam splitter.

Because all three fluorescence channels were used to measure synthetic circuit components throughout, we used the bright-field channel for segmentation. This was achieved with the automated deep learning-based cell segmentation software DeLTA (Lugagne et al., 2020; O’Connor et al., 2022). The U-Net deep learning models used for channel identification and cell segmentation were trained with a large dataset built from 5 mother machine experiments performed in the Potvin Lab that were generated with the same microscopy and image acquisition setup. In these training datasets, a bright RFP is expressed, and a thresholding segmentation pipeline is performed on the RFP channel images to obtain segmentation masks. The U-Net deep learning models from DeLTA are trained with these segmentation masks in combination with bright-field images from the training dataset as described in Lugagne et al., 2020; O’Connor et al., 2022. The training dataset was segmented with the same thresholding pipeline used in previously described procedures (Potvin-Trottier et al., 2016; Norman et al., 2013). We also trained a DelTA model to use the RFP channel for segmentation and obtained similar results on the strains with bright RFP fluorescence. In Appendix 1, we include videos of a mother machine with segmentation boundaries obtained from our DeLTA models, along with the single-cell growth trajectories.

We estimated the YFP, CFP, and RFP single-cell concentrations as the average fluorescence intensities of all pixels in the segmentation mask of each cell. Cell area is measured as returned by opencv’s contourArea() over the segmentation mask. Cell length is measured by fitting a rotated bounding box to the segmented cell. We only analyzed data from the mother cells trapped at the top of the growth chambers.

For a given cell chamber with the Region of Interest identified by the U-Net model, a single-cell trace was constructed by selecting the segmented cell in each frame with highest average vertical pixel location. This, in effect, keeps track of the mother cells trapped at the top of the growth chambers.

Each single-cell trace went through a manual purging process. Traces with cell areas that are not growing and dividing were removed, as they are expected to correspond to dead cells. Traces with growing and dividing cell areas with very low constant fluorescence correspond to cells that lost the inserted plasmid due to random partitioning at cell division. These traces were also purged, but were analyzed separately to measure the cell autofluorescence as described in the next section. Moreover, segmentation and tracking errors were reduced by purging any parts of the traces that exhibited a clear non-biological anomaly.

Cell divisions were identified by a sudden decrease in the cell area. A division is called whenever the cell lengths dropped to less than 80% of their previous value.

Focal drift was reduced automatically with a 850 nm infrared laser (Zeiss Definite Focus 2). However, we used an oil objective, which can cause temporal drift from spreading of the oil over time. We thus applied an additional correction to the data as follows.

In a given mother machine experiment, we computed the mean YFP, CFP, and RFP signals of all the cells at each time frame. For example, the mean YFP time-trace is given by

where ${\text{YFP}}_{avg}^i(t)$ is the average YFP intensity taken over the segmented area of the $i^{\mathit{th}}$ cell at time frame $t$, and $N_{\text{cells}}(t)$ is the number of surviving cells at time frame $t$ . The mean time-traces of each fluorophore are fitted with a second-degree polynomial according to a least-squares fit. To correct for the drift, we multiply a measured intensity with the reciprocal of the respective fitted polynomial (see Appendix 1—figure 13). For example, if $f(t)$ is the obtained fit of $⟨{\text{YFP}}_{avg}{⟩}_t$ from Equation 3, we correct all the ${\text{YFP}}_{\mathit{avg}}$ measurements as follows

where $f(0)$ is the fitted function $f(t)$ taken over the first time frame.

A single image frame from our setup encompasses ∼30 cell trenches with total spatial extension of ∼ 100 μm. This results in uneven illumination onto the cells. We use the linear gain model (Peng et al., 2017; Smith et al., 2015) to correct for uneven illumination and background fluorescence:

where $I^m(x)$ and $I^{true}(x)$ are the measured and true intensities at horizontal pixel position $x$ respectively, the multiplicative term $S(x)$ models the uneven illumination, and the additive term $D(x)$ is any background noise that is present when no light is incident on the sensor. For our imaging setup, we decompose $I^{\mathit{true}}(x)$ as follows

where $I^{\mathit{FP}}(x)$ is the intensity from the fluorescent proteins, $I^a$ is the autofluorescence of the cell, $I^{media}$ is any fluorescence originating from the media, and $I^b$ is any fluorescence originating from the PDMS background of the microfluidic device. The total background is given by $b(x)=I^b(x)\times S(x)+D(x)$ , in which case the imaging model becomes

We use a ‘prospective’ approach to correct for $S(x)$ and $b(x)$ as follows. For each segmented cell in an image, we compute the average intensity in a box of dimensions 50 by 50 pixels located 15 pixels above the cell trench (see Appendix 1—figure 15A). This in effect estimates $b(x)$ at the $x$ position of just above each cell trench. In each image, this $b(x)$ is removed from the $I^{meas}(x)$ measurements of each segmented cell. To estimate the uneven illumination gain function $S(x)$, we pool all of the single cell fluorescence measurements according to their horizontal pixel position $x$. Data is binned into bins of size 50 pixels and averages are taken at each bin (see Appendix 1—figure 15B). The resulting curve gives an estimate of $⟨I^{FP}+I^{af}+I^{media}⟩\times S(x)$ . To estimate $S(x)$, we fit the curve with a 3rd degree polynomial $p_S(x)$ . To correct for the uneven illumination gain function, we multiply the fluorescence measurements by the reciprocal of $p_S(x)$ :

where $p_{S,max}$ is the max value of the fit $p_S(x)$ over $x$.

The preceding section described how the background and the uneven illumination were corrected using the linear gain imaging model of Equation 6. Even with the background correction, the fluorescence profile across growth chambers empty of cells is not negligible compared to the profiles of some cell-containing chambers, see Appendix 1—figure 15. This indicates that media fluorescence is not negligible in experiments with strains that produce low levels of fluorescence. Here, we show how we estimated the autofluorescence $I^a$ and the media fluorescence $I^{\mathit{media}}$ to obtain the sought-after fluorescent protein fluorescence $I^{\mathit{FP}}$.

In each mother machine experiment, there was a subset of cells that lost the synthetic circuit plasmid due to random partitioning at cell division. These cells were used to estimate $I^a+I^{media}$ . Cells that lost the plasmid were identified manually by observing single-cell traces: When a cell loses the plasmid, the CFP and YFP fluorescence decay rapidly and reach a constant low signal level for the remainder of the experiment, see Appendix 1—figure 14. The segments of the time traces at constant low signal level were cut and saved, with the autofluorescence and the media fluorescence taken as the average fluorescence of the saved traces. Around 10–30 cases of plasmid loss occur for each strain in a mother machine experiment.

Note that $I^a$ and $I^{\mathit{media}}$ do not need to be corrected in the RFP channel measurements. This is because RFP is taken as the component either affected or not affected by X. For the systems in which RFP is not affected by X, any function of RFP (like adding autofluorescence and media fluorescence) will also not be affected by X and the invariant of Equation 2 will be satisfied. Alternatively, if X affects RFP, then it will generally also affect any typically expected function of RFP. We did not correct for $I^a$ and $I^{\mathit{media}}$ in the RFP channel because a fraction of circuits have RFP located on the chromosome and not the plasmid.

The corrections from the preceding sections rely on the distribution of single-cell measurements. As a result, sampling error affects the accuracy of the corrections, along with the final estimators of the normalized covariances. We applied two methods to estimate normalized covariance confidence intervals by taking into account sampling error and its effect on the corrections.

First, we use bootstrapping, where the data corrections and the normalized covariances are computed over many samples of the data, allowing for replacement in sampling. If an experiment produces N single-cell traces of cells with plasmid and M single-cell traces of cells that lost the plasmid, random samples of size N and M of each respective ensemble are taken. The corrections from the previous sections are then applied using the samples, with the corrected data pooled into a distribution from which the normalized covariances are computed. This is repeated until 100 normalized covariances ${\eta }_{\mathit{sample}}$ have been computed for each sample. The final reported normalized covariance corresponds to the average over the ${\eta }_{\mathit{sample}}$ , with confidence intervals given by twice the standard deviation. See Appendix 1—figure 16 and Appendix 1 Section 17 for details.

Second, we use sampling without replacement, where the single-cell traces from each experiment are divided into 7 to 10 disjoint sets, and the corrections and the normalized covariances are computed for each set. The reported normalized covariances correspond to the average of those computed from the sets, with error bars being the standard error of the mean (see Appendix 1 Section 17 for details).

A single mother machine experiment typically produced 100–500 single-cell traces of cells with plasmid and 7–30 single-cell traces of cells that lost the plasmid due to random partitioning at division. The computed normalized covariances from the bootstrapping method are shown in Figure 3D,E, while those from the alternative splitting method are shown in Appendix 1—figure 6. The two methods give similar results.

A complete list of the synthetic circuits can be found in Appendix 1—figures 17–20 . A complete list of strains is in Appendix 1—table 2.

The MG1655 E. coli strain used (SKA360 from Aoki et al., 2019) for the microfluidic experiments with circuits 1, 2, 3, 6, 7, 8, 9, 10, and 11 had gene attB::lacYA177C and the following deletions: ΔaraCBAD, ΔlacIZYA, ΔaraE, ΔaraFGH, ΔrhaSRT, and ΔrhaBADM. The strain used for circuit 5 was identical with an added attn7::PRNA1-mCherry-mKate2-hybrid chromosomal gene. This last strain was also used in Appendix 1—figure 10 for the strain with no plasmid. The MG1655 E. coli strain used (NDL162 from Luro et al., 2020) for the microfluidic experiments with circuits 12, 14, and 15 had chromosomal gene glmS::PRNA1-mCherry-mKate2-hybrid and a motility deletion ΔmotA. The strain used for circuit 13 was identical with the additional ΔrpoS::FRT-Kan-frt. The strain used for circuit 4 was identical to NDL162 but without the chromosomal gene glmS::PRNA1-mCherry-mKate2-hybrid.

The standard heat shock procedure was used to transform plasmid DNA into chemically competent E. coli. During cloning procedures, E. coli strains were grown in LB medium (1% tryptone, 0.5% yeast extract, 1NaCl) with aeration at 250 rpm or LB agar plates at 37°C with appropriate antibiotics at standard concentrations: Ampicillin 100 μg/ml (Amp) and Kanamycin 50 μg/ml (Kan) (antibiotics and LB medium sourced from Fisher Scientific).

Prior to a mother machine experiment, cells were grown overnight in LB medium with appropriate antibiotics from glycerol stocks. At ∼3-4 hr prior to the experiment start time, the overnight cultures were diluted 1:100 in imaging media consisting of M9 salts, 10% (v/v) LB, 0.2% (w/v) glucose, 2 mM MgSO4, 0.1 mM CaCl2, 1.5 μM thiamine hydrochloride and 0.85 g/L Pluronic F-108 (Sigma Aldrich, surfactant to prevent cells sticking to the surface of a microfluidic device). For experiments with circuits 2, 7, 8, and 9 in which IPTG is used, the respective concentrations of IPTG shown in Appendix 1—figures 17–19 are added to the imaging media given to the cells at this point and throughout the remainder of the experiment.

A complete list of plasmids can be found in Appendix 1—table 2. Snapgene files are available on GitHub. Primers and PCR fragments are labeled in the Snapgene files.

Plasmids on the left side of Appendix 1—figures 17–20 were constructed using Gibson assembly (NEBuilder HiFi DNA Assembly Master Mix, NEB, MA, USA). Reactions were performed according to the manufacturers protocols, with Gibson assembly reactions incubated for 15 min at 50°C.

Primers were ordered from Thermo Fisher Scientific and PCR were performed with Accuprime Pfx (Thermo Fisher) or Phusion polymerase (New England BioLabs) according to manufacturers protocols. All plasmids were verified by DNA sequencing, either Sanger sequencing by Quintara Biosciences (MA, USA) or whole plasmid nanopore sequencing with Plasmidsaurus certified by Oxford Nanopore Technologies.

Various fragments from plasmids built in Potvin-Trottier et al., 2016 were amplified by PCR with the ordered primers. Most notably is pLPT119, which formed the backbone for all the plasmid constructs. It consists of the repressilator with no degradation tags on the lacI, cI, and tetR genes (Potvin-Trottier et al., 2016).

To assemble pEJS1 and pEJS10 by phusion reaction, two fragments were amplified from pLPT119 and made up the three repressors from the repressilator, SCFP3A with PLlacO-1 promoter was amplified from pLPT107, and mVenus NB was amplified from pTP85.

To assemble pEJS3 and pEJS4, SCFP3A with PLlacO-1 promoter was amplified from pLPT107, mCherry-mKate2-hybrid with PLtetO-1 promoter was amplified from pLPT27, mVenus fusion from pEJ1, and repressilator fragments were amplified from pLPT119.

To assemble pEJS5, SCFP3A with PLlacO-1 promoter was amplified from pLPT107, mCherry-mKate2-hybrid with PRNA1 promoter was amplified from pLPT26, mVenus fusion from pEJS1, and repressilator fragments were amplified from pLPT119.

For the experiment shown in Figure 4D, we used a pSC101 plasmid with Kanamycin resistance and a gfpmut2 reporter gene with gadX promoter taken from Zaslaver et al., 2006 (ordered from the E. coli Promoter Collection from Horizon Discovery).

The tetR-yfp fusion was made by removing the stop codon from the tetR sequence and replacing it with a short linker sequence followed by the mVenus NB sequence. The first nucleotide from the mVenus NB start codon is not added. The linker sequence is GGGTCTGACATCCTCGAGT.

All repressilator gene sequences and fluorescent protein sequences are followed by transcription terminator T1 from the E. coli rrnB gene.

We let X correspond to the transcript of a gene of interest, with $x$ corresponding to the number of X molecules. Functional interactions in this network are specified by the dependencies of the reaction rates in the network. We say that X affects another component Z if a path can be drawn from X to Z in the topology of the network, see Appendix 1—figure 21A. The set of all components in the network that are affected by X or Y is denoted by ${\mathbf{z}}_{\text{aff}}$. All other components in the cell are denoted as ${\mathbf{z}}_{\text{aff}}^c$, see Appendix 1—figure 21B. We let the X number dynamics be described by a stochastic birth-death process, with production and degradation reactions of X to be given by

where the production rate $R$ can depend in any way on the sets of unspecified components ${\mathbf{z}}_{\text{aff}}$ and ${\mathbf{z}}_{\text{aff}}^c$, and the degradation rate $\beta$ can depend in any way on all components not affected by X or Y, as long as any possible realization of ${\beta }_t$ does not decay to zero over time (i.e., $⟨\beta {⟩}_t=\underset{T\to \mathrm{\infty }}{lim}\frac{1}{T}{\int }_0^T{\beta }_{t}dt>0$). If X is the transcript of geneX, then the rate $R$ corresponds production rate corresponds to the transcription rate of geneX.

We next introduce into the system a passive reporter Y, engineered to read out the same signal as X that degrades with the same rate

where $\alpha$ is an arbitrary proportionality constant. We extend the definition of ${\mathbf{z}}_{\text{aff}}$ to include all components affected by X or Y. These components can in turn depend in arbitrary ways on the number of X and Y molecules. The dynamics of all other cellular components in the cell thus remain unspecified, and the resulting dynamics need not be Markovian or ergodic in X and Y. Note that the sets of unknown components can include an unbounded number of nonphysical mock variables, which can model the different spatial compartments that X can find itself in as well as any complex history-dependent dynamics. We only specify that the total X abundance undergoes first-order degradation, where the degradation rate can depend in an arbitrary way on any of the variables not affected by X or Y.

We can also consider a different class of systems in which X denotes the protein of a gene of interest that is fused to a fluorescent protein. This fluorescent protein must undergo a maturation step before it can fluoresce. We denote the immature fluorescent protein as X, and the transcript for this protein as $M_x$. The mRNA abundances are described by the following stochastic birth-death processes modeling transcription and mRNA degradation

The ‘dark’ fluorescent proteins X_d are translated from the mRNA and then mature into ‘bright’ fluorescent proteins X, which are then degraded:

where $\lambda ({\mathbf{z}}_{\text{aff}}^c)$ is an unspecified translation rate, ${\tau }_{mat}$ is the average maturation time, and ${\beta }_p({\mathbf{z}}_{\text{aff}}^c)$ is an unspecified protein degradation factor. The translation rate $\lambda$ and degradation rate ${\beta }_p$ can fluctuate and can depend in arbitrary ways on all the components not affected by X or $M_x$ or Y or $M_y$. We assume that any possible realization of ${\beta }_p$ does not decay to zero over time (i.e. $⟨{\beta }_p{⟩}_t=\underset{T\to \mathrm{\infty }}{lim}\frac{1}{T}{\int }_0^T{\beta }_{p,t}dt>0$). The fluorescent proteins are assumed to undergo first-order maturation. In Balleza et al., 2018, the authors show that approximately a third of a large set of common fluorescent proteins undergo first-order maturation kinetics in bacteria.

We next introduce into the system a passive reporter Y corresponding to a second fluorescent protein reading out the same transcription rate:

where the transcription and translation rates are proportional to those of geneX with unknown proportionality constants $\alpha$ and ${\alpha }_p$. The average maturation and degradation times are equal to those of X.