Stochastic models of gene transcription with upstream drives: exact solution and sample path characterization

Abstract

Gene transcription is a highly stochastic and dynamic process. As a result, the mRNA copy number of a given gene is heterogeneous both between cells and across time. We present a framework to model gene transcription in populations of cells with time-varying (stochastic or deterministic) transcription and degradation rates. Such rates can be understood as upstream cellular drives representing the effect of different aspects of the cellular environment. We show that the full solution of the master equation contains two components: a model-specific, upstream effective drive, which encapsulates the effect of cellular drives (e.g. entrainment, periodicity or promoter randomness) and a downstream transcriptional Poissonian part, which is common to all models. Our analytical framework treats cell-to-cell and dynamic variability consistently, unifying several approaches in the literature. We apply the obtained solution to characterize different models of experimental relevance, and to explain the influence on gene transcription of synchrony, stationarity, ergodicity, as well as the effect of time scales and other dynamic characteristics of drives. We also show how the solution can be applied to the analysis of noise sources in single-cell data, and to reduce the computational cost of stochastic simulations.

Keywords: chemical master equation; gene expression; noise; non-stationarity; stochastic models; transcription.

Figure 1.

Single-cell gene transcription under upstream drives. The transcription of each cell i takes place under particular cellular drives and , representing time-varying transcription and degradation rates. Both cellular drives are combined into the upstream effective drive , which dictates the long-term probability distribution describing the stochastic gene expression within each cell (2.10). When there is cell-to-cell variability in the population, the cellular drives are described by processes M and L leading to the upstream effective drive X. The probability distribution of the population corresponds to the mixture of the upstream process X and the Poissonian downstream transcriptional component, as given by (2.14). Increased synchrony in the population implies decreased ensemble variability of the random variables M_t, L_t, X_t and N_t. (Online version in colour.)

Figure 2.

Gene transcription under the RE model (3.1) with constant degradation rate and transcription rates entrained to an upstream sinusoidal signal with ω = 1. Each cell has a random phase offset ϕ drawn from a distribution. (a) The synchronous population corresponds to identical phases across the population. In this case, the transcription reflects the time variability of the upstream drive mixed with the stochasticity due to the downstream Poisson process of transcription. When the random phases ϕ are uniformly distributed on an interval of range (b) π and (c) 2π, the population becomes increasingly asynchronous. For all three cases, we show (top to bottom): sample paths of the effective drive, X; its density given by equation (3.2); sample paths of the number of mRNAs, N; and the full solution of the ME P(n, t). (Online version in colour.)

Figure 3.

Asynchronous stochastic promoter switching models correspond to upstream stochastic processes. The promoter cycles between the discrete states, transitioning stochastically with rates as indicated: (a) the standard 2-state random telegraph model; (b) the 3-state refractory promoter model. (Online version in colour.)

Figure 4.

(a) Modulated random telegraph model: each cell switches asynchronously between ‘ON’ and ‘OFF’ states, but the magnitude of the ‘ON’ transcription rate is modulated by the function ρ(t; ϕ), a sinusoid representing an upstream periodic process. The phase ϕ represents the cell-to-cell variability and leads to the varying degree of synchrony across the population. (b) Sample paths μ_i(t) and solution of the probability distribution P(n, t) of the MRT model for synchronous (left) and asynchronous (right) modulation. In the asynchronous case, the upstream drive has a random phase across the cells with distribution . (Online version in colour.)

Figure 5.

Analysis of single-cell temporal transcription of a gene in response to an upstream oscillatory cAMP signal, motivated by recent experiments [29]. Individual single-cell time courses ν(t) of mRNA counts are highly variable with no clear entrainment to the driving signal, whereas the time-dependent ensemble average oscillates with the same frequency as the external drive. This is consistent with equations (4.5)–(4.6), which also show that the total phase lag is the resultant of the signal transduction and transcription lags. For a signal with period T = 5 min and a gene with degradation rate λ = 0.04 min⁻¹ [29], the transcription phase lag is , which corresponds to a delay of . Given a measured total mean lag of 9π/10, this implies that the signal transduction introduces a phase lag , equivalent to a delay of . (Online version in colour.)

Figure 6.

Noise characteristics of the Kuramoto promoter model (3.4). (a) Numerical simulations for C = 100 oscillatory cells and different coupling parameters: K = 0.002, 0.1, 0.4 ((i), (ii), (iii), respectively). For each coupling, the sample paths of the upstream effective drive X and mRNA counts N are shown. The mean, variance and ensemble Fano factor of N were calculated from the sample paths of N (blue lines) and, more efficiently, from the sample paths of X (black lines). The last row shows the Kuramoto order parameter r(t) measuring the cell synchrony, signalled by . (b) Ensemble Fano factor (averaged over the simulated time courses) against coupling parameter . As K is increased, the oscillators become synchronized and the ensemble Fano factor decreases towards the Poisson value of unity. (c) Scatter plot of the ensemble Fano factor against the order parameter r(t) (both averaged over the simulated time courses). As the oscillators become synchronized (), the ensemble Fano factor also approaches 1, signifying that the distribution is Poisson at all times. (Online version in colour.)

Figure 7.

‘Burn-in’ transient in the RT model. (a) Sample paths of the transcription rate M, the effective upstream drive X and the number of mRNAs N for an initial condition P(0, 0) = 1 with all cells initialized in the inactive state [16]. (b) The full solution of the RT model for this initial probability distribution exhibits an exponential decay as the system approaches its stationary distribution. The delta distribution at t = 0 is omitted for scaling purposes. (Online version in colour.)

Figure 8.

Ergodic transcription models under periodic and stochastic upstream drives. (a) We consider gene transcription under three drives : a sinusoidal wave with period T (yellow); a square wave with period T (red); a RT process with expected waiting time T/2 in each state (blue). For such ergodic systems, the distribution computed over time corresponds to the stationary distribution. (b) The distribution presents distinct features as the period T is varied. (Online version in colour.)

Figure 9.

The temporal Fano factor. (a) A sample path of mRNA counts from the (leaky) RT model. The time periods when the gene is in the active state are shaded. (b) The temporal Fano factor (5.4), , computed over a time window W of fixed length indicated by the horizontal bars at each t. When W extends over a stationary section of the sample path, TFF is close to unity, corresponding to the Poisson distribution (black dashed line). (c) Heat map of the cTFF (5.5), , defined only for t ≥ t₁. Note the marked step pattern corresponding to the switching times, indicated by dashed lines as a guide to the eye. (Online version in colour.)

Figure 10.

Efficient sampling of the full distribution P(n, t) for transcription with upstream cellular drives. We consider upstream drives governed by the Kuramoto promoter model (3.4) for C = 10 000 coupled oscillatory cells. Sample paths of N are simulated directly with the Gillespie algorithm to approximate P(n, t) at time t₁ (bottom, blue). Alternatively, sample paths of X are used to estimate , which is then mixed by performing the numerical integration (2.14) to obtain P(n, t) (top, red). The latter sampling through X is more regular and far less costly: CPU time via N is ≈36000 s, whereas CPU time via X is ≈0.1 s. (Online version in colour.)