### Expressing mRNA tagged with MS2d-GFP fusion protein in *E. coli* DH5α-PRO

The method of RNA detection and quantification was proposed in [22] and characterized in E. coli DH5α-PRO [19]. It exploits the ability of bacteriophage MS2 coat protein to tightly bind specific RNA sequences. High resolution detection of single RNA transcripts with 96 tandem repeats of the MS2 binding sites was demonstrated in *E. coli* by using dimeric MS2d fused to GFPmut3 (MS2d-GFP fusion protein) as a detection tag [18]. The method uses the controlled expression of two genetic constructs: a medium-copy vector that expresses MS2d-GFP fused protein, whose promoter (tetO1) is regulated by tetracycline repressor, and a single copy F-based vector, with a lac/ara promoter controlling the production of the transcript target, mRFP1 followed by a 96 MS2 binding site array. Constructs were generously provided by Ido Golding (University of Illinois).

Cells with both MS2d-GFP and transcript target plasmids were grown in Miller LB medium, supplemented by antibiotics according to the specific plasmids. For full induction of protein and RNA, cells were grown in overnight at 37°C with aeration, diluted into fresh medium to maintain exponential growth until reaching an optical density of OD600 ≈ 0.3-0.5. Inducer aTc (100 ng/ml) was added to get full induction of MS2d-GFP production. Approximately 60 min incubation allows sufficient production for RNA detection. After, expression of target RNA is induced (see below).

Following induction, cells are placed on a microscopic slide between a cover slip and 0.8% LB-agarose gel pad set, and visualized by fluorescence microscopy, using a Nikon Eclipse (TE2000-U, Nikon, Tokyo, Japan) inverted C1 confocal laser-scanning system with a 100× Apo TIRF (1.49 NA, oil) objective. GFP fluorescence is measured using a 488 nm laser (Melles-Griot) and a 515/30 nm detection filter. Images of cells are taken from each slide using C1 with Nikon software EZ-C1, approximately 7 min after induction, one per minute, for approximately 2 hours. Measurements under the microscope were made at room temperature (~ 24°C).

Maximum induction of target RNA is achieved with 1 mM of IPTG and 6.7 mM of arabinose [18]. Besides maximum induction, in one case we induced using 5% of the concentrations needed for maximum induction (weak induction), and in another with 15% (medium induction). At maximum induction we observed approximately 4 RNA/cell/hour, in agreement with previous reports and qPCR measurements [18].

We measured the relative changes in mean mRNA numbers with induction strength with quantitative real time PCR. Target RNA was induced with low and high concentrations of inducers. From isolated RNA, complementary DNA was prepared and used for expression analysis [18]. 16S rRNA was used as an internal control. The Livak method [31] was used to confirm the relative gene expression changes. The following primer pairs were used to amplify the mRFP1 region of the target RNA:

Forward: 5' TAC GAC GCC GAG GTC AAG 3'

Reverse: 5' TTG TGG GAG GTG ATG TCC A 3'

and for 16S rRNA:

Forward: 5'CGT CAG CTC GTG TTG TGA A 3'

Reverse: 5' GGA CCG CTG GCA ACA AAG 3'

For details and results of qPCR measurements see Additional File 2.

### Segmentation of cells, MS2-GFP-RNA spots in cells, RNA molecules from spots, and intervals between transcription events

We detect cells from raw images as in [32]. This method divides a greyscale image in three classes: background, cell border and cell region. It then exploits an iterative cell segmentation process that identifies and segments clumped cells based on size and edge information. To avoid degradation of performance of detection in regions where cells are clumped we apply a threshold based on cell size and discard cells whose size goes beyond the threshold.

The automatic spot detection method segments the MS2d-GFP-RNA spots with the kernel density estimation method for spot detection as in [33]. This method estimates the probability density function over the image from local information, and processes the image by filtering it with a desired kernel. We used a Gaussian kernel and then applied Otsu's thresholding method [34] to segment MS2d-GFP-RNA spots from the kernel density estimated image, highlighting the spots (Figure 2). Finally, the number of RNA molecules in each spot was quantified by normalizing the MS2d-GFP-RNA spot intensity distribution as in [18]. This approach, here named the "slicing approach", consists of estimating the number of tagged transcripts in the cell by dividing a spot's intensity by the intensity of the first peak in the histogram of spot intensities (Additional file 2, supplementary figure 3). An example of a distribution of spot intensities, obtained from the images of multiple cells is shown in Additional File 2.

By performing this analysis for each frame, it is possible to determine when new RNA molecules appear in the cell. From that, we calculate intervals between the productions of consecutive RNAs in individual cells. For a detailed description of this analysis, as well as examples, refer to [18] and Additional files.

Finally, we only count intervals between consecutive RNA molecules which are produced in the same cell. If a cell division occurs, the interval between the last RNA produced in the mother cell and the first RNA produced in a daughter cell is not included in the counts of intervals between consecutive RNA molecules.

### Fitting the model to a d-step model, each step with an exponentially distributed duration

Given the distribution of time intervals between consecutive transcription events, obtained from multiple cells subject to the same induction, it is possible to determine the maximum likelihood fit of a model with *d* statistically independent steps, whose time lengths each follow and exponential distribution. For such a d-step model with parameters μ = [μ_{1}, μ_{2,}..., μ_{d}], and given N measured intervals between transcription events, Δt_{k}, where k goes from 1 to N, the log-likelihood is:

L\left(\mu \right)=\sum _{k=1}^{N}log{\pi}_{d}\left(\Delta {t}_{k};\mu \right)

(3)

where π_{d} is the probability density function for a sum of d exponential random variables with means μ_{d}. The probability density function for the sum can be found by the convolution of the probability density functions of the individual exponential random variables. The density functions for d = 1,..3 are:

{\pi}_{1}\left(\Delta {t}_{k};{\mu}_{1}\right)=\frac{{e}^{\frac{-x}{{\mu}_{1}}}}{{\mu}_{1}}

(4)

{\pi}_{2}\left(\Delta {t}_{k};{\mu}_{1},{\mu}_{2}\right)=\frac{{e}^{\frac{-x}{{\mu}_{1}}}}{{\mu}_{1}-{\mu}_{2}}+\frac{{e}^{\frac{-x}{{\mu}_{2}}}}{{\mu}_{2}-{\mu}_{1}}

(5)

\begin{array}{c}{\pi}_{3}\left(\Delta {t}_{k};{\mu}_{1},{\mu}_{2},{\mu}_{3}\right)=\frac{{\mu}_{1}{e}^{\frac{-x}{{\mu}_{1}}}}{\left({\mu}_{1}-{\mu}_{2}\right)\left({\mu}_{1}-{\mu}_{3}\right)}\\ +\frac{{\mu}_{2}{e}^{\frac{-x}{{\mu}_{2}}}}{\left({\mu}_{2}-{\mu}_{1}\right)\left({\mu}_{2}-{\mu}_{3}\right)}+\frac{{\mu}_{3}{e}^{\frac{-x}{{\mu}_{3}}}}{\left({\mu}_{3}-{\mu}_{1}\right)\left({\mu}_{3}-{\mu}_{2}\right)}\end{array}

(6)

The values of μ = [μ_{1}, μ_{2} ... μ_{d}] are the expected means and standard deviations of the durations of each of the steps composing the intervals between production events. We use this procedure to find the values of μ that provide the highest log-likelihood for d = 1,...,4. No significant improvement of fit was observed for values of d > 2 (Table 1).

We note that the singularities of the probability density functions, formulas (5) and (6), were not problematic since the maximum likelihood estimate of the μ's differed from the second decimal onward. Furthermore, the singularities can be removed. For example, in (5), if μ_{1} = μ_{2}, the singularity can be removed by various means (e.g. L'Hôpital rule), so that π_{2} equals the density function of the gamma distribution with parameters k = 2 and θ = μ_{1} = μ_{2}.

The goodness of fit of the models can be assessed by comparison. For that, we perform a likelihood-ratio test between pairs of models to reject a null model in favour of the alternative. Finally, we verified that the method reliably distinguishes the duration of each step, when they differ by ~25% in duration, from 200 intervals sampled from a model of gene expression.