Precision of readout at the hunchback gene
The simultaneous expression of the hunchback gene in the numerous nuclei of the developing fly embryo gives us a unique opportunity to study how transcription is regulated in living organisms. A recently developed MS2-MCP technique for imaging nascent messenger RNA in living Drosophila embryos allows us to quantify the dynamics of the developmental transcription process. The initial measurement of the morphogens by the hunchback promoter takes place during very short cell cycles, not only giving each nucleus little time for a precise readout, but also resulting in short time traces of transcription. Additionally, the relationship between the measured signal and the promoter state depends on the molecular design of the reporting probe. We develop an analysis approach based on tailor made autocorrelation functions that overcomes the short trace problems and quantifies the dynamics of transcription initiation. Based on live imaging data, we identify signatures of bursty transcription initiation from the hunchback promoter. We show that the precision of the expression of the hunchback gene to measure its position along the anterior-posterior axis is low both at the boundary and in the anterior even at cycle 13, suggesting additional post-transcriptional averaging mechanisms to provide the precision observed in fixed embryos.
During development the different identities of cells are determined by sequentially expressing particular subsets of genes in different parts of the embryo. Proper development relies on the correct spatial-temporal assignment of cell types. In the fly embryo, the initial information about the position along the anterior-posterior (AP) axis is encoded in the exponentially decaying Bicoid gradient. The simultaneous expression of the Bicoid target gene hunchback in the multiple nuclei of the developing fly embryo gives us a unique opportunity to study how transcription is regulated and controlled in a living organism Jaeger (2011); Gregor et al. (2014). Despite many downstream points where possible mistakes can be corrected Driever and Nuesslein-Volhard (1988); Jaeger (2011); Tikhonov et al. (2015), the initial mRNA readout of the maternal Bicoid gradient by the hunchback gene is remarkably accurate and reproducible between embryos Porcher and Dostatni (2010a); Little et al. (2013): it is highly expressed in the anterior part of the embryo, quickly decreasing in the middle and not expressed in the posterior part. This precision is even more surprising given the very short duration of the cell cycles (6-15 minutes) during which the initial Bicoid readout takes place and the intrinsic molecular noise in transcription regulation Elowitz et al. (2002); Ozbudak et al. (2002); Raser and O’Shea (2004).
Even though most of our understanding of transcription regulation in the fly embryo comes from studies of fixed samples, gene expression is a dynamic process. The process involves the assembly of the transcription machinery and depends on the concentrations of the maternal gradients Crauk and Dostatni (2005). Recent studies based on single-cell temporal measurements of a short lived luciferase reporter gene under the control of a number of promoters in mouse fibroblast cell cultures Suter et al. (2011); Zoller et al. (2015) and experiments in E. Coli and yeast populations Taniguchi et al. (2010); Kandhavelu et al. (2012); Muthukrishnan et al. (2012); Chong et al. (2014) have quantitatively confirmed that mRNAs are generally produced in bursts, which result from periods of activation and inactivation. In early fly development, what are the dynamical properties of transcription initiation that allow for the concentration of the Bicoid gradient and other maternal factors to be measured in these short intervals between mitoses?
In order to quantitatively describe the events involved in transcription initiation, we need to have a signature of this process in the form of time dependent traces of RNA production. Recently, live imaging techniques have been developed to simultaneously track the RNA production in all nuclei throughout the developmental period from nuclear cycle 11 to cycle 14 Lucas et al. (2013); Garcia et al. (2013). In these experiments, an MS2-binding cassette is placed directly under the control of an additional copy of a proximal hunchback promoter. As this reporter gene is transcribed, mRNA loops are expressed that bind fluorescent MCP proteins. Their accumulation at the transcribed locus gives an intense localized signal above the background level of unbound MCP proteins (Fig. 1C) Ferraro et al. (2016a). By monitoring the developing embryo, we obtain for each nucleus a time dependent fluorescence trace that is indicative of the dynamics of transcription regulation at the hunchback promoter (Fig. 1B, D and F).
However the fluorescent time traces inevitably provide an indirect observation of the transcription dynamics. The signal is noisy, convoluting both experimental and intrinsic noise with the properties of the probe: the jitter in the signal is not necessary indicative of actual gene switching but could simply result from a momentarily decrease in the recording of the intensity. To obtain a sufficiently strong intensity of the signal to overcome background fluorescence, a long probe with a large number of loops is needed, which introduces a buffering time. In the current experiments the minimal buffering time is the time needed to transcribe a fluorescent probe made of loops. It is s and it prevents direct observation of the activation events Ferraro et al. (2016a).
To understand the details of the regulatory process that controls mRNA expression we need to quantify the statistics of the activation and inactivation times, as has been performed in cell cultures Suter et al. (2011); Zoller et al. (2015); Kandhavelu et al. (2012); Muthukrishnan et al. (2012). However the very short duration of the cell cycles (6-15 minutes for cell cycles 11-13) in early fly development prevents accumulation of statistics about the inactivation events and interpretation of these distributions. Direct observation of the traces suggests that transcription regulation is not static but displays bursts of activity and inactivity. However the eye can often be misleading when interpreting stochastic traces. In this paper we develop a statistical analysis of time dependent gene expression traces based on specially designed autocorrelation functions to investigate the dynamics of transcription regulation. This method overcomes the analysis difficulties resulting from naturally short traces caused by the limited duration of the cell cycles that make it impossible to infer the properties of regulation directly from sampling the activation and inactivation time statistics. Combining our analysis technique with models of transcription initiation, estimates of the precision of the transcriptional readout and high resolution microscopy imaging of the MS2 cassette under the control of the hunchback promoter in heterozygous flies, we find evidence suggesting bursty transcription initiation in cell cycles 12-13. For the switching timescales we observe experimentally, the autocorrelation function analysis alone is not able to reliably distinguish between different models for promoter activation and we use information about the precision of the transcriptional readout to conclude that transcription is most likely bursty. Based on the analysis of the time traces, we show that the precision of the transcriptional readout in each cell cycle is relatively imprecise compared to the expected precision of the mRNA measurement obtained from fixed samples, both in terms of cell-to-cell variability Porcher and Dostatni (2010a) and embryo-to-embryo variability Little et al. (2013). We discuss the limitations of the inference for models of different complexity in different parameter regimes.
Characterizing the time traces
Before we present our results, we first analyze the traces and present a new analysis technique. We study the transcriptional dynamics of the hunchback promoter (depicted in Fig. 1A and B) by generating embryos that express an MS2 reporter cassette under the control of the proximal hunchback promoter (Fig. 1C), using previously developed tools Lucas et al. (2013); Garcia et al. (2013), with an improved MS2 cassette Lucas and et al (2016) (see Materials and Methods for details). The MS2 cassette was placed towards the 3’ end of the transcribed sequence and contained MS2 loop motifs. While the gene is being transcribed, each newly synthesized MS2 loop binds MCP-GFP (expressed at low levels and freely diffusing in the embryo). In each nucleus, where transcription at this reporter is ongoing, we observe a unique bright fluorescent spot, which corresponds to the accumulation of several MCP-GFP molecules at the locus (Fig. 1C). We assume that the fluorescent signal from a labeled mRNA disappears from the recording spot when the RNAP reaches the end of the transgene. With this setup we image the total signal in four fly embryos using confocal microscopy, simultaneously in all nuclei (Fig. 1D) from the beginning of cell cycle (cc) 11 to the end of cell cycle 13. In each nucleus we obtain a signal that corresponds to the temporal dependence of the fluorescence intensity of the transcriptional process, which we refer to as the time trace of each spot. A cartoon representation of such a trace resulting from the polymerase activity (Fig. 1E) dictated by the promoter dynamics (Fig. 1B) is shown in Fig. 1F. We present examples of the traces analyzed in this paper in Fig. 8 and the signal preprocessing steps in the Materials and Methods and Section A.
To characterize the dynamics of the hunchback promoter we need to describe its switching rates between ON states, when the gene is transcribed by the polymerase at an enhanced rate and the OFF states when the gene is effectively silent with only a small basal transcriptional activity (Fig. 1A and B). Estimating the ON and OFF rates directly from the traces is problematic due to the buffering time and to the high background fluorescence levels coming from the unbound MCP-GFP proteins that make it difficult to distinguish real OFF events from noise. To overcome this problem, we consider the autocorrelation function of the signal. To avoid biases from different signal strengths from each nucleus, we first subtract the mean of the fluorescence in each nucleus, and then calculate the steady state connected autocorrelation function of the fluorescence signal (equivalent to a normalized auto-covariance), , at two time points separated by a delay time , and , normalized by the variance of the signal over the traces, according to Eqs. 12 and 13 in Materials and Methods. We limit our analysis to the constant expression part of the interphase (which we call ”steady state” – we discuss this assumption at the end of the Simulated data Results section) by taking a window in the middle of the trace to avoid the initial activation and final deactivation of the gene between the cell cycles (see Materials and Methods and Fig. 9). We will always work with the connected autocorrelation function, which indicates that the mean of the signal is subtracted from the trace. The autocorrelation function is a powerful approach since it averages out all temporally uncorrelated noise, such as camera shot noise or the instantaneous fluctuations of the fluorescent probe concentrations.
Fig. 2A compares the normalized connected autocorrelation functions calculated for the steady state expression in the anterior of the embryo (excluding the initial activation and final deactivation times after and before mitosis) in cell cycles 12 and 13 of varying durations: and minutes. Fig. 2B shows the same functions for traces that have been curtailed to all have equal length. The steady state signal from cell cycle 11 did not have enough time points to gather sufficient statistics to calculate the autocorrelation function. As expected, the functions decay showing a characteristic correlation time, then reach a valley at negative values before increasing again. Since the number of data points separated by large intervals is small the uncertainty increases with . Autocorrelation functions calculated for very long time traces have neither the negative valley nor the increase at large . For example, the long-time connected autocorrelation functions calculated from the simulated traces (Fig. 2C) of the process described in Fig. 1 that are shown in Fig. 2D, differ from the short time connected autocorrelation function in Fig. 2E calculated from the same trace (see Section G for a description of the simulations). As the traces get longer the connected autocorrelation function approaches the longtime results (Fig. 11). The connected autocorrelation function of a finite duration trace of a simple correlated brownian motion (an Ornstein-Uhlenbeck process) displays the same properties (see Fig. 12). The dip is thus an artifact of the finite size of the trace. We also see that the autocorrelation functions shift to the left for short cell cycles (Fig. 2A), resulting, for earlier cell cycles, in shorter directly read-off correlation times, defined as the value of at which the autocorrelation function decays by . However, calculating the autocorrelation functions for time traces of equal lengths for all cell cycles (Fig. 2B) shows that the shift was also a bias of the finite trace lengths, and after taking it into account, the transcription process in all the cell cycles has the same dynamics (although we note that the dynamics directly read out from this truncated trace is not the true long time dynamics).
This preliminary analysis shows that to extract information about the dynamics of transcription initiation we will need to account for the finite time traces. Additionally, a direct readout of even effective rates from the correlation time is difficult, because the autocorrelation coming from the underlying gene regulatory signal (Fig. 1B) is obscured by the autocorrelation due to the timescale needed for the transcription of the sequence containing the MS2 cassette (Fig. 1D) – the gene buffering time, . The observed time traces are a convolution of these inputs (Fig. 1F). The analysis is thus limited by the buffering time of the signal ( in our system), given as the length of the transcribed genomic sequence that carries the fluorescing MS2 loops divided by the polymerase velocity. A direct readout of the switching rates is only possible if the autocorrelation time of the promoter is larger than the buffering time.
The form of the autocorrelation function and our ability to distinguish signal from noise also depends on the precise positioning and length of the fluorescent gene Ferraro et al. (2016a). A construct with the MS2 transgene placed at the 3’ end of the gene (Fig. 3A) gives a differentiable readout of the promoter activity even for two sets of fast switching rates between the active and inactive states. However, in this case the weak signal is hard to distinguish from background fluorescence levels. Conversely, a 5’ positioning of the transgene (Fig. 3B) is insensitive to background fluorescence. However it only differentiates autocorrelation functions calculated from very slow switching processesFerraro et al. (2016a). In summary, a construct with the MS2 placed at the 3’ end of the gene allows for a direct readout of the transcriptional kinetics in a much wider range of switching rates than a 5’ construct, although the autocorrelation function of a 3’ construct is more sensitive to background fluorescence.
Method development – Promoter switching models
The promoter activity we are interested in inferring can in principle be described by models of varying complexity (see Fig. 1A). We consider and compare three types of models in this paper. We note this is a small subset of possible models. In particular, we do not consider models with multiple levels of transcription as was considered in Bothma et al. (2014) or reversible promoter cycles. In the simplest case, the gene is consecutively yet noisily expressed. The RNAP starts transcribing following a Poisson distribution of discrete ON-activation (or firing) events – this has previously been called a static promoter (not represented in Fig. 1A). After the polymerase binds, the next polymerase cannot bind before the promoter is cleared (a timescale estimated to be in our experiments). The effective firing rate of this model is the Poisson rate, , shifted by a deterministic s, , and we call this discrete time model a Poisson-like promoter. Although the promoter dynamics would be uncorrelated in this case, the gene buffering would still produce a finite correlation time (see Section F). Alternatively, the promoter could have two well defined expression states: an ON state during which the polymerase is transcribing at an enhanced level and an OFF state when it transcribes at a basal level. This situation can be modeled by stochastic switching between the two states with rates and (left panel in Fig. 1A and Materials and Methods). However, as was previously observed in both eukaryotic and prokaryotic cell cultures Suter et al. (2011); Zoller et al. (2015); Kandhavelu et al. (2012); Muthukrishnan et al. (2012), once the gene is switched off the system may have to progress through a series of OFF states before the gene can be reactivated. Recently these kinds of cycle models have been discussed for the hunchback promoter Estrada et al. (2016). The intermediate states can correspond to, for example, the assembly of the transcription initiation complex, opening of the chromatin or transcription factor cooperativity. These kinds of situations can either be modeled by a promoter cycle (middle panel in Fig. 1A and Materials and Methods), with a number of consecutive OFF states, or by an effective two state model that accounts for the resulting non-exponential, but gamma function distribution of waiting times in the OFF state (right panel in Fig. 1A and Materials and Methods). The time the polymerase spends transcribing the DNA does not dependent on the promoter model.
In both the two-state and promoter cycle model the gene switches from the ON to the OFF state with exponentially distributed waiting times described by a rate (Fig. 1A). In the two-state model the jumps from the OFF to the ON state are also exponentially distributed with a switching rate (Fig. 1A). In the three state cycle model considered in this paper, an inactive gene can be in two different OFF states. The gene leaves these states with different switching rates, and , respectively. The ordering of and is impossible to detect in the current experiment. In the three state cycle model we can define an effective on-switching rate . corresponds to the inverse of the average waiting time in the overall OFF state, and the waiting times for exiting this effective OFF state are not exponentially distributed. The gamma function distributed switching time is an approximation of this effective rate. We present our method for all of these models and consider all but the gamma function distributed switching time model to learn about the dynamics of hunchback promoter dynamics.
Method development – Autocorrelation approach
To infer the transcription dynamics from the data we built a mathematical model that calculates the autocorrelation functions accounting for the experimental details of the probes, incorporating the MS2 loops at various positions along the gene and correcting for the finite length of the time traces. The basic idea behind our approach is that while the initiation of transcription is stochastic and involves switching between the ON and possibly a number of OFF states ( in Fig. 1B denotes the binary gene expression state), if we assume a constant elongation velocity the obscuring of the signal by the probe design is completely deterministic Coulon et al. (2014); Garcia et al. (2013), which results in the random variable that describes the presence or absence of the polymerase at position at time (Fig.1D). We count the progression of the polymerase in discrete time steps, where one time step corresponds to the time it takes the polymerase to cover a distance of base pairs equal to its own length (Fig. 4A). The promoter dynamics can thus be learned from the noisy autocorrelation function of the fluorescence intensity normalized by the intensity coming from one MS2 loop, (Fig.1F), even for switching timescales smaller than the fluorescent probe buffering time , provided the parameters of the probe design encoded in the loop function (positioning of the probe etc.) are known (Fig. 1C) and the intensity signal is calibrated knowing the fluorescence intensity coming from one MS2 loop Garcia et al. (2013).
Broadly, our model assumes that once the promoter is in an ON state the polymerase binds and deterministically travels along the gene producing MS2 loops containing mRNA that immediately bind MCP and result in a strong localized fluorescence (Fig. 4). The presence or absence of a polymerase at position at time , is simply a delayed readout of the promoter state at time , where is measured in polymerase time steps (Fig. 1B). We assume that polymerase is abundant and that at every time step a new polymerase starts transcribing, provided the gene is in the ON state (Fig. 1B and D). The amount of fluorescence produced by the gene at one time point is determined by the number of polymerases on the gene (Fig. 4A). The amount of fluorescence from one polymerase that is at position on the gene depends on the cumulated number of loops that the polymerase has produced , where . corresponds to the maximum number of polymerases that can transcribe the gene at a given time and corresponds to one loop fluorescing, as depicted in the cartoon in Fig. 4B. The known loop function depends on the build and the position of the MS cassette on the gene, it is input to the model and does not necessarily take an integer value since the polymerase length and the loop length do not coincide (Fig. 4B). Given the average fraction of time the transcription initiation site is occupied by the polymerase, , the average fluorescence in the steady state is:
Since we assume the polymerase moves deterministically along the gene, seeing a fluorescence signal both at time and position , and at time and position means the gene was ON at time and , which is determined by how many loops ( and ) the polymerase has produced. Taking the earlier of these times, we need to calculate the probability that the gene is also ON at the later time. The autocorrelation function of the fluorescence can thus be written as:
where is the probability that the gene is ON at time given that it was ON at time , and time is expressed in polymerase steps. The precise form of , and depends on the type of the promoter switching model. We assume that the polymerase moves at constant speed along the gene and that there is no splicing throughout the transcription process. We give explicit expressions for all the models used in the Materials and Methods section and the Supplementary Sections. Knowing the design of the construct (length of the probe and number of loops that have been transcribed at each position) and calibrating the signal, we use Eq. 1 to directly learn from the data. In the two and multi-state models provides us with the ratio of switching rates and we then use Eq. 2 to obtain their particular values (see Materials and Methods).
To avoid biases coming from nucleus to nucleus variability, we calculated the normalized connected correlation function defined in Eqs. 12 and 13 in Materials and Methods. The theoretically calculated connected autocorrelation function, (Eq. 14, which corresponds to the longtime correlation function in Fig. 2C and D), differs from the empirically calculated connected autocorrelation function from the traces, (Eqs. 12 and 13 in Materials and Methods, which correspond to the short time correlation function in Fig. 2C and E), due to finite size effects coming from spurious correlations between the empirical mean and the data points. Since by definition the mean of a connected autocorrelation function is zero (see Eqs. 12 and 13 in Materials and Methods), the area under the autocorrelation function must be zero. For short traces this produces the artificial dip discussed in Fig. 2, which for long traces is not visible, as it is equally distributed over long times. To compare our theoretical and empirical correlation functions we explicitly calculate the finite size correction and include this correction in our analysis (Materials and Methods and Section H and I).
In this paper, we have analyzed data from fly embryos with 3’ promoter constructs only, limiting ourselves to the steady state part of the trace (see Fig. 9). However the method can also be applied to non-steady state systems (see Section C) and to other constructs, including cross-correlation functions calculated from signals of different colors inserted at different positions along the gene (see Section J). We use simulated data to show that prediction and inference are possible for cross-correlation functions of a two-colored signal (see Fig. 13), but that the accuracy of inference is limited by the use of the 5’ probe.
To check that the inference method correctly infers the parameters of the model, we first tested the autocorrelation based inference on simulated short-trace data with underlying molecular models with different levels of complexity ( Fig. 3C) for a construct with the MS2 probe placed at the 3’ end of the gene (Fig. 3B). In Fig. 3D we compare autocorrelation functions for the three state model for constructs with the MS2 loops positioned at the beginning of the transcribed region (5’, Fig. 3B) and at the end of the transcribed region (3’, Fig. 3A), and the cross-correlation function calculated from a two-colored probe construct (Fig. 3E). The analytical model correctly calculates the short trace autocorrelation function and is able to infer the dynamics of promoter switching for all models. It can also be adapted to infer the promoter switching parameters for any intermediate MS2 construct position, given the limitations of each of the constructs discussed above Ferraro et al. (2016a).
The autocorrelation function based inference reproduces the underlying parameters of the dynamics with great accuracy not just for switching timescales longer than the gene buffering time, , that obscures the signal (Fig. 3F), but also for smaller timescales that are within an order of magnitude of the gene buffering time. In Fig. 3F we show the results of the inference for the 3’ two state model for different values of the ON and OFF rates, and . For switching timescales much shorter than the gene buffering time, the autocorrelation function coming from the length of the construct dominates the signal and the precision of the inference goes down. For very fast switching rates (), increasing the length of the traces or the number of nuclei (red vs blue curve for values of larger than in Fig. 3F) does not help estimate the properties of transcription. In this regime, the inferred value of disagrees with the true parameters even when the inference uses long time traces and a large number of nuclei. For intermediate switching rates (), increasing the trace length or increasing the number of nuclei extends the inference range (black and green dashed lines vs blue solid line in Fig. 3F), and in all cases increasing the number of nuclei decreases the uncertainty as can be seen from the smaller error bars (shown only for the red and blue lines for figure clarity).
Using two colored probes attached at different positions along the gene gives two measurements of transcription and allows for an independent measurement of the speed of the polymerase – one of the parameters of the model that currently must be taken from other experiments. While the estimates of polymerase speed in the fly embryo are reliable Garcia et al. (2013), this parameter has been pointed out as a confounding factor in other correlation analyses Coulon and Larson (2016).
The autocorrelation approach also correctly infers the parameters of transcriptional processes when applied to traces that are out of steady state (see Section C). However, since the process is no longer translationally invariant more traces are needed to accumulate sufficient statistics. For this reason, in the current analysis of fly embryos we do not analyze the transient dynamics at the beginning and end of each cycle and we restrict ourselves to the middle of the interphase assuming steady state is reached (see Fig. 9 for details). We do not know whether the underlying dynamics is completely in steady state. We limit our analysis to a time frame window where the intensity of the fluorescence signal plateaus (see Fig. 9 for an example). We can motivate the steady state assumption a posteriori: the inferred switching timescales (smaller than s) are small enough for the system to relax to steady state within one cell-cycle. However we cannot fully rule out other mechanisms that could keep the system out of steady state (such as changes in the Bicoid concentration).
Fly trace data analysis
We divided the embryo into the anterior region, defined as the region between and of the egg length (the position at of the egg length marks the embryo midpoint), where hunchback expression is high, and the boundary region, defined as the region between and egg length, where hunchback expression decreases. The mean probability for the gene to be ON during a given cell cycle (restricted to the times excluding the initial activation and deactivation of the gene, which we will call the steady state regime), given by Eq. 1, is consistent between the four embryos in cell cycle 12 and 13, both in the anterior region and at the boundary (Fig. 5A). The probability for the gene to be ON is over three fold higher in the anterior region than in the boundary and does not change with the cell cycle. in the anterior indicates that in each nucleus the polymerase spends about half the steady state expression time transcribing the observed gene. At the boundary the gene is transcribed on average during about of the steady state part of the cell cycle. The estimates for in the earlier cell cycles were not reproducible between the four embryos, likely because the time traces were too short to gather sufficient statistics to accurately calculate the maximum and average of the signal. We concentrated on cell cycle 12 and 13 for the remainder of the analysis.
Initial comparison of the promoter models
Based on the different behavior at the boundary and in the anterior, we separately inferred the transcriptional dynamics parameters in the two regimes, using the autocorrelation approach that corrects for finite time traces for different models. The Poisson-like promoter model, the two and three state cycle models all provide reasonably good fits to all the traces in both regions (see Fig. 5B for an example and Fig. 10 for the fits in both regions in all embryos). Although the fit of the Poisson-like promoter model (red line) only captures the short time behavior of the measured autocorrelation function, there is not enough statistical evidence from the autocorrelation analysis to exclude this model. The two and three state model fits are indistinguishable (Fig. 5B) and the two state fit is reproducible between cell cycles and embryos (Fig. 5C and D). The form of the autocorrelation function in the Poisson-like promoter model is completely determined by the autocorrelation signal of the fluorescent construct (the loop function ), since the random firing is itself uncorrelated. In the two and three state models, the autocorrelation signal from the fluorescent construct is convoluted with the autocorrelation signal of the promoter. The fact that the Poisson-like promoter model fits the data so well, indicates that the autocorrelation time of the promoter is comparable to the autocorrelation time of the fluorescent construct. In Fig. 10 we plot the autocorrelation functions for simulated two state models with different correlation times () and a Poisson-like promoter with the same . For short correlation times, the Poisson-like promoter model and two state model have indistinguishable autocorrelation functions, just like in the analyzed data. For long autocorrelation times, the difference between the two models is clear.
The three state fit is reproducible at the level of the sum of the effective ON and OFF rates (same fit as shown for the two state model in Fig. 5C), but gives fluctuating values for , the parameter ratio determining how well it is approximated by a two state model (see Fig. 15, describes one fast reaction between the OFF states, effectively giving a two state model, while gives equal weights to the two reactions, clearly distinguishing two OFF states). Since the two state model is reproducible, and has lesser complexity we will further consider only the two state and Poisson-like promoter models.
Discussion of the two state model
For the two state model, the inference procedure independently fits the characteristic timescale of the process from the autocorrelation function, defined as the inverse of the sum of two rates, (Fig. 5C) , and then uses an independent fit of the probability of the gene to be ON, , to disentangle the two rates (Fig. 5D). Examples of the simulated promoter state over time with the rates’ inferred values are shown in Fig. 5E (for the anterior region) and Fig. 5F (for the boundary region). Assuming the two state model we find that the characteristic timescale, , in most embryos is slightly shorter at the boundary () than in the anterior region () and the variability between the two cell cycles is comparable to the embryo-to-embryo variability (Fig. 5C). Both timescales are much larger than the polymerase blocking time, , during which a second polymerase cannot bind because the first one has not cleared the binding site (shown as the gray dashed line in Fig. 5D), which sets a natural scale for the timescales we can infer. We find that in the anterior region of the embryo the two switching rates, and , show variability from embryo to embryo (between to – see Table I and II), but always scale together, which gives the observed one-half probability of the gene to be ON in a given nuclei during the steady state part of the interphase. Since the polymerase in the anterior on average spends half the steady state interphase window transcribing the gene, the inferred rates suggests a clear bursting behavior of the transcription process, with switching between an identifiable active and inactive state of the promoter if the two state promoter is correct, and rare firing if the Poisson-like promoter model is correct.
At the boundary is much smaller than in the anterior with very little embryo to embryo variability, while has a similar range in the anterior and at the boundary. This behavior is expected since high Bicoid concentrations in the anterior upregulate the transgene whereas lower concentrations at the boundary result in smaller activation rates. The ratio of the average rates at the boundary and anterior is , which can be compared to the fold decrease expected from pure Bicoid activation, assuming the Bicoid gradient decays with a length scale of Gregor et al. (2007a) and comparing the activation probabilities in the middle of the anterior and boundary regions. Given the crudeness of the argument stemming from the variability of the Bicoid gradient in the boundary region and the uncertainty of the inferred rates, these ratios are in good agreement and suggest that a big part of the difference in the transcriptional process between the anterior and boundary is due to the change in Bicoid concentration. Of course other factors, such as maternal Hunchback, could also affect the promoter, leading to discrepancies between the two estimates.
Discriminating between two and three state models
The current data coming from four embryos and nuclei in each region with trace lengths of does not make it possible to distinguish between the two and three state models. We asked whether having longer traces or more nuclei could help us better characterize the bursty properties. We performed simulations with characteristic times similar to those inferred from the data () assuming a two state (, Fig. 6A) and three state cycle model (Fig. 6B). We then inferred the sum of the ON and OFF rates () and the ratio of the two OFF rates (). If the two OFF rates are similar (), we infer a three state model. If one of the rates is much faster (, which implies ), we infer a two state model. We find that having more nuclei, which corresponds to collecting more embryos, would not significantly help our inference. However looking at longer traces would allow us to disambiguate the two scenarios, if the traces were times longer, or minutes long. Since cell cycle 14 lasts minutes, analyzing these traces could inform us about the effective structure of the OFF states. However in cell cycle 14, several other direct Bicoid targets (mostly transcription factors) are likely to be expressed, so additional regulatory elements could be responsible for the observed transcriptional dynamics compared to cell cycle 12 and 13. Our results suggest that with our current trace length we should be able to identify a two state model with large certainty, but we could not clearly identify a three state model. Our data may point towards a more complex model than two state, but different kinds of multistate models or a two state model obscured by other biases cannot be ruled out.
The error bars for the autocorrelation functions (see Fig. 5B and Fig. 3) describe the variability between nuclei coming from both natural variability and measurement imprecision. While the autocorrelation function is insensitive to white noise, it does depend on correlated noise. The noise increases for large time differences , as the number of pairs of time points that can be used to calculate the autocorrelation function decreases and in our inference we reweigh the points according to their sampling, so that the noise does not impair the precision of our inference. The error bars on the inferred parameters (e.g. Fig. 5A, C and D) are due to variability between nuclei and are obtained from sampling different subsets of the data in each region and cell cycle. Additionally to the inter-nuclei and experimental noise, there is natural variability between embryos. Since each nucleus transcribes independently and we assume similar Bicoid concentrations in each of the regions, the inter-embryo variability is of a similar scale as the inter-nuclei variability (Fig. 5C), as one expects given that the Bicoid gradient is reproducible between embryos Gregor et al. (2007a). Additonally, variability between Bicoid gradients in different embryos, for example due to their different lengths Cheung et al. (2014), could also contribute to the observed variability.
Accuracy of the transcriptional process
At the boundary, neighboring nuclei have dramatically different expression levels of the Hunchback protein. From measurements of the Bicoid gradient, Gregor and collaborators estimated that for two neighboring nuclei to make different readouts, they must be able to distinguish Bicoid concentrations that differ by Gregor et al. (2007b). Following the Berg and Purcell Berg and Purcell (1977) argument for receptor accuracy, and using measurements of diffusion constants for Bicoid proteins from cell cycle 14, the authors showed that, based on protein concentrations, the hunchback gene is not able to read-out the differences in the concentrations of Bicoid proteins to the required accuracy in the time that cell cycle 14 lasts. Even considering the revised higher values of Bicoid’s diffusion coefficient measured in a subsequent study Porcher and Dostatni (2010a), the precision of the Bicoid gradient read-out remains difficult to explain. The authors invoked spatial averaging of Hunchback proteins as a possible mechanism that achieves this precision. Spatial averaging can increase precision, but it can also smear the boundary. Erdmann et al calculated the optimal diffusion constant Hunchback proteins must have for the averaging argument to work Erdmann et al. (2009) and showed it is consistent with experimental observations Porcher and Dostatni (2010a); Gregor et al. (2007a). However precision can already be established at the mRNA level and, using measurements on fixed embryos, Little and co-workers found that the relative intrinsic nuclei-to-nuclei variability of the mRNA transcribed from a hunchback locus is Little et al. (2013). Measurements of cytoplasmic mRNA reduced this variability to Little et al. (2013).
Here we go one step further and use our direct measurements of transcription from the hunchback gene to directly estimate the precision with which the hunchback promoter makes a readout of its regulatory environment in a given cell cycle in a given region of the embryo, . is the relative error of the probability of the gene to be ON averaged over the steady state part of a cell cycle. Since the total number of mRNA molecules produced in a given cycle is proportional to (shown in Fig. 16A as a function of embryo length), the precision at the level of produced mRNA in a given cycle is equal to the precision in the expression of the gene, . The accuracy of transcription activation is encoded in the stochasticity of gene activation.
In the two state model, the gene randomly switches between two states: active and inactive, making a measurement about the regulatory factors in its environment and indirectly inferring the position of its nucleus. Since no additional information is provided by a measurement that is strongly correlated to the previous one, the cell can only base its positional readout on a series of independent measurements. Two measurements are statistically independent, if they are separated by at least the expectation value of the time it takes the system to reset itself:
where in a two state model and . A more detailed estimate obtained by computing the variance of the time spent ON by the gene during the interphase (see Section K) shows that Eq. 3 underestimates the time needed to perform independent measurements. We find that for a two state model the accuracy of the readout of the total mRNA produced is limited by the variability of a two state variable divided by the estimated number of independent measurements within one cell cycle:
where is the duration of the cell cycle and the factor is a prefactor correction to the naive estimate. Eq. 4 is valid in the limit of (the exact result if given in Section K). Using the rates inferred from the autocorrelation analysis (Fig. 5D) we see that the precision of the gene readout is much lower at the boundary than in the anterior, does not change with the cell cycle and is reproducible between embryos (blue and red points on the ordinate in Fig. 7A). In the anterior part of the embryo it reaches , while at the boundary, it is very large, , even at cell cycle 13.
In the Poisson-like promoter model, to calculate the relative error in the total mRNA produced, the polymerase arrival times are described by an effective firing rate of . Within this model, the fraction of the total time the polymerase cannot bind, because the binding site is occupied is (see Section F). The total produced mRNA is then proportional to the time the gene is transcribed, and the relative error in the total produced mRNA depends on the relative error of the firing times of this modified Poisson process and the number of independent measurements, (see Section K):
Using the rates for inferred from the data (Fig. 5A), the relative error in the total mRNA produced (green and purple points on the ordinate in Fig. 7A) is slightly higher in the boundary region () than in the anterior of the embryo () and does not change with the cell cycle.
We can compare both of these theoretical estimates with direct estimates of the relative error of the total mRNA produced during a cell cycle, , from the data. We divide the embryo into anterior and boundary strips, as we did for the inference procedure and calculate the mean and variance of . These empirical estimates of the gene measurement precision agree with the theoretical estimates (Fig. 7A) for the two state model, but disagree with the predictions of the Poisson-like promoter model. For completeness, we also calculated the relative error for the three state model (see Fig. 16B), which shows better agreement than the Poisson-like promoter model but slightly worse than a two state model. We verified that our conclusions about the scale of our empirical estimates are the same for all embryos (Fig. 16C) and do not depend on the definition of the boundary and anterior regions (Fig. 16D). Since the predicted relative error for the Poisson-like promoter model is lower than the relative error calculated directly from the data, we could imagine that the data estimate is more susceptible to additional sources of experimental noise. However the very large disagreement between the Poisson-like promoter prediction and the data at the boundary suggests the Poisson-like promoter model is not an accurate description. This higher experimental variability also cannot be explained by variable expression levels within the regions we are considering: Fig. 16D shows that the experimental relative error does not significantly decrease if we take smaller windows and taking the in the boundary region to range from to (Fig. 3A) would translate into a, at best, two fold increase in the relative error predicted from the Poisson-like promoter model, which is not enough to reach the experimentally observed relative error. While we are unable to rule out the Poisson-like promoter model based on the fit to the autocorrelation function, a different statistic - the relative error in the produced mRNA - suggests that the promoter is most likely well described by a two state model, and possibly a three state cycle.
To see whether temporal integration of the mRNA produced can increase precision, we compared the empirical estimate of the steady state mRNA production (red line in Fig. 7B) to the relative error of the total mRNA produced in cell cycle 13 (blue line in Fig. 7B) and the total mRNA produced from cell cycle 10 to 13 (green line in Fig. 7B) averaged over embryos. Assuming that the mRNA molecules are equally divided between daughter cells during division, and they are all kept in the cell throughout cell-cycles (which is incorrect but provides a best case estimate), then each nuclei has the total mRNA produced in cell cycle 13, of the total mRNA produced by its mother in cell cycle 12, of the mRNA produced by its grandmother in cell cycle etc. While we see about a increase in the precision at the boundary from integrating the mRNA produced in different cell cycles, the estimate in the anterior region is not helped by integration over the cell cycles.
For completeness of the discussion of the relative errors in the different models, we calculated the relative error assuming the same for a three state cycle () as for a two state model () for different values of and (Fig. 7C). We found that the relative error is always lower for the three state cycle model and the error decreases regardless of the duration of the cell cycle. As expected from Eq. 4, the relative error is decreased by increasing and decreasing . However the increase in precision from a three state cycle model in the parameter regime we inferred for the two state model in the the fly embryo is relatively modest (from for the two state model to for the three state model). Similarly, in Fig. 7D we compared the prediction for the relative errors for the Poisson-like promoter model to two state models with the same probability of the gene to be transcribed, , but different switching rates between the two states ( and ). Faster switching increases the precision of the two state promoter, since the number of independent measurements increases. The Poisson-like promoter is always more accurate than the two state promoter.
Many previous analysis of precision from static images calculated the relative error of the distribution of a binary variable, which in each nucleus was if the nucleus expressed mRNA in the snapshot, and if it did not expressPorcher and Dostatni (2010b); Perry et al. (2012). We analyzed our data using this definition of activity (see Fig. 16E for mean activity as a function of position) and found that for most embryos the relative error in the anterior drops to zero (Fig. 16F), indicating that all nuclei in a given region show the same expression state, but at the boundary the precision is still , in agreement with previous reports about the total mRNA in the nucleus Little et al. (2013). This provides additional evidence for the bursty nature of transcription in the anterior of the embryo, in agreement with previous results that showed a relationship between Bicoid concentration and transcriptional burst of downstream genes He et al. (2011).
In contrast to previous studies Lucas et al. (2013); Garcia et al. (2013), including ours, which failed to show evidence for bursty switching of the hunchback promoter, by developing more advanced analysis techniques we show that the promoter has distinct periods of enhanced polymerase transcription followed by identifiable periods of basal polymerase activity. Our conclusions are based on combining a new autocorrelation based analysis approach, applied to live imaging MS2 data to infer switching parameters, with an analysis of the precision of readout of promoter. The data we used in this paper was generated with a modified MS2 cassette Lucas and et al (2016) (see the Experimental procedures section in Materials and Methods) compared to the previously published data Lucas et al. (2013). However the difference in our conclusions mainly comes from a detailed analysis of the traces.
Quantification of transcription from time dependent fluorescent traces in prokaryotes and mammalian cell cultures has shown that the promoter states cycle through at least three states Suter et al. (2011); Zoller et al. (2015). In one of these states the polymerase transcribes at enhanced levels, while in most of the remaining states the transcription machinery gets reassembled or the chromatin remodels. We find that in the living developing fly embryo, the hunchback promoter also cycles through at least two states, although based on the parameter inference alone we cannot conclusively rule out the possibility of a static promoter with a Poisson-like firing rate or of a more complex promoter with more effective states when the gene is inactive. Only a combination of the inferred parameters using the autocorrelation function with another statistic (the relative error in the produced mRNA) allows us to favor the two state model (or more complex models) over the other considered mode of transcription activation.
The main impediment to distinguishing different types of transcriptional models comes from the very short durations of the interphase in the early cell cycles when the hunchback gene is expressed. We showed using simulations that increasing the number of embryonic samples would not help us distinguish between two and three state models, however looking at longer time traces would be informative (Fig. 6). Since cell cycle 14 lasts about 45 minutes, our analysis shows that the steady state part of the interphase provides enough time to gather statistics that can inform us about the detailed nature of the bursts. Unfortunately, other transcription factors, such as the other gap genes regulating hunchback expression in cell cycle 14, could possibly change the nature of the transcriptional dynamics in a time dependent manner. We showed that the transcriptional dynamics is constant and reproducible in the earlier cell cycles (12-13) (Fig. 2), so independently of the question of the nature of the bursts, it would be very interesting to see whether and how it changes when the nature of regulation changes.
In the parameter regime of relatively fast switching that we inferred from our data, the autocorrelation function for the two state model and the Poisson-like promoter model are very similar. In this parameter regime, the form of the autocorrelation function is governed by the autocorrelation of the fluorescent probe. So while the autocorrelation function approach is able to disentangle the real promoter switching from the buffering of the construct to determine the parameters assuming an underlying model, we cannot conclusively discriminate between these two models, without looking at other statistics. Using simulations (Fig. 14), we showed that for promoters with slower switching characteristics, this discrimination task is possible and the autocorrelation function approach alone can reliably discriminate between different models. In the parameter regime inferred for the hunchback promoter, having longer traces would not be helpful for this discrimination task and we have to look for other statistics (Fig. 17). However using new constructs with MS2 binding sites that have higher binding affinity to MCP and decrease the noise from the binding/unbinding of MCP to the RNA would make it possible to use shorter MS2 cassettes without increasing background fluorescence. These cassettes would decrease the buffering time and extend the parameter regime in which we can distinguish between the Poisson-like and two state promoter models.
Alternatively to focussing on longer traces, a construct with two sets of MS2 loops placed at the two ends of the gene that bind different colored probes could be used to learn more about transcription dynamics Fukaya et al. (2016). We do not have access to data coming from such a promoter , but our analysis approach can be extended to calculate the cross-correlation function between the intensities of the two colored probes. Such cross-correlation analysis has previously been used to study transcription in cell cultures Wang et al. (2011), transcriptional noise Lin et al. (2015) and regulation in bacteria Dunlop et al. (2008); Munsky et al. (2012) . Our theoretical prediction for such a cross-correlation function agrees with simulation results (Fig. 3C). Unfortunately, the cross-correlation function with one set of probes inserted at the 5’ end and the other at the 3’ shares the same problems of a 5’ construct. For fast switching rates, such a cross-correlation function suffers from the large buffering time ( inGarcia et al. (2013)) drawback of the 5’ design and can only be used for inferring large switching rates Ferraro et al. (2016b) (see Fig. 13). Similarly, the cross-correlation function cannot discriminate between a two state and Poisson-like promoter for relatively fast switching. However, it does gives us access into dynamical parameters of transcription such as the speed of polymerase and it is able to characterize whether mRNA transcription is in fact deterministic and identify potential introns. Possibly, cross-correlations from two colored probes both inserted closer to the 3’ end could be optimal designs.
Our method requires knowing the design of the experimental system (number and position of the loops), the speed of polymerase as input and calibrating the maximal fluorescence from one gene. While the polymerase speed is an important parameter and erroneous assumption could influence the inference, we have shown that our inference is relatively insensitive to polymerase speeds (see Fig. 18). In the current experiments we do not have an independent calibration of the maximal fluorescence coming from one gene, which could introduce potential errors in our analysis. However the reproducibility of our results suggests that these potential errors are small.
We assumed an effective model that describes the transcription state of the whole gene and does not explicitly take into account the individual binding sites. As a result all the parameters we learn are effective and describe the overall change in the expression state of the gene and not the binding and unbinding of Bicoid to the individual binding sites. For concreteness we presented our model assuming a change in the promoter state and constitutive polymerase binding, but our current model does not discriminate between situations where the transcriptional kinetics are driven by polymerase binding and unbinding and promoter kinetics. The presented formalism can be extended to more complex scenarios that describe the kinetics of the individual binding sites and random polymerase arrival times. Since we already have little resolution power to discriminate between these effective models, we chose to interpret the results of only these effective models. The exact contribution of the individual transcription binding sites could be inferred from the activity of promoters with mutated binding sites. Similarly, other more complex models, such as a reversible three state model, or a model with many ON states, have not been ruled out by our current analysis but are possible within the current framework.
The time traces we had to analyze are very short and finite size effects are pronounced. Unlike in cell culture studies, where long time traces are available, we could not collect enough ON and OFF time statistics to characterize the promoter dynamics from the waiting time distributions. In this paper we show that simple statistics, the auto- and cross-correlation functions are powerful general tools that can be used in these kinds of challenging circumstances. To reach our final conclusion we had to combine different kinds of statistics, which is also a useful strategy when limited by data.
The approach we propose is a general method that can be used for any type of time trace analysis. However it becomes very useful when studying in vivo biological processes, where the biology naturally limits the available statistics. In our case the number of ON and OFF events is naturally limited by the short duration of the cell cycles. Our method explicitly calculates correlation functions for short traces, correcting for the finite size effects, and can be also used without making steady state assumptions about the dynamics (although this requires collecting sufficient statistics about two time points, which may be hard for short traces). With these corrections we see that while an effective two state model of the underlying dynamics of transcription regulation holds in the anterior and boundary regions of the embryo in all of the early cell cycles, the rates are different in the boundary and anterior regions, showing a strong dependence on position dependent factors such as Bicoid or maternal and zygotic Hunchback concentrations. More statistics will make it possible to build more explicit models of Bicoid dependent activation.
While our method is able to deconvolute the effects of the fluorescent probe and infer rates below the buffering limit of the probe (in our case s, see Fig. 3F), in all cases, the rates that we can infer from time dependent traces are naturally limited by the timescales at which the polymerase leaves the promoter, which in our case is estimated to be . If the switching rates are faster than this scale, even a perfect, noiseless and infinitely accurate sampling of the dynamics will not be able to overcome this natural limit.
The inferred rates are reproducible between nuclei and embryos and the inter-embryo variability is similar to the inner-embryo variability (Fig. 5A, C and D). The embryo-to-embryo variability can come from Bicoid variability, which is Gregor et al. (2007b), so we do not expect the observed expression variability to be less, variability in growth rate and RNAP availability and external environmental factors. Additional sources of noise are experimental noise and most importantly problems with data calibration of what is the maximal level of fluorescence intensity.
We used the obtained results to estimate the precision of the transcriptional process from the hunchback promoter. We found that even in the anterior region, the variability in the mRNA produced in steady state by the different nuclei is large, with a relative error of about (Fig. 7A). This variability further increases to of the mean mRNA produced at the boundary. These empirical estimates are completely explained for a two state promoter model by theoretical arguments, which treat the gene as an independent measuring device that samples the environment, correcting for the number of independent measurements during a cell cycle. In both cases, the precision at the level of the gene readout is not sufficient to form the precise Hunchback boundary up to half a nuclear width Gregor et al. (2016). Even extending our argument to the total mRNA produced in the early cell cycles (Fig. 7B) does not help. Having an irreversible promoter cycle could increase the theoretical precision, but only slightly in the parameter regime we have inferred and it would not change the quantitative conclusions about low precision backed by the empirical results. A Poisson-like promoter, while not compatible with the observed error rates, does have a significantly smaller error.
The construct we used here was limited to the bp of the proximal hunchback promoter, which is known to recapitulate the hunchback endogenous expression observed in Fluorescent In Situ Hybridization (FISH) Lucas and et al (2016). It is possible that the boundary phenotype is recovered by averaging of mRNAs and proteins produced by the real gene or the transgenes in other nuclei. In the latter case, this would point towards a robust ”safety” averaging mechanism that relies on the population. Alternatively, we have to be aware that the sharp boundaries were only detected on fixed samples and that having access to the dynamics of the transcription process likely provides a more accurate view on the process. We calculated and estimated from the data the precision of the gene readout based on the variability of the transcription process between nuclei. We find that the transcriptional process at a given position is quite noisy. Previous estimates of precision were based on data from fixed samples and did not consider the probability of the gene to be ON, but assumed a binary representation where each nuclei is either active or inactive. By analyzing the full dynamic process we show that the gene is bursty and the transcriptional process itself is much more variable. Reducing the information contained in our traces to binary states, we find precise expression in the anterior, but still large variability at the boundary, similarly to previous results from Fluorescent In Situ Hybridization (FISH) aiming to detect all mRNAs Little et al. (2013).
Assuming that the precision in determining the position of the nuclei is encoded in the precision of the gene readout, a gene with the dynamics characterized in this paper needs to measure the signal times longer at the boundary to achieve the observed precision. A gene in the anterior would need to integrate only times longer. These results again suggest that the precision in determining the position of the nuclei is not only encoded in the time averaged gene readout, but probably relies either on spatial averaging mechanisms Gregor et al. (2007b); Okabe-Oho et al. (2009); Erdmann et al. (2009) or more detailed features of the temporal information encoded in the full trace He et al. (2011).
In summary, the early developing fly embryo provides a natural system where we can investigate a functional setting the dynamics of transcription in a living organism. In our data analysis we are confronted with the same limitations that natural genes face: an estimate of the environmental conditions must be made in a very short time. Analysis of dynamical traces suggests that transcription is a bursty process with relatively large inter-nuclei variability, suggesting that simply the templated one to one time-averaged readout of the Bicoid gradient is unlikely. Comparing mutant experiments can shed light on exactly how the decision to form the sharp hunchback mRNA and protein boundary is made.
Materials and Methods
For live monitoring of hb transcription activity in Drosophila embryos, we used the MS2-MCP system, which allows fluorescent labeling of RNAs as they are being transcribed Bertrand et al. (1998); Lucas et al. (2013); Ferraro et al. (2016b). To implement the reporter system in embryos, we generated flies transgenic for single insertions of a P-element carrying the hb proximal promoter upstream of an iRFP-MS2 cassette carrying 24 MS2 repeats Janicki et al. (2004); Lucas et al. (2013), from which Zelda binding motifs have been removed Lucas and et al (2016). The flies also carry the PmRFP-Nup107.K Katsani et al. (2008) transgenic insertion on the 2 chromosome and the Pw[+mC]=Hsp83-MCP-GFP transgenic insertion on the 3rd chromosome. These allow the expression of the Nucleoporin-mRFP (mRFP-Nup) for the labeling of the nuclear envelopes and the MCP-GFP required for labeling of nascent RNAs Bertrand et al. (1998). All stocks were maintained at 25C.
Embryo collection, dechorionation and imaging have been done as described in Lucas et al. (2013). Image stacks (19Z 0.5m, 2m pinhole) were collected continuously at 0.197m XY resolution, 8bits per pixels, 1200x1200 pixels per frame. A total of 4 movies capturing 4 embryos from nuclear cycle 10 to nuclear cycle 13 were taken. Each movie had different scanned fields along the embryo’s width, which results in different time resolutions: 13.1 s, 10.2 s, 5.1 s and 4.3 s.
Nuclei segmentation, tracking and MS2-MCP loci analysis were performed as in Lucas et al. (2013) and recapitulated here. All steps were inspected visually and manually corrected when necessary. Nuclei segmentation and tracking were done by analyzing, frame by frame, the maximal Z–projection of the movies’ mRFP-Nup channel. Each image was fit with a set of nuclei templates, disks of adjustable radius and brightness comparable with those of raw nuclei, from which the nuclei’s positions were extracted. During the cycle’s interphase, each nucleus was tracked over time with a simple minimal distance criterion. For MS2-MCP loci detection and fluorescent intensity quantification, the 3D GFP channel (MS2-MCP) were masked with the segmented nuclear images obtained in the previous step. This procedure also helps associating spots to nuclei. We then applied a threshold equal to times the background signal to the masked images and selected only the connected regions with an area larger than 10 pixels. The spot positions were set as the position of the centroids of the connected regions. The intensity of each spot was calculated by summing up the total pixel intensity in the vicinity of the centroids (a region of 1.5m x 1.5m x 1m) subtracted from the background intensity, which was extracted from the nearby region but excluding the spots. The flies were heterozygous for the MS2 reporter, so one spot was visible at a time. In the (rare) case of multiple spots detected per nucleus, the biggest spot was selected.
For each nucleus, we collected the nucleus’ position and the spot intensity over time (here referred to as ”traces”). The traces were then classified according to their respective embryos (out of 4 embryos), cell cycle (10 to 13) and position along the AP axis (either anterior or boundary). See Section A and Fig. 8 for examples of traces.
Before the autocorrelation function can be calculated the traces need to be preprocessed. To ensure that the data captures the dynamics of gene expression in its steady state, for each embryos and each cell cycle, we observed the spot intensity only in a specific time window. The beginning and the end of this window is determined as the moment the mean spot intensity over time of all traces (both at the anterior and the boundary) reaches and leaves an expression plateau (see example in Fig. 9).
The two state model
The detailed form of the autocorrelation function in Eq. 2 depends on the underlying gene promoter switching model. For the two state – telegraph switching model (left panel in Fig. 1A) the jumping times between the two states are both exponential and the dynamics is Markovian. The mean steady state probability for the promoter to be ON is , which combined with Eq. 1 gives the form of the mean fluorescence . The probability that the gene is ON at time given that it was on at time is , where . The steady state connected correlation function depends only on the time difference (see Section B):
The Cycle model
In the cycle model (also called the three state model in the text) (center panel in Fig. 1A) the OFF period is divided into different sub-steps that correspond to intermediate states with exponentially distributed jumping times from one to the next. The transition matrix encodes the rates of this irreversible chain. The probability of the promoter to be in the ON state is:
and that the steady state connected autocorrelation function is (see Section D):
where is counted in polymerase steps , Id is the identity matrix of dimension , and the two unit vectors are of dimension . In the simple case of a two state model Eq. The Cycle model reduces to Eq. 6.
The waiting time model
An alternative description of a promoter cycle relies on a reduced description to an effective two state model where we use the fact that the transitions between the states are irreversible. The distribution of times spent in the effective OFF state , is no longer exponential, as it was in the two state model, but it has a peak at nonzero waiting times, which can be approximated by a Gamma distribution
with mean , where is the scale parameter, is the shape parameter and is the gamma function. The true distribution of waiting times in a cycle model approaches the distribution if the rates are all the same and . In this limit , and describes the number of intermediate OFF states. In the more general case it correctly captures the effective properties of the process. The mean probability of the promoter to be in the ON state in the waiting time model is given by
The autocorrelation function cannot be computed directly analytically. The steady state Fourier transform of the steady state autocorrelation is (see Section E):
Finite cell cycle length correction to the connected autocorrelation function
Due to the short duration of the cell cycle, the theoretical connected correlation functions need to be corrected for finite size effects when comparing them to the empirically calculated correlation functions. When analyzing the data we calculate the autocorrelation function from traces of the same length , . We calculate the connected autocorrelation function for each trace and normalize it to at the second time point to avoid spurious nucleus to nucleus variability:
and then average over all traces to obtain the final connected autocorrelation function:
We define – the steady state true theoretical average of the random fluorescence intensity over random realizations of the process, and – the true theoretical second moment of the fluorescence signal. When the average over time points is equal to the theoretical average, . Using time invariance in steady state the autocorrelation function becomes:
where is an average over random realizations of the process. Eq. 14 corresponds to the limit we calculated in the theoretical models. To account for the finite size effects that arise due to short time traces, we need to correct for the fact that for short traces and , instead both the mean and the variance are functions of . We note that for short traces the definitions of autocorrelation and autocovariance differ:
In practice for the analyzed dataset we found that the finite size effects for the variance can be neglected, however the mean over time points is a bad approximation to the ensemble mean. We present the finite size correction to the mean below. For completeness we include the finite size correction for the variance in Section I, although we do not use it in the analysis due to its numerical complexity and small effect.
If the variance of the normalized fluorescence intensity over random realizations of the process is well approximated by the average over the time points, we can replace the denominator in Eq. 12 by and in steady state evaluate the mean connected autocorrelation function (see Section H for details):
where is the infinite-size steady state non-connected correlation function of the process (given in Eq. 6 for the two state model, Eq. The Cycle model for the cycle model and as the Fourier transform of Eq. The waiting time model for the model) and the average is over random realizations of the process. The mean and variance of the signal, and , provide a normalization factor that is constant for all time differences . We normalize the autocorrelation function setting the second term to and these terms are not needed for the inference. If then is proportional to .
The inference proceeds in three steps:
Step 1. Signal calibration. The intensity of the measured signal depends on a constant trace-dependent offset value , . To calibrate this offset we take the maximum expression to be the mean of the maximun expression over all traces in a given region . We take the mean of the maxima of the intensities rather than the absolute maximum of all signals to avoid errors from overestimating the maximum. The calibrated fluorescence signal used in the analysis is then . is directly calculated using Eq. 1. is a known function.
Step 2. Estimating parameter ratios. For the two state, three state cycle and models, the ratios of the rates can be estimated directly from the steady state mean fluorescence values using Eqs. 7 and 10. The Poisson model is uniquely defined by and does not require further parameter inference beyond Step 1.
Step 3. Estimating parameters. Using the estimate for the ratio of the rates, the ON and OFF rates are found by minimizing the mean squared error between the autocorrelation function calculated from the data (Eq. 13), and the model (Eq. Finite cell cycle length correction to the connected autocorrelation function with the theoretical prediction for the appropriate model: Eq. 6 for the two state, Eq. The Cycle model for the cycle model and the Fourier transform of Eq. The waiting time model for the model).
We thank T. Mora for helpful discussions. This work was supported by a Marie Curie MCCIG grant No. 303561(AMW), PSL IDEX REFLEX (ND, AMW, MC), ARC PJA20151203341 (ND), ANR-11-LABX-0044 DEEP Labex (ND), ANR-11-BSV2-0024 Axomorph (ND and AMW) and PSL ANR-10-IDEX-0001-02.
- Jaeger (2011) J. Jaeger, Cellular and Molecular Life Sciences 68, 243 (2011), ISSN 1420-9071.
- Gregor et al. (2014) T. Gregor, H. G. Garcia, and S. C. Little, Trends in Genetics 30, 364 (2014), ISSN 0168-9525.
- Driever and Nuesslein-Volhard (1988) W. Driever and C. Nuesslein-Volhard, Cell 54, 95 (1988).
- Tikhonov et al. (2015) M. Tikhonov, S. C. Little, and T. Gregor, Royal Society Open Science 2 (2015).
- Porcher and Dostatni (2010a) A. Porcher and N. Dostatni, Current Biology 20, R249 (2010a), ISSN 1879-0445.
- Little et al. (2013) S. C. Little, M. Tikhonov, and T. Gregor, Cell 154, 789 (2013), ISSN 1097-4172.
- Elowitz et al. (2002) M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain, Science 297, 1183 (2002), ISSN 1095-9203.
- Ozbudak et al. (2002) E. M. Ozbudak, M. Thattai, I. Kurtser, A. D. Grossman, and A. van Oudenaarden, Nature Genetics 31, 69 (2002), ISSN 1061-4036.
- Raser and O’Shea (2004) J. Raser and E. O’Shea, Science 304, 1811 (2004).
- Crauk and Dostatni (2005) O. Crauk and N. Dostatni, Current Biology 15, 1888 (2005), ISSN 0960-9822.
- Suter et al. (2011) D. M. Suter, N. Molina, D. Gatfield, K. Schneider, U. Schibler, and F. Naef, Science 332, 472 (2011), ISSN 0036-8075.
- Zoller et al. (2015) B. Zoller, D. Nicolas, N. Molina, and F. Naef, Molecular Systems Biology 11, 823 (2015), ISSN 1744-4292.
- Taniguchi et al. (2010) Y. Taniguchi, P. J. Choi, G.-W. Li, H. Chen, M. Babu, J. Hearn, A. Emili, and X. S. Xie, Science 329, 533 (2010), ISSN 1095-9203.
- Kandhavelu et al. (2012) M. Kandhavelu, J. Lloyd-Price, A. Gupta, A.-B. Muthukrishnan, O. Yli-Harja, and A. S. Ribeiro, FEBS Letters 586, 3870 (2012), ISSN 1873-3468.
- Muthukrishnan et al. (2012) A.-B. Muthukrishnan, M. Kandhavelu, J. Lloyd-Price, F. Kudasov, S. Chowdhury, O. Yli-Harja, and A. S. Ribeiro, Nucleic Acids Research 40, 8472 (2012), ISSN 1362-4962.
- Chong et al. (2014) S. Chong, C. Chen, H. Ge, and X. S. Xie, Cell 158, 314 (2014), ISSN 1097-4172.
- Lucas et al. (2013) T. Lucas, T. Ferraro, B. Roelens, J. D. L. H. Chanes, A. Walczak, M. Coppey, and N. Dostatni, Current Biology 23, 2135 (2013), ISSN 09609822.
- Garcia et al. (2013) H. G. Garcia, M. Tikhonov, A. Lin, and T. Gregor, Current Biology 23, 2140 (2013), ISSN 1879-0445.
- Ferraro et al. (2016a) T. Ferraro, T. Lucas, M. Clémot, J. De Las Heras Chanes, J. Desponds, M. Coppey, A. M. Walczak, and N. Dostatni, Wiley interdisciplinary reviews: Developmental Biology (2016a), ISSN 1759-7692.
- Lucas and et al (2016) T. Lucas and et al (2016).
- Bothma et al. (2014) J. P. Bothma, H. G. Garcia, E. Esposito, G. Schlissel, T. Gregor, and M. Levine, Proceedings of the National Academy of Sciences of the United States of America 111, 10598 (2014), ISSN 1091-6490.
- Estrada et al. (2016) J. Estrada, F. Wong, A. DePace, and J. Gunawardena, Cell 166, 234 (2016), ISSN 1097-4172.
- Coulon et al. (2014) A. Coulon, M. L. Ferguson, V. de Turris, M. Palangat, C. C. Chow, and D. R. Larson, eLife 3, 1 (2014), ISSN 2050-084X, eprint arXiv:1011.1669v3.
- Coulon and Larson (2016) A. Coulon and D. R. Larson, Methods in Enzymology 03, 1 (2016).
- Gregor et al. (2007a) T. Gregor, E. F. Wieschaus, A. P. McGregor, W. Bialek, and D. W. Tank, Cell 130, 141 (2007a), ISSN 0092-8674.
- Cheung et al. (2014) D. Cheung, C. Miles, M. Kreitman, and J. Ma, Development (Cambridge, England) 141, 124 (2014), ISSN 1477-9129.
- Gregor et al. (2007b) T. Gregor, D. W. Tank, E. F. Wieschaus, and W. Bialek, Cell 130, 153 (2007b), ISSN 0092-8674.
- Berg and Purcell (1977) H. C. Berg and E. M. Purcell, Biophysical Journal 20, 193 (1977), ISSN 0006-3495.
- Erdmann et al. (2009) T. Erdmann, M. Howard, and P. R. Ten Wolde, Physical Review Letters 103, 2 (2009), ISSN 00319007, eprint arXiv:0907.4289v1.
- Porcher and Dostatni (2010b) A. Porcher and N. Dostatni, Current Biology 20, 249 (2010b), ISSN 09609822.
- Perry et al. (2012) M. W. Perry, J. P. Bothma, R. D. Luu, and M. Levine, Current Biology 22, 2247 (2012), ISSN 09609822.
- He et al. (2011) F. He, J. Ren, W. Wang, and J. Ma, PloS one 6, e19122 (2011), ISSN 1932-6203.
- Fukaya et al. (2016) T. Fukaya, B. Lim, and M. Levine, Cell pp. 1–11 (2016), ISSN 1097-4172.
- Wang et al. (2011) Y. Wang, K. C. Tu, N. P. Ong, B. L. Bassler, and N. S. Wingreen, Biophysical Journal 100, 3045 (2011), ISSN 1542-0086.
- Lin et al. (2015) Y. Lin, C. H. Sohn, C. K. Dalal, L. Cai, and M. B. Elowitz, Nature 527, 54 (2015), ISSN 0028-0836.
- Dunlop et al. (2008) M. J. Dunlop, R. S. Cox, J. H. Levine, R. M. Murray, and M. B. Elowitz, Nature Genetics 40, 1493 (2008), ISSN 1546-1718.
- Munsky et al. (2012) B. Munsky, G. Neuert, and A. van Oudenaarden, Science 336, 183 (2012).
- Ferraro et al. (2016b) T. Ferraro, E. Esposito, L. Mancini, S. Ng, T. Lucas, M. Coppey, N. Dostatni, A. M. Walczak, M. Levine, and M. Lagha, Current Biology 26, 212 (2016b), ISSN 1879-0445.
- Gregor et al. (2016) T. Gregor, G. Tkacik, and et al (2016).
- Okabe-Oho et al. (2009) Y. Okabe-Oho, H. Murakami, S. Oho, and M. Sasai, PLoS Computational Biology 5, e1000486 (2009), ISSN 1553-7358.
- Bertrand et al. (1998) E. Bertrand, P. Chartrand, M. Schaefer, S. M. Shenoy, R. H. Singer, and R. M. Long, Molecular Cell 2, 437 (1998), ISSN 10972765.
- Janicki et al. (2004) S. M. Janicki, T. Tsukamoto, S. E. Salghetti, W. P. Tansey, R. Sachidanandam, K. V. Prasanth, T. Ried, Y. Shav-Tal, E. Bertrand, R. H. Singer, et al., Cell 116, 683 (2004), ISSN 00928674.
- Katsani et al. (2008) K. R. Katsani, R. E. Karess, N. Dostatni, and V. Doye, Molecular Biology of the Cell 19, 3652 (2008).
- Gillesple (1977) D. T. Gillesple, 93555, 2340 (1977).
- Bratsun et al. (2005) D. Bratsun, D. Volfson, L. S. Tsimring, and J. Hasty, PNAS 102, 14593 (2005).
Appendix A Autocorrelation analysis
a.1 Basic setup and data preprocessing
The raw data produced experimentally is a fluorescent signal measured at discrete times corresponding to the sampling time frame of the movie (see Fig. 13 for examples of traces). At each locus and at each time point it is the sum of the background signal and a number of fluorescent molecules attached to loops formed by the mRNA. Each loop contributes to the signal by a constant . This constant is unknown and can vary from trace to trace due to noise in the experimental setup and the variability in the locations of the nuclei in the embryo. All models are written for the renormalized signal .
Because the fluorescent signal is produced by discrete polymerases that travel down the gene, we divide the gene into chunks of base pairs, a length that corresponds to the irreducible space occupied by a polymerase on the gene (Fig. 4 in the main text). The positions the polymerase can occupy on the gene are labeled by an index . The number of MS2 loops that have been formed by a polymerase that has reached a given position depends only on the MS2 gene construct and we define a deterministic function for the whole length of the gene that describes the number of MS2 loops that have been produced by a polymerase at position . In practice the exact number of loops is not an integer and varies from base pair to base pair so we take as the average number of loops at this polymerase position (see Fig. 4 in the main text).
When the gene is fully loaded with polymerases (the number of polymerases is equal to the length of the gene divided by bp), the fluorescence intensity is . Assuming that the maximum of the signal over the whole trace is a good approximation for the fully loaded value we can determine and renormalize the data. In practice, since we see variability in the expressed signal in different nuclei at the same position, we are not sure the fully loaded polymerase scenario occurs in each nuclei, so we take the mean of the maximum intensity values in the anterior. We use this renormalized fluorescence signal to infer the parameters of the dynamics.
The experimental data is analyzed assuming the system is in steady state and does not take into account the initial activation period after mitosis, and the end of the trace when the gene is deactivated before mitosis. We take only a window of the traces where the mean spot intensity in all traces is stable (see Fig. 9 for an example). As the duration of the interphase differs slightly between embryos, we use a different steady state window for each embryo as summarized in Table 1.
In all models based on a stochastic gene switching (so all models except the Poisson-like model) we assume that the gene can be in several states with only two effective transcription rates: a non zero transcription rate in the ON state and an basal production rate equal to zero in the OFF state. When the gene is ON the polymerase loads at a maximal rate set by clearing of the binding site by the previous polymerase, which is one polymerase every seconds (calculated as the irreducible polymerase length along the gene bp divided by the polymerase speed, ). The state of the gene is described by a stochastic process that is equal to when the gene loads polymerase (i.e is ON) and when the gene is OFF (see Fig. 1B in the main text). Once the polymerase is loaded its path is assumed to be deterministic with constant speed.
The gene can be described by the locations where there is a polymerase: we define as a function of time and position that is equal to if there is polymerase at position at time and otherwise (see Fig. 1D in the main text). The fluorescence signal is then a convolution of the polymerase position, , and the details of the loop design of the MS2 construct, :
and the polymerase position can easily be translated back to the gene state through the deterministic relation, (see Fig. 3D in the main text for the form of ). This disruption is exact for a system with a discrete regulatory process and a discrete time step equal to the polymerase time step. Unfortunately, the moments in time when the gene switches are not necessarily multiples of the natural coarse graining steps of the system (the polymerase time step and its equivalent length) so it is necessary to introduce a continuous time in the system. We will present results for both the discrete and continuous time models. The continuous description is valid in the limit where the typical time spent by the gene in each state is long compared to the polymerase step or equivalently the gene switching constants are small compared to s. See Section A.2 for a more detailed argument.
|beginning (s)||end (s)||interphase duration (s)|
|Embryo 1 - cc 12||278||391||652|
|Embryo 1 - cc 13||318||546||796|
|Embryo 2 - cc 12||281||393||570|
|Embryo 2 - cc 13||278||484||695|
|Embryo 3 - cc 12||297||408||676|
|Embryo 3 - cc 13||330||587||824|
|Embryo 4 - cc 12||261||515||616|
|Embryo 4 - cc 13||321||529||751|
a.2 The two-state model
In this section we derive the equations required for the inference of the dynamics under the assumption that the gene can be in two states: ON or OFF, represented by a two dimensional vector . is the probability of the gene to be ON and is the probability of the gene to be OFF. is the average over traces of the random variable depicted in Fig. 1B of the main text. We assume that the switching times between the two are exponentially distributed:
The steady state probability to be ON is , where is the duration of the steady state window in Fig. 9, and is:
We learn from Eq. 1 in the main text:
and use it to obtain the ratio of the switching rates from Eq. 19.
The autocorrelation function is:
where the brackets are an average over traces (different realizations of the random process). We define – the probability that the polymerase is at position and time given that there was a polymerase at position at time (here we assume that ). Using the deterministic relation between the polymerase position at a given time and the probability to be on at an earlier time , is equivalent to the probability that the gene in ON at time given that is was ON at time :
Plugging the expression into Eq. 21 we obtain Eq. 2 in the main text:
In steady state the system is translationally invariant and for brevity we will denote it as - the probability that the gene is ON at time , given that it was ON at time . To find we need to solve for :
where is given by Eq. 18 and calculate the expectation value that the gene is ON at time given in was ON initially:
Eq. 25 is correct in a continuous time model. Its discrete time equivalent is
The eigenvalues of are , where with corresponding eigenfunctions:
The transition matrix is
In steady state and the connected autocorrelation is:
Since we already know the ratio of the rates from , inferring using Eq. 31 determines and .
a.3 Computing out of steady state
The autocorrelation approach can be generalized to a case when the system is out of steady state, when the autocorrelation function explicitly depends on the two time points and not only on their difference. During mitosis the gene is OFF and then gets turned ON in early interphase. Motivated by the hunchback expression we will present the calculation assuming the gene is initially OFF, but it is generalizable to any other initial condition. Assuming , we want to calculate the probability that the polymerase is at position at time , given that it was at position at time . Since the gene is initially OFF, we need to calculate the probability that the gene is ON at time . The autocorrelation function of the polymerase position is:
Using Eq. 30 and
a.4 Multiple off states
The calculations presented in Appendix A.2 can be extended to models that include more OFF or ON states as long there are only two production states for the mRNA: one enhanced and one basal production state. The transition matrix will then be of higher dimension and in practice should be (and has to be for dimensions larger than 3) diagonalized numerically. The exact analytical solution for the autocorrelation function is still valid written in terms of the powers of .
a.5 Generalized multi step model
A gene with many OFF states can also be described using a reduced model with two effective gene expression states ON and OFF, where the times of transitions between these two state are not exponential but follow a peaked distribution approximated by a Gamma distribution. The Gamma distribution describes an effective transition over many irreversible transitions between a series of OFF states:
where is the scale parameter, is the shape parameter, and is the gamma function. The mean time spent in the OFF state is , so the probability for the gene to be in the ON state is:
This model has three parameters, regardless of the number of OFF states, and using Eq. 1 of the main text reduces the number of parameters to two, which greatly simplifies the inference. The remaining two parameters are learned from the autocorrelation function in Eq. 21, which formally has the same form as Eq. 23:
but is now the two-point correlator of a non-Markovian process. We limit our presentation to the steady state, but the calculation generalizes to out of steady state systems.
We cannot solve the problem in real space, but we compute the Fourier transform of the autocorrelation function of the fluorescence signal:
which using Eq. 21
we reduce to calculating
We decompose into a sum over full cycles of the gene turning from ON to OFF, with the constraint that at time the gene is ON: