One of the greatest problems of molecular biology is how single, undifferentiated cells give rise to many different cell types, all being genetically identical yet performing different functions. Since the pioneering work of Jacob and Monod [1] it is known that this multiplicity of behaviors is possible since the protein production depends not only on the existence of a gene but also on the quantities of regulatory molecules that can change with the cell type and environmental conditions. The protein production takes place in two steps [2]. First, during a process called transcription, the genetic information from DNA is copied to the messenger RNA (mRNA). Then, during a process called translation, the mRNA is used as template for protein production by the ribosomes. The amount of transcribed mRNA and translated protein, namely the gene expression, can vary from one cell to another for various reasons. One of these reasons is that the gene expression is not a property of a single gene but is a property of a set of interacting genes; a gene network [3]. As per usual in complex systems, the whole is different from the sum of its parts. In this case, a gene network can have many different stable expression levels that correspond to different network attractors [4]. Another reason for expression variability, is the fact that transcription and translation are stochastic [5–7]. A single cell transcription is often intermittent, periods of strong mRNA production being followed by periods of silence when transcription is stopped [7–10]. The periods of inactivity may correspond to paused RNA polymerase, incompletely formed activation complexes, transient presence of regulatory proteins and complexes on the DNA, or to transitions between open and condensed chromatin [11]. The durations of the periods of activity and inactivity are random. Furthermore, snapshots of the cell population at various times show a heterogeneous picture in which a significant proportion of the cell population express much lower or much higher than the average protein levels [6, 12]. The observed distributions of the gene expression can deviate significantly from a Gaussian by having skewness [12], heavy tails, or multimodality [11].

The gene expression stochasticity can have important biological consequences. In infections and tumors, a part of the pathogen population (bacteria, viruses, cancer) can escape drug treatment when they adopt behaviors very different from the average; for instance they may stop growth [13, 14], remodel their metabolism [15], or stop absorbing, or expel drugs [16]. Latency phenomena [17] can be responsible for resistance to treatment and relapse of the disease once the treatment is stopped.

Gene expression stochasticity is also important in synthetic biology. Synthetic biology devices are supposed to have a well-defined biological function for given combinations of the input conditions. Therefore, in the absence of error-correction, the precision of biological devices depends critically on the degree of stochasticity [18, 19].

For all these reasons, one needs mathematical tools for predicting quantitatively the amplitude and distribution of gene expression fluctuations in gene networks.

Models of stochastic gene networks represent molecular interactions as discrete events (biochemical reactions) separated by random waiting times. This modeling approach, introduced by Delbruck [20] and further developed by Renyi [21] and Bartholomay [22] covers practically all the aspects of gene-gene interactions but is mathematically and computationally challenging. Indeed, the underlying models are continuous time Markov processes with an infinite or extremely large number of states. The corresponding master equation can be solved exactly only in a limited number of cases, corresponding to single genes without or with much simplified feed-back control [23–26]. The Monte-Carlo simulation algorithms (Gillespie algorithm [27]) can be very inefficient for simulating the full process, as the number of individual reactions that have to be simulated for a significant change of the system's state can be tremendously large.

Several types of approximation were used in order to simplify stochastic biochemical networks in order to reduce their simulation time and to facilitate their analysis. These simplifications were possible because the stochastic gene networks have heterogeneous variables and multiple time scales [28]. This heterogeneity comes from the fact that some variables X_D such as DNA/regulatory proteins and complexes/polymerase states are discrete and other variables X_C such as protein and mRNA copy numbers are continuous. The species dichotomy leads to a partition of the biochemical reactions. Although literature provides a number of different ways for reaction partitioning, the four set partition R_D, R_C, R_DC, R_CD (Figure 1) seems to us quite natural. A very similar partitioning was used for rigourously justified approximations in Crudu et al. [29, 30].

Two main timescales have to be taken into account in order to find the appropriate approximation. As shown in the Figure 1, the discrete variables switch between a number of discrete states. The characteristic time of this process was called switching time, τ_S [28]. The trajectories of continuous variables are smooth only in the average, but for these variables, the average is a good approximation. The characteristic time of fluctuations of continuous species around their average was named discreteness time, τ_D [28]. The discreteness time scales like 1/N where N is the copy number of the continuous species.

There are two main classes of approximated models (Figure 2):

In this paper we consider the situation when the swiching time is relatively large. In biology, this corresponds to the so-called “random telegraph” [36], “bursting” [8, 12], or “multi-scale bursting” [10] fluctuations. Simply speaking there are some variables (DNA, regulatory complexes, and/or RNA polymerase states) that have ON/OFF or multiple state Markov chain dynamics. The discrete variables control the ODE dynamics of the continuous variables (mRNA, proteins). The underlying approximation is piecewise-deterministic, because the ODEs change each time that the discrete variables perform a transition (see Figure 1).

In our PDP models of gene networks each gene promoter is described as a finite state Markov chain. The transition rates between states of the promoter depend on the expression levels of the regulating genes. The promoter triggers synthesis of gene products with intensities depending on its state. We represent the state of a gene network as a random vector whose components indicate the promoter state and the copy numbers of all the mRNAs and proteins in the network. The dynamics of the gene network will be represented as the time-dependent multivariate distribution of this random vector. These distributions contain information about the randomness of each gene but also about correlation among genes and promoter states. We present three methods to compute time dependent multivariate distributions for such models: the direct simulation of the PDP process (the Monte-Carlo method); the numerical solutions of the Liouville-master equation, and the push-forward method. The first method is stochastic, whereas the last two are entirely deterministic (do not use random number generators). The PDP Monte-Carlo method combines the numerical simulation of the Markov chain of promoter states with symbolic solutions of the deterministic ODEs and is much faster than the direct Gillespie method [27].

This paper is structured as follows: In section 2 we recall the justification of piecewise-deterministic approximation using the partial Ω expansion. In section 3 we introduce a class of piecewise-deterministic models covering gene networks. In section 4.1 we discuss the Monte-Carlo methods for piecewise-deterministic models. In section 4.2 we discuss the Liouville-master equation approach. In section 4.3 we discuss the push-forward method. In section 5 we briefly discuss possible applications of these methods to extracting information from mRNAs' and proteins' spectra of fluctuations.

The dynamics of stochastic biochemical networks can be described by a pure jump Markov process [20–22] (For a general presentation of Markov jump processes see Gikhman and Skorokhod [37]). The state of the network is a vector X ∈ ℕⁿ whose components represent copy numbers of molecules of various species. Each biochemical reaction is defined by a stoichiometric (jump vector) γi∈ℤn,i∈R where R is the set of biochemical reactions in the model. The occurrence of the reaction i is represented as a jump in the system's state X → X + γ_i. Finite, discrete states of gene promoters can be represented by extending the meaning of species to include “places” with finite values of the copy numbers. For instance, on/off gene promoters can be represented by using two places P_on and P_off with two possible occupancies 0 or 1. The transitions from OFF to ON and back can be represented by a reversible reaction that consumes a particle on the place P_off and produces a particle on the place P_on.

In the Markov pure jump representation the dynamics of the system is a series of jumps separated by exponentially distributed waiting times [35, 37]. The number of reactions of the type i occurring in the average per unit time is given by the propensity (or rate) function X → V_i(X; μ), where μ are kinetic parameters. The time between successive reactions is exponentially distributed with an average (∑i∈RVi)-1 and the next reaction is i with probability Vi/(∑i∈RVi). We are interested in the multivariate distribution p(X, t) representing the probability to be in state X at time t. This satisfies the time dependent master equation:

The van Kampen Ω (system size) expansion [33] or, equivalently, the central limit theorem [31, 38] lead, in the first order, to deterministic (ODE) approximations and to diffusion (Fokker Planck) approximations, in the second order. The partial omega expansion consists in applying the Ω expansion only to species that are produced in sufficiently large copy numbers [29, 30]. The large copy numbers species denoted X_C are called continuous, because the biochemical reactions change their values gradually, by the accumulation of a large number of small steps. Other species are present only in a few copies per cell (these include the “places” describing promoter states). We denote these species X_D and call them discrete. This leads to a decomposition of the state vector as X = (X_D, X_C) and of all stoichiometric vectors as γi=(γiD,γiC), corresponding to discrete and continuous species coordinates. The interactions among discrete and continuous species are suitably described by a partition of the reactions in four sets R = R_D ∪ R_DC ∪ R_CD ∪ R_C [28–30].

The reactions R_D act on X_D (the corresponding γ_i have non-zero coordinates on discrete species, γiD≠0) and have propensities depending on X_D only. The reactions R_C act on X_C (the corresponding γ_i have non-zero coordinates on continuous species, γiC≠0) and have propensities depending on X_C only. The reactions R_DC, R_CD act on X_C and X_D, respectively, and their propensities depend on both X_D and X_C.

In this paper we consider gene network models. For each gene, we model the transitions between promoter states, as well as other processes such as transcription, translation, protein folding, and protein and mRNA degradation. We will consider that the mRNA molecules and proteins are in sufficiently large copy numbers to justify continuous approximations. The only discrete variables are in this case the promoter states. The set R_D contains transitions between discrete promoter states whose rates do not depend on regulatory proteins. The set R_CD contains transitions between promoter states whose rates depend on concentrations of regulatory proteins. The set R_C contains translation, maturation (folding), degradation reactions. The set R_DC contains transcription initiation reactions that depend on the promoter state.

We further consider that the copy numbers of continuous species X_C and the propensities of reactions in the sets R_C and R_DC are “extensive,” in other words, scale with the system size Ω, X_C = Ωx_c, V_i = Ωv_i(x_c, X_D), for i ∈ R_DC ∪ R_C, whereas the propensities of reactions in R_D and R_CD are considered “intensive” and do not scale with Ω. The Ω dependence is not only a useful mathematical tool, but has also a biological meaning. In proliferating cells, the protein concentrations are important for biochemical reactions and should be maintained by synthesis reactions. This is only possible if the propensities of synthesis reactions (including mRNA synthesis) scale with the size. Rates of monomolecular reactions consuming continuous reactants (for instance, degradation reactions) are proportional to the reactant copy number that scale with the size. Rates of switching reactions between discrete promoter states do not scale with size, unless they are proportional to copy numbers of activator or repressor proteins. For a more complete discussion of these scaling relations we refer to Crudu et al. [29].

In rescaled variables, the master equation reads

Using the first order Taylor series expansion vi(XD,xc-γic/Ω;μ)p(XD,xc-γiC/Ω,t)-vi(X,xc;μ)p(XD,xc,t)≈-γicΩ∂vi(X,xc;μ)p(XD,xc,t)∂xc we obtain the Liouville-master equation

The Liouville-master equation has been used in various fields; in statistical physics of quantum systems [39–41] or with a different name in statistics and operations research [42]. It describes the time-dependent distribution of a piecewise-deterministic model. Indeed, conditionally on X_D (considering that the discrete variables are fixed) the continuous variables x_c satisfy deterministic, ODEs:

The discrete variables follow a pure jump Markov dynamics defined by the stoichiometric vectors γiD and the propensities V_i(X_D, x_c; μ), i ∈ R_D ∪ R_CD. The noise in the system is produced by the stochastic transitions of the discrete variables. The probability distribution of the continuous variables is advected (transported by the flow defined by the ODEs, see Figure 3) by the deterministic flow Φ(X_D, x_c; μ) in the continuous variables space; this explains the advection term -∂∂xc[Φ(XD,xc;μ)p(XD,xc,t)] in the Liouville-master equation.

In this section we introduce a family of PDP models that can be used to represent gene networks. They are a special, simplified case of the class of models defined by (2.3). The main simplification is that the propensities of reactions in R_DC depend on X_D and do not depend on x_C.

As an example, we consider a gene circuit with dichotomous noise. This model is made of n_g genes, each one controlled by a promoter with 2 states; ON and OFF. The continuous variables are x_i and y_i; the protein and mRNA level for each gene i, respectively. The discrete variables are two values variables s_i ∈ {0, 1} representing the promoter states (0 stands for OFF and 1 stands for ON). In this model, the set R_DC consists of transcription initiation reactions. Their rates depend on the promoters states (ON or OFF) but do not depend on the protein or mRNA levels. The set R_CD ∪ R_D consists of reversible transitions between the states ON and OFF. The corresponding transition rates are constant when the gene i is constitutive; namely, these rates are f_i from ON to OFF and h_i from OFF to ON. When the gene i is activated by the gene j the transition rate from OFF to ON is f_ix_j, whereas when the gene i is repressed by the gene j the transition rate from ON to OFF is h_ix_j. The discrete variables' dynamics are thus a Markov chain with the state set M={0,1}ng. It is convenient to relabel the states from 1 to 2ng using the lexicographic order on M. For instance, two gene circuits have, in order, the states 1:(0, 0), 2:(0, 1), 3:(1, 0), 4:(1, 1). Also, instead of using reaction propensities, we equivalently provide a transition rate matrix S whose elements are the transition probabilities per time unit between states. For instance, for a two gene circuit where the first gene is constitutive and the second gene is activated by the first, we have:

The mRNA and protein variables follow ODE dynamics

where ρ_i, b_i, a_i, i ∈ [1, n_g] represent mRNA degradation, translation, and protein degradation rates for the gene i, respectively, and s = (s₁, s₂, …, s_ng) is the state of the Markov chain, such that

For the sake of illustration let us consider the simple model of a single constitutive gene controlled by a two state (ON/OFF) promoter. We denote the states of the promoter by 1 and 0, respectively. The transition rate from 0 to 1 is f and from 1 to 0 is h. The protein and mRNA concentrations are x and y, respectively. The transcription initiation rate in the state 1 is k₁ and in the state 0 is k₀ < < k₁. The translation rate is b. The mRNA and protein degradation rates are ρ and a, respectively. The master-Liouville equation reads

The protein and mRNA concentrations follow the ODEs

The probability distribution of the promoter state s results from the dynamics of the two state Markov chain

where p₀ = ℙ[s = 0] = ∫p(0, x, y)dxdy, p₁ = ℙ[s = 1] = ∫p(1, x, y)dxdy.

We also define the asymptotic occupancy probabilities p0=hh+f and p1=fh+f, representing the probabilities, at steady state, that the promoter state is OFF and ON, respectively.

The single constitutive gene model and the advection fluxes of the Liouville-master equation are illustrated in Figure 3. More complex, two gene circuits models are represented in Figure 4 and their Liouville-master equations are given in the Appendix 1.

In this section we compare several numerical methods for solving the Liouville-master equation in the context of gene networks models. In order to quantify the relative difference between methods we use the L¹ distance between distributions. More precisely, if p(x) and p~(x) are probability density functions to be compared, the distance between distributions is

In this paper we used mathematical models to predict mRNA and protein “intrinsic” fluctuations (by “intrinsic” we understand fluctuations generated by the stochastic gene network dynamics). An important question in this context is if “intrinsic” gene expression fluctuations can be used to reverse engineer gene networks. “Extrinsic” fluctuations (by “extrinsic” we understand perturbations of gene networks coming from their environment) of the transcriptome were extensively used in the past to reconstruct gene networks using correlation or mutual information (for a popular method see [48]). A few reverse-engineering studies obtained gene regulation parameters from intrinsic gene expression fluctuations [12, 49]. Quantifying intrinsic fluctuations requires single cell measurements of several genes. A variety of technologies are now ready to take this challenge—single cell sequencing [50], single-cell RNA microscopy [10], various versions of time-lapse microscopy [51], fluorescence correlation microscopy [12]. It is therefore important to test the propensity of expression fluctuations to discriminate between various gene network architectures.

The predictions of our PDP models are provided as multivariate distributions of mRNAs and proteins copy numbers for one or several genes. These predictions can be directly compared with results obtained with molecular biology and biophysics experimental methods.

First, one would like to test if the differences between distributions predicted with various gene circuit models are significant and therefore can be used to discriminate between gene circuit models. To this end, we generated bivariate proteins and mRNA distributions for five different gene circuits like in Figure 4. The visual inspection of Figure 8, suggests that mRNA distributions resulting from five different gene circuits, with gene-gene interactions that differ by their signs, are very similar in the same regime of switching (fast or slow) of the genes. The mRNA distributions discriminate between model parameters (fast or slow switching) but do not discriminate between circuit architectures. The protein bivariate distributions are shown in Figure 9. They differ strongly from mRNA distributions and discriminate between both parameters and architectures.

As visual colormap differences may be judged subjective, we developed a quantitative test for the discriminant power. This test is based on the distance d defined by (4.1). We have computed d pairwise, for mRNA and for protein distributions produced from different gene circuits at steady state. In order to test if these distances are significantly large we compared them to distribution of distances between random samples generated by the same gene circuit model for a fixed number of cells. The result of this comparison is shown in Figure 10 for the two gene circuits G₁→G₂ and G₁⊣G₂ that differ by the sign of the interaction; one can notice that the protein fluctuation based distance is significant, whereas the mRNA fluctuation distance is not, both for slow/slow and fast/fast genes.

We have also tested the significance of the correlation computed from bivariate mRNA or protein distributions. A simple gene reconstruction method is to consider that genes interact if the correlation coefficient is significantly different from zero. We have computed the Bravais-Pearson correlation coefficient from bivariate mRNA and protein samples generated with our PDP model, at steady state and for increasing numbers of cells N_c. For both models G₁→G₂ and G₁⊣G₂ a significant (upper tail probability p < 5%) protein-protein correlation is obtained for moderate cell populations (N_c > 100 for p < 5%, see Figure 10). In order to obtain significant mRNA/mRNA correlation one has to use very large numbers of cells (N_c > 1, 000 for p < 5%, see Figure 10). This is possible for single cell sequencing and flow cytometry (with the drawback of the lack of precision in estimating the mRNA copy numbers) but is very difficult for techniques such as MS2 tagging microscopy, or time-lapse microscopy.

We have discussed three methods to compute time-dependent distributions of mRNA and protein copy numbers generated by gene networks. All the three methods are much faster than the Gillespie exact chemical master equation. The finite difference Liouville-master equation method is precise and fast for small models. Simple (non-adaptive) finite difference schemes, however, are demanding in terms of space and time resolution leading to computer memory limitations. In future implementations of the Liouville-master equation method we plan to use spectral methods for bypassing these limitations. The push-forward method is not as precise as the Liouville-master equation, but it is much more stable, even at low resolution. The differences between execution times of various methods are illustrated in the Table 1 indicating that the push-forward method is the fastest. The performance of the push-forward method results from the reduced cardinality of the discrete phase space (2 states for one gene, 4 states for 2 genes) which is an improvement with respect to the initial version in Innocentini et al. [47]. Further improvements, lifting the mean-field approximation for gene coupling will be presented elsewhere.

Piecewise-deterministic models are valuable tools for understanding stochastic gene expression in a wide spectrum of regimes, covering both slow and fast switching. The source of stochasticity in such models are the random transitions between discrete promoter states; a phenomenon usually associated with transcriptional bursting. In this paper we have only discussed dichotomous noise (ON/OFF promoters), however, as seen in section 2, our methods work also for promoters with more than two states.

As application of our numerical methods we tested the capacity of mRNA and protein copy numbers fluctuations to unravel differences between gene circuit architectures. We showed that protein fluctuations are sensitive to differences of architecture and that protein-protein correlation reveals gene-gene interactions for moderate cell population sizes (100 cells). In contrast, mRNA fluctuations are much less sensible to differences in circuit architecture and mRNA-mRNA correlation is small, even for interacting genes. This reinforces the already well established conclusion that proteome contains much more information than the transcriptome. In the past we used the spectrum of protein copy number fluctuations to extract information about promoter repression mechanisms [12]. The difference in behavior of the mRNA and protein fluctuations can be explained by the fact that the mRNA half life is usually much shorter than the protein half life. Gene-gene interactions are mediated by proteins that slowly modulate the gene switching times. Proteins follow these slow modulations, which results in significant protein-protein correlation. mRNA dynamics are permanently submitted to the faster (uncorrelated) gene switching, which explains the lower correlation of steady state mRNA fluctuations. This suggests that reconstruction of gene networks from mRNA intrinsic fluctuations is submitted to severe limitations. More general results concerning these limitations will be presented elsewhere.

OR and GI conceived the research. OR implemented the Monte-Carlo method. GI implemented the push-forward method. AH implemented the Liouville-master method. All authors were involved in writing the paper. All the authors have read and approved the final manuscript.