Skewness and Kurtosis of mRNA Distributions in Stochastic Gene Transcription with Promoter Switching

1. Introduction

Gene expression is an inherently stochastic process, arising from the small copy numbers of DNA templates, transcription factors, and mRNA molecules involved in transcriptional regulation. Early theoretical and experimental studies have established that intrinsic molecular fluctuations can lead to pronounced cell-to-cell variability, even in genetically identical populations under homogeneous environmental conditions [1,2,3,4]. With the advent of single-cell and single-molecule experimental techniques, such variability has been directly observed at the mRNA level, revealing substantial heterogeneity in transcriptional output across individual cells [5,6,7]. Quantifying and interpreting this variability has therefore become a central problem in systems biology and mathematical biology.

Variance and the Fano factor are among the most widely used statistical measures for characterizing transcriptional noise. These quantities effectively describe the overall dispersion of gene expression levels and have played an important role in elucidating the relationship between mean expression and noise [8,9,10,11]. However, variance-based measures provide only limited information about the geometric structure of the underlying mRNA copy number distributions. In particular, they are insensitive to distributional asymmetry and tail behavior, which naturally arise from transcriptional bursting driven by stochastic promoter switching dynamics [12,13,14,15]. A detailed understanding of the shape of mRNA distributions is therefore essential for uncovering the molecular mechanisms governing gene regulation and for accurately interpreting experimental data.

In principle, the full probability mass function of mRNA copy numbers contains complete information about distributional geometry. For simple stochastic gene expression models, exact or approximate stationary distributions can sometimes be obtained [15,16]. However, for biologically realistic models with explicit promoter switching, feedback regulation, or time-dependent synthesis rates, analytical expressions for the full distribution are often unavailable or unwieldy [17,18]. Even when closed-form solutions exist, they do not necessarily yield intuitive or easily comparable summaries for experimental analysis. This motivates the introduction of reduced statistical descriptors that capture essential geometric features of the distribution while remaining analytically tractable and experimentally accessible.

In this work, we focus on two such higher-order statistical measures: skewness and kurtosis. Skewness quantifies the degree of asymmetry of the mRNA distribution and indicates whether transcriptional fluctuations are biased toward rare, high-expression events or toward low-expression states. Kurtosis characterizes the sharpness of the distribution peak and the heaviness of its tails relative to a Gaussian distribution, thereby reflecting the prevalence of extreme transcriptional events. Recent studies have demonstrated that these higher-order moments provide sensitive indicators of transcriptional bursting and gene expression heterogeneity, particularly in large single-cell datasets [6,8].

The skewness is defined as the third standardized moment

$$\begin{array}{c}\hfill \mathrm{Skew}\left[X\right]=\mathbf{E}\left[{\left({\displaystyle \frac{X-\mu }{\sigma }}\right)}^3\right],\end{array}$$

(1)

where $\mu$ and $\sigma$ denote the mean and standard deviation of the random variable X, respectively. The kurtosis is defined as the fourth standardized moment

$$\begin{array}{c}\hfill \mathrm{Kurt}\left[X\right]=\mathbf{E}\left[{\left({\displaystyle \frac{X-\mu }{\sigma }}\right)}^4\right]-3.\end{array}$$

(2)

By construction, these quantities are dimensionless and therefore enable direct comparison across genes, experimental conditions, and modeling frameworks.

The remainder of this article is organized as follows. In Section 2, we introduce a two-state stochastic gene transcription model with explicit promoter switching and present its mathematical formulation. Section 3 is devoted to the analytical derivation of skewness and kurtosis for the model. In Section 4, we perform numerical simulations and apply our theoretical results to open-source experimental data to investigate how skewness, kurtosis, and the coefficient of variation depend on key kinetic parameters. Finally, concluding remarks and biological implications are discussed in Section 5.

2. The Model and Expressions of Skewness and Kurtosis

2.1. Two-State Gene Transcription Model

We consider a stochastic two-state gene transcription model, which has been widely used to characterize promoter-mediated transcriptional bursting and has also been adopted in our previous studies. A schematic representation of the model is shown in Figure 1.

In this framework, the promoter of a gene stochastically switches between two discrete states: an inactive (OFF) state and an active (ON) state. Transitions from the OFF state to the ON state occur at a constant rate $\lambda$, whereas transitions from the ON state back to the OFF state occur at a rate $\gamma$. As a result, the residence times in both promoter states are exponentially distributed. When the promoter is in the ON state, mRNA molecules are synthesized at a constant transcription rate $\nu$, while no transcription occurs when the promoter is in the OFF state. Independent of the promoter dynamics, each mRNA molecule is degraded at a constant rate $\delta$. This minimal model captures the essential features of transcriptional bursting induced by stochastic promoter switching.

Let $M\left(t\right)$ denote the total number of mRNA transcripts present in a cell at time t, where $M\left(t\right)$ takes values in the set of nonnegative integers. We introduce an indicator variable $I\left(t\right)$ to represent the promoter state at time t, with $I\left(t\right)=0$ corresponding to the gene OFF state and $I\left(t\right)=1$ corresponding to the gene ON state. The joint probability distributions are defined as

$$\begin{array}{cc}\hfill P_0(m,t)& =\mathrm{Prob}\left\{M\left(t\right)=m,I\left(t\right)=0\right\},\hfill \\ \hfill P_1(m,t)& =\mathrm{Prob}\left\{M\left(t\right)=m,I\left(t\right)=1\right\}.\hfill \end{array}$$

The time evolution of the joint probabilities $P_0(m,t)$ and $P_1(m,t)$ is governed by the following system of chemical master equations:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{P}_0(m,t)}{dt}}=& \gamma P_1(m,t)-\lambda P_0(m,t)+\delta \left[(m+1)P_0(m+1,t)-m{P}_0(m,t)\right],\hfill \\ \hfill {\displaystyle \frac{d{P}_1(m,t)}{dt}}=& \lambda P_0(m,t)-\gamma P_1(m,t)+\delta \left[(m+1)P_1(m+1,t)-m{P}_1(m,t)\right]\hfill \end{array}$$

(2.1)

$$\begin{array}{cc}\hfill & +\nu \left[P_1(m-1,t)-P_1(m,t)\right],\hfill \end{array}$$

(2.2)

for $m\ge 0$, with the convention that $P_1(-1,t)=0$.

The marginal probability mass function of the mRNA copy number is then given by

$$P(m,t)=P_0(m,t)+P_1(m,t),$$

which represents the probability that m mRNA molecules are present in the cell at time t.

The terms in Eqs. (2.1)–(2.2) account, respectively, for stochastic switching of the promoter state, mRNA synthesis in the active state, and first-order degradation of mRNA molecules. This coupled system of master equations fully characterizes the stochastic dynamics of the two-state transcription process.

In principle, the probability mass function $P(m,t)$ provides a complete description of the distribution of mRNA copy numbers across a cell population. However, for transcription models of even moderate complexity, obtaining explicit analytical expressions for $P(m,t)$ is often intractable. Even for the two-state model considered here, closed-form solutions are available only under restrictive assumptions and are typically expressed in terms of integrals or special functions. Consequently, rather than relying exclusively on the full probability mass function, we focus on higher-order statistical descriptors—such as skewness and kurtosis—to characterize the geometric features of the mRNA copy-number distribution. These quantities provide tractable and informative measures of transcriptional variability when explicit distributional forms are unavailable.

2.2. Moments, Skewness, and Kurtosis of Transcript Distributions

Let ${\mu }_k\left(t\right)$ denote the kth moment of the mRNA transcript number at time t. By definition,

$$\begin{array}{c}\hfill {\mu }_k\left(t\right)=\mathbb{E}\left[M^k\left(t\right)\right]=\sum _{m=0}^{\infty }{m}^{k}P(m,t),\end{array}$$

(2.3)

where $P(m,t)$ is the marginal probability mass function of the mRNA copy number.

The first moment ${\mu }_1\left(t\right)$ corresponds to the mean mRNA copy number and characterizes the average transcription level across the cell population. The second central moment defines the variance,

$$\begin{array}{c}\hfill {\sigma }^2\left(t\right)=\mathbb{E}\left[{(M\left(t\right)-{\mu }_1\left(t\right))}^2\right]={\mu }_2\left(t\right)-{\mu }_1^2\left(t\right),\end{array}$$

(2.4)

which quantifies the magnitude of transcriptional fluctuations around the mean.

To further characterize the geometric properties of the transcript number distribution, we consider higher-order standardized moments. The skewness is defined as

$$\begin{array}{c}\hfill \mathrm{Skew}\left(t\right)={\displaystyle \frac{\mathbb{E}\left[{(M\left(t\right)-{\mu }_1\left(t\right))}^3\right]}{{\sigma }^3\left(t\right)}},\end{array}$$

(2.5)

which measures the degree of asymmetry of the distribution. The excess kurtosis is defined by

$$\begin{array}{c}\hfill \mathrm{Kurt}\left(t\right)={\displaystyle \frac{\mathbb{E}\left[{(M\left(t\right)-{\mu }_1\left(t\right))}^4\right]}{{\sigma }^4\left(t\right)}}-3,\end{array}$$

(2.6)

and quantifies the heaviness of the distribution tails relative to a Gaussian distribution.

Equations (2.5) and (2.6) show that both skewness and kurtosis are fully determined by the first four moments of the transcript number distribution. Consequently, analyzing the dynamics of these moments provides direct insight into the temporal evolution of the shape of mRNA distributions arising from stochastic gene transcription.

2.3. Moment Expressions and Higher-Order Statistics

We now derive explicit expressions for the moments of the transcript number distribution. In the following analysis, we count only newly synthesized mRNA molecules and assume that the promoter initially resides in the gene OFF state at time $t=0$. Accordingly, the initial condition is

$$\begin{array}{c}\hfill M\left(0\right)=0,I\left(0\right)=0.\end{array}$$

(2.7)

Condition (2.7) implies

$$\begin{array}{c}\hfill P_0(0,0)=1,P_0(m,0)=0\mathrm{for}m\ge 1,P_1(m,0)=0\mathrm{for}m\ge 0.\end{array}$$

(2.8)

Consequently, all moments vanish at the initial time,

$$\begin{array}{c}\hfill {\mu }_k\left(0\right)=0,k\ge 1.\end{array}$$

(2.9)

Using calculations analogous to those in previous studies, the time-dependent mean mRNA copy number is obtained as

$$\begin{array}{c}\hfill {\mu }_1\left(t\right)={\displaystyle \frac{\nu \lambda }{\delta (\lambda +\gamma )}}-{\displaystyle \frac{\nu \lambda }{\delta (\lambda +\gamma -\delta )}}{e}^{-\delta t}+{\displaystyle \frac{\nu \lambda }{(\lambda +\gamma )(\lambda +\gamma -\delta )}}{e}^{-(\lambda +\gamma )t}.\end{array}$$

(2.10)

The last two terms in (2.10) decay to zero exponentially as $t\to \infty$. Hence, as $t\to \infty$, the mean mRNA copy number converges to the steady-state value

$${\mu }_1={\displaystyle \frac{\nu \lambda }{\delta (\lambda +\gamma )}}.$$

Similarly, the second moment of the transcript number distribution at steady state is given by

$$\begin{array}{c}\hfill {\mu }_2={\mu }_1\left[1+{\displaystyle \frac{\nu (\delta +\lambda )}{\delta (\delta +\lambda +\gamma )}}\right].\end{array}$$

(2.11)

The variance of the transcript number distribution then follows as

$$\begin{array}{c}\hfill {\sigma }^2={\mu }_2-{\mu }_1^2={\mu }_1\left[1+{\displaystyle \frac{\nu \gamma }{(\lambda +\gamma )(\delta +\lambda +\gamma )}}\right],\end{array}$$

(2.12)

where $\sigma$ denotes the standard deviation.

The squared coefficient of variation (noise) and the Fano factor are therefore given by

$$\begin{array}{c}\hfill {\mathrm{CV}}^2={\displaystyle \frac{{\sigma }^2}{{\mu }_1^2}}={\displaystyle \frac{1}{{\mu }_1}}+{\displaystyle \frac{\delta \gamma }{\lambda (\delta +\lambda +\gamma )}},\Phi ={\displaystyle \frac{{\sigma }^2}{{\mu }_1}}=1+{\displaystyle \frac{\nu \gamma }{(\lambda +\gamma )(\delta +\lambda +\gamma )}}.\end{array}$$

(2.13)

These quantities characterize the dispersion of transcript numbers across single cells but do not fully describe the shape of the underlying distribution.

To further elucidate the distributional geometry of mRNA copy numbers, we consider the skewness and kurtosis. By definition,

$$\begin{array}{cc}\hfill \mathrm{Skew}& ={\displaystyle \frac{{\mu }_3-3{\mu }_1{\mu }_2+2{\mu }_1^3}{{\sigma }^3}},\hfill \end{array}$$

(2.14)

$$\begin{array}{cc}\hfill \mathrm{Kurt}& ={\displaystyle \frac{{\mu }_4-4{\mu }_1{\mu }_3+6{\mu }_1^2{\mu }_2-3{\mu }_1^4}{{\sigma }^4}}-3,\hfill \end{array}$$

(2.15)

where ${\mu }_3$ and ${\mu }_4$ denote the third and fourth moments of the transcript number distribution, respectively. Explicit expressions for these higher-order moments are given by

$$\begin{array}{cc}\hfill {\mu }_3& ={\mu }_1\left[1+{\displaystyle \frac{3\nu (\delta +\lambda )}{\delta (\delta +\lambda +\gamma )}}+{\displaystyle \frac{{\nu }^2(\delta +\lambda )(2\delta +\lambda )}{{\delta }^2(\delta +\lambda +\gamma )(2\delta +\lambda +\gamma )}}\right],\hfill \\ \hfill {\mu }_4& ={\mu }_1\left[1+{\displaystyle \frac{7\nu (\delta +\lambda )}{\delta (\delta +\lambda +\gamma )}}+{\displaystyle \frac{6{\nu }^2(\delta +\lambda )(2\delta +\lambda )}{{\delta }^2(\delta +\lambda +\gamma )(2\delta +\lambda +\gamma )}}\right.\hfill \\ \hfill & \left.+{\displaystyle \frac{{\nu }^3(\delta +\lambda )(2\delta +\lambda )(3\delta +\lambda )}{{\delta }^3(\delta +\lambda +\gamma )(2\delta +\lambda +\gamma )(3\delta +\lambda +\gamma )}}\right].\hfill \end{array}$$

These expressions explicitly reveal how skewness and kurtosis depend on transcriptional kinetics, promoter switching dynamics, and mRNA degradation rates.

Biological interpretation of skewness and kurtosis. Skewness and kurtosis provide biologically meaningful information about the shape of the mRNA copy number distribution beyond what is captured by the mean and variance. A positive value of $\mathrm{Skew}$ indicates a pronounced right tail, corresponding to rare cells exhibiting exceptionally high transcript abundance. Such asymmetry is a characteristic feature of stochastic gene expression and is commonly associated with transcriptional bursting.

Kurtosis quantifies the weight of the distribution tails relative to a Gaussian distribution. A large positive value of $\mathrm{Kurt}$ implies an elevated probability of extreme transcriptional events, reflecting strong cell-to-cell heterogeneity. Together, skewness and kurtosis offer intuitive and quantitative measures of how transcriptional noise manifests at the population level and how strongly the mRNA distribution deviates from symmetric, light-tailed behavior.

3. Discussion

The explicit moment expressions derived above enable a systematic asymptotic analysis of the skewness and kurtosis in biologically relevant limits. In this section, we summarize the leading-order behavior of these higher-order statistics in three canonical parameter regimes that are commonly encountered in stochastic gene expression models. All kinetic parameters are nondimensionalized by the mRNA degradation rate $\delta$. Under this rescaling, both the numerical values and qualitative features of the moments remain unchanged. Therefore, without loss of generality, we set $\delta =1$ in the following analysis.

The asymptotic structure of the skewness and kurtosis is primarily governed by the relative timescales of promoter switching and mRNA degradation. Depending on whether promoter switching is fast or slow compared to transcript turnover, qualitatively distinct transcriptional behaviors emerge, each associated with characteristic distributional signatures [19,20,21,22]. Accordingly, we classify the asymptotic regimes into fast promoter switching, slow promoter switching, and a bursty transcription limit. For each regime, we derive the leading-order asymptotic expressions for skewness and kurtosis and interpret their implications for transcriptional variability and distributional shape.

Fast promoter switching regime.      Assume

$$\lambda +\gamma \gg 1,$$

so that promoter switching is much faster than mRNA degradation [23]. In this limit, the promoter rapidly equilibrates on the timescale of mRNA dynamics, and the system effectively reduces to a single birth–death process with mean

$${\mu }_1={\displaystyle \frac{\nu \lambda }{\lambda +\gamma }}.$$

Expanding the higher-order moments in powers of ${(\lambda +\gamma )}^{-1}$ yields, to leading order,

$$\begin{array}{cc}\hfill {\mu }_2& \approx {\mu }_1^2+{\mu }_1,\hfill \\ \hfill {\mu }_3& \approx {\mu }_1^3+3{\mu }_1^2+{\mu }_1,\hfill \\ \hfill {\mu }_4& \approx {\mu }_1^4+6{\mu }_1^3+7{\mu }_1^2+{\mu }_1.\hfill \end{array}$$

Consequently, the skewness and excess kurtosis scale asymptotically as

$$\begin{array}{cc}\hfill \mathrm{Skew}& \sim {\displaystyle \frac{1}{\sqrt{{\mu }_1}}},\hfill \\ \hfill \mathrm{Kurt}& \sim {\displaystyle \frac{1}{{\mu }_1}}.\hfill \end{array}$$

Thus, both skewness and excess kurtosis decrease as the mean expression level increases, and the mRNA copy number distribution approaches a Poisson-like form with weak asymmetry and light tails.

Slow promoter switching regime.      Assume

$$\lambda +\gamma \ll 1,$$

so that mRNA molecules degrade rapidly relative to promoter switching [24,25]. In this regime, transcription occurs in well-separated episodes, leading to pronounced cell-to-cell heterogeneity. Retaining the leading-order terms in $(\lambda +\gamma )$, we obtain the asymptotic expressions

$$\begin{array}{cc}\hfill {\mu }_2& \sim {\mu }_1^2\left(1+{\displaystyle \frac{\gamma }{\lambda }}\right),\hfill \\ \hfill {\mu }_3& \sim {\mu }_1^3\left(1+3{\displaystyle \frac{\gamma }{\lambda }}+2{\displaystyle \frac{{\gamma }^2}{{\lambda }^2}}\right),\hfill \\ \hfill {\mu }_4& \sim {\mu }_1^4\left(1+6{\displaystyle \frac{\gamma }{\lambda }}+11{\displaystyle \frac{{\gamma }^2}{{\lambda }^2}}+6{\displaystyle \frac{{\gamma }^3}{{\lambda }^3}}\right).\hfill \end{array}$$

Consequently, the skewness and excess kurtosis scale as

$$\begin{array}{cc}\hfill \mathrm{Skew}& \sim 2\sqrt{{\displaystyle \frac{\gamma }{\lambda }}},\hfill \\ \hfill \mathrm{Kurt}& \sim 6{\displaystyle \frac{\gamma }{\lambda }}.\hfill \end{array}$$

Both quantities increase with the ratio $\gamma /\lambda$, reflecting strong distributional asymmetry and heavy-tailed behavior driven by long inactive periods and rare promoter activation events. Although the distribution may exhibit bimodality in this regime, the large excess kurtosis originates from intermittent transcription bursts and long OFF periods, which generate strong deviations from Poisson-like statistics.

Bursty transcription regime.      Assume

$$\lambda \ll \gamma ,\nu \gg 1,$$

while keeping the burst size

$$B={\displaystyle \frac{\nu }{\gamma }}$$

finite. In this limit, transcription occurs in rare but intense bursts. The mean mRNA copy number is

$${\mu }_1=\lambda B.$$

The higher-order moments are asymptotically given by

$$\begin{array}{cc}\hfill {\mu }_2& \sim {\mu }_1^2(1+B),\hfill \\ \hfill {\mu }_3& \sim {\mu }_1^3(1+3B+2B^2),\hfill \\ \hfill {\mu }_4& \sim {\mu }_1^4(1+6B+11B^2+6B^3).\hfill \end{array}$$

From these, the skewness and excess kurtosis scale as

$$\begin{array}{cc}\hfill \mathrm{Skew}& \sim {\displaystyle \frac{2+B}{\sqrt{{\mu }_1}}},\hfill \\ \hfill \mathrm{Kurt}& \sim {\displaystyle \frac{6+6B+B^2}{{\mu }_1}}.\hfill \end{array}$$

These expressions show that increasing the burst size amplifies both skewness and excess kurtosis, even when the mean expression level is fixed, providing a natural mechanism for generating highly asymmetric, heavy-tailed mRNA distributions.

The asymptotic behavior of the mean, variance, skewness, and kurtosis in different promoter switching and transcription regimes is summarized in Table 1.

Table 1. Asymptotic expressions for steady-state mRNA statistics in different promoter switching and transcription regimes, including the squared coefficient of variation (CV²).

Regime	Mean ${\mathbf{\mu }}_{\mathbf{1}}$	Variance ${\mathbf{\sigma }}^{\mathbf{2}}$	CV2	Skewness	Kurtosis
Fast switching	${\displaystyle \frac{\nu \lambda }{\lambda +\gamma }}$	${\mu }_1$	$\sim {\mu }_1^{-1}$	$\sim {\mu }_1^{-1/2}$	$\sim {\mu }_1^{-1}$
Slow switching	${\displaystyle \frac{\nu \lambda }{\lambda +\gamma }}$	$\sim {\mu }_1^2\left(1+{\displaystyle \frac{\gamma }{\lambda }}\right)$	$\sim 1+{\displaystyle \frac{\gamma }{\lambda }}$	$\sim 2\sqrt{{\displaystyle \frac{\gamma }{\lambda }}}$	$\sim 6{\displaystyle \frac{\gamma }{\lambda }}$
Bursty transcription	$\sim \lambda B$	$\sim {\mu }_1(1+B)$	$\sim {\displaystyle \frac{1+B}{{\mu }_1}}$	$\sim {\displaystyle \frac{2+B}{\sqrt{{\mu }_1}}}$	$\sim {\displaystyle \frac{6+6B+B^2}{{\mu }_1}}$

Table 1. Asymptotic expressions for steady-state mRNA statistics in different promoter switching and transcription regimes, including the squared coefficient of variation (CV²).

Regime	Mean ${\mathbf{\mu }}_{\mathbf{1}}$	Variance ${\mathbf{\sigma }}^{\mathbf{2}}$	CV2	Skewness	Kurtosis
Fast switching	${\displaystyle \frac{\nu \lambda }{\lambda +\gamma }}$	${\mu }_1$	$\sim {\mu }_1^{-1}$	$\sim {\mu }_1^{-1/2}$	$\sim {\mu }_1^{-1}$
Slow switching	${\displaystyle \frac{\nu \lambda }{\lambda +\gamma }}$	$\sim {\mu }_1^2\left(1+{\displaystyle \frac{\gamma }{\lambda }}\right)$	$\sim 1+{\displaystyle \frac{\gamma }{\lambda }}$	$\sim 2\sqrt{{\displaystyle \frac{\gamma }{\lambda }}}$	$\sim 6{\displaystyle \frac{\gamma }{\lambda }}$
Bursty transcription	$\sim \lambda B$	$\sim {\mu }_1(1+B)$	$\sim {\displaystyle \frac{1+B}{{\mu }_1}}$	$\sim {\displaystyle \frac{2+B}{\sqrt{{\mu }_1}}}$	$\sim {\displaystyle \frac{6+6B+B^2}{{\mu }_1}}$

Connection to negative binomial distributions. In bursty transcription limits, the steady-state mRNA copy number distribution converges to a negative binomial distribution with mean ${\mu }_1$ and variance

$${\sigma }^2={\mu }_1(1+B),$$

where B denotes the burst size. The negative binomial probability mass function (PMF) can be written as

$$\begin{array}{c}\hfill P\left(m\right)=\left({\displaystyle \genfrac{}{}{0pt}{}{m+r-1}{m}}\right)p^r{(1-p)}^m,\end{array}$$

(18)

with parameters

$$\begin{array}{c}\hfill r={\displaystyle \frac{{\mu }_1^2}{{\sigma }^2-{\mu }_1}}={\displaystyle \frac{{\mu }_1}{B}},p={\displaystyle \frac{{\mu }_1}{{\sigma }^2}}={\displaystyle \frac{1}{1+B}}.\end{array}$$

(19)

This negative binomial distribution reproduces the asymptotic expressions for skewness and excess kurtosis derived above:

$$\begin{array}{c}\hfill \mathrm{Skew}={\displaystyle \frac{2+B}{\sqrt{{\mu }_1}}},\mathrm{Kurt}={\displaystyle \frac{6+6B+B^2}{{\mu }_1}}.\end{array}$$

(20)

This explicit correspondence provides a theoretical justification for the widespread use of negative binomial models in fitting single-cell RNA-sequencing data, where transcriptional bursting and promoter switching naturally generate transcript count distributions that are overdispersed relative to a Poisson process. The relative timescale of promoter switching and mRNA degradation critically determines transcriptional variability and distributional shape. In the fast-switching regime, the promoter rapidly equilibrates between active and inactive states, so that transcription effectively proceeds at an averaged rate [20,26]. Consequently, transcriptional bursts are weak and the mRNA copy number distribution approaches a Poisson-like form with reduced skewness and kurtosis. In contrast, slow promoter switching gives rise to long inactive periods punctuated by rare activation events that generate intense transcriptional bursts [7,14,27]. This temporal clustering leads to highly asymmetric and heavy-tailed mRNA distributions, characterized by large positive skewness and elevated kurtosis. Importantly, these higher-order features persist even when mean expression levels are fixed, demonstrating that slow promoter switching encodes mechanistic information beyond that captured by mean–variance statistics alone.

4. Numerical Simulations

In this section, all parameters are nondimensionalized by the mRNA degradation rate. We define

$$\nu =200{\min }^{-1},\lambda =0.5{\min }^{-1},\gamma =0.5{\min }^{-1},\delta =1{\min }^{-1}.$$

This normalization does not affect the values or qualitative properties of the mean, skewness, or kurtosis of the mRNA distribution.

To identify the dominant regulatory factors underlying the skewness and kurtosis of mRNA copy number distributions, we examined their joint dependence on the promoter activation rate $\lambda$ and the mRNA degradation rate $\gamma$ using three-dimensional representations. In Figure 2, the three-dimensional surfaces display the values of the skewness (A) and kurtosis (B) of the steady-state mRNA distribution as functions of $\lambda$ and $\gamma$, both varying within the range $[0.05,5]$. Red dots mark the locations of the minimum values on each surface. Slow promoter switching, corresponding to small $\lambda$ or $\gamma$, gives rise to large skewness and kurtosis due to pronounced transcriptional bursting, whereas faster switching reduces higher-order variability. Interestingly, when $\lambda$ and $\gamma$ are comparable in magnitude and both are small, a global minimum of both skewness and kurtosis is observed. These surfaces provide a clear visualization of how promoter kinetics shape mRNA variability across a biologically relevant parameter space.

To further assess the impact of promoter activation dynamics, we compare two representative regimes: a slow-switching promoter ($\lambda =0.5$) and a fast-switching promoter ($\lambda =5$), while keeping all other parameters fixed. Figure 3 shows the steady-state distributions of mRNA copy numbers obtained from exact Gillespie simulations of the two-state gene transcription model.

For slow promoter activation ($\lambda =0.5$, Figure 3A), the gene spends long periods in the inactive state, interrupted by rare but extended activation events. During each active period, mRNA molecules are produced rapidly, giving rise to pronounced transcriptional bursts. As a consequence, the distribution is strongly right-skewed. Moreover, the coexistence of low- and high-expression subpopulations leads to a bimodal structure, in which probability mass is distributed between two modes rather than concentrated in extreme tails, resulting in a reduced kurtosis compared to burst-dominated, heavy-tailed distributions.

In contrast, for fast promoter activation ($\lambda =5$, Figure 3B), frequent switching between the active and inactive states effectively averages transcriptional activity over time. This suppresses burstiness and yields a narrower and more symmetric mRNA distribution with reduced skewness and kurtosis.

In both regimes, the vertical lines denote the mean mRNA copy number. Analytical predictions derived from the exact expressions of the first four moments are in quantitative agreement with the simulation results, validating the theoretical analysis of mean expression levels and higher-order fluctuations.

5. Conclusions

Stochastic gene expression is most often characterized using low-order statistical descriptors such as the mean expression level, variance, or related noise measures including the Fano factor. While these quantities effectively quantify the overall magnitude of transcriptional variability, it has long been recognized that they are insufficient to capture the full geometric structure of mRNA copy number distributions generated by promoter switching and transcriptional bursting [3,11]. In this work, we demonstrate that skewness provides essential information about distributional asymmetry that cannot be inferred from variance alone.

Positive skewness emerges naturally from bursty transcription dynamics, in which rare promoter activation events give rise to intense transcriptional episodes [7,14,27]. Our analytical results show that the skewness $\mathrm{Skew}$ depends explicitly on promoter switching kinetics and burst size, increasing with the ratio $\gamma /\lambda$ and with the burst size $B=\nu /\gamma$, even when the mean expression level is held fixed. This establishes skewness as a mechanistically informative statistic that reflects the underlying regulatory architecture rather than merely the overall noise magnitude. From a biological perspective, skewness quantifies the contribution of rare, high-expression cells within a population. Such cells can play a disproportionate role in downstream regulatory processes, including threshold-dependent responses, stress adaptation, and phenotypic diversification [25,28]. Variance-based measures treat fluctuations symmetrically around the mean and therefore obscure the statistical weight of these rare transcriptional events, whereas skewness directly captures the asymmetry induced by transcriptional bursting and slow promoter switching.

Beyond skewness, kurtosis further characterizes the sharpness and tail heaviness of mRNA distributions. Our results indicate that low mRNA abundance regimes are associated with elevated skewness and kurtosis, reflecting highly asymmetric and heavy-tailed distributions dominated by rare transcriptional bursts. As mRNA abundance increases, both skewness and kurtosis decrease, and the distribution approaches a Gaussian-like form. This behavior is consistent with effective averaging over many transcription events and aligns with previous theoretical analyses of noise–mean relationships in gene expression models [8,29].

Finally, our analysis provides a theoretical foundation for extending negative binomial–based analyses of single-cell transcriptomic data beyond conventional mean–variance fitting [6]. Distinct promoter dynamics may generate identical mean expression levels and variances while exhibiting markedly different skewness and kurtosis, a limitation that has been highlighted in recent studies on model identifiability and model reduction [26,30]. Incorporating higher-order statistics therefore offers a principled route toward improved discrimination between competing transcriptional mechanisms and establishes a tighter quantitative link between statistical inference and the underlying molecular regulation of gene expression.

Limitations and future directions. The present study focuses on a minimal two-state promoter model with constant kinetic parameters, which, while analytically tractable, represents an idealized description of transcriptional regulation. In particular, real promoters often exhibit multiple inactive or active states, transcription factor–dependent switching rates, and memory effects that can give rise to non-exponential waiting times between transcriptional events. Extending the analytical framework developed here to multi-state promoter architectures or non-Markovian switching processes would provide a more comprehensive characterization of transcriptional heterogeneity.

Moreover, our analysis is restricted to intrinsic noise arising from stochastic transcription and promoter switching and does not explicitly account for extrinsic sources of variability such as cell-to-cell differences in transcription factor abundance, cell cycle stage, or global physiological state. Incorporating extrinsic fluctuations, either through additional stochastic variables or through hierarchical modeling approaches, represents an important direction for future work.

Finally, while the present results establish skewness and kurtosis as mechanistically informative statistics at the mRNA level, an important open question concerns how these higher-order features propagate through translation and post-transcriptional regulation to shape protein-level variability. Systematic analysis of higher-order moments in coupled mRNA–protein models, together with direct comparisons to single-cell transcriptomic and proteomic data, will be essential for further bridging stochastic gene expression theory with experimental inference.

Appendix A.1. Promoter State Probabilities

Summing over all transcript numbers m, define

$$P_0\left(t\right)=\sum _{m=0}^{\infty }{P}_0(m,t),P_1\left(t\right)=\sum _{m=0}^{\infty }{P}_1(m,t),$$

which give the probability that the promoter resides in the OFF and ON states, respectively. Since these form a complete set of promoter states,

$$P_0\left(t\right)+P_1\left(t\right)=1,\forall t\ge 0.$$

Differentiating and using the CME, we obtain

$$P_0^{\prime }\left(t\right)=\gamma P_1\left(t\right)-\lambda P_0\left(t\right),P_1^{\prime }\left(t\right)=\lambda P_0\left(t\right)-\gamma P_1\left(t\right).$$

Solving with initial conditions $P_0\left(0\right)=1$, $P_1\left(0\right)=0$, we get

$$\begin{array}{cc}\hfill P_0\left(t\right)& ={\displaystyle \frac{\gamma }{\lambda +\gamma }}+{\displaystyle \frac{\lambda }{\lambda +\gamma }}{e}^{-(\lambda +\gamma )t},\hfill \\ \hfill P_1\left(t\right)& ={\displaystyle \frac{\lambda }{\lambda +\gamma }}-{\displaystyle \frac{\lambda }{\lambda +\gamma }}{e}^{-(\lambda +\gamma )t}.\hfill \end{array}$$

The steady-state probabilities are

$$P_0={\displaystyle \frac{\gamma }{\lambda +\gamma }},P_1={\displaystyle \frac{\lambda }{\lambda +\gamma }}.$$

Appendix A.2. First Moment (Mean)

Define the k-th moment of transcript numbers as

$${\mu }_k\left(t\right)=\sum _{m=0}^{\infty }{m}^{k}P(m,t),\mathrm{with}P(m,t)=P_0(m,t)+P_1(m,t).$$

Decomposing by promoter state,

$${\mu }_k\left(t\right)={\mu }_{k0}\left(t\right)+{\mu }_{k1}\left(t\right),{\mu }_{k0}\left(t\right)=\sum _{m=0}^{\infty }{m}^k{P}_0(m,t),{\mu }_{k1}\left(t\right)=\sum _{m=0}^{\infty }{m}^k{P}_1(m,t).$$

For $k=1$ (the mean),

$$\begin{array}{cc}\hfill {\mu }_{10}^{\prime }\left(t\right)& =\gamma {\mu }_{11}\left(t\right)-(\lambda +\delta ){\mu }_{10}\left(t\right),\hfill \end{array}$$

(A1)

$$\begin{array}{cc}\hfill {\mu }_{11}^{\prime }\left(t\right)& =\nu P_1\left(t\right)+\lambda {\mu }_{10}\left(t\right)-(\gamma +\delta ){\mu }_{11}\left(t\right).\hfill \end{array}$$

(A2)

Step 1: Steady-state solution.

Setting derivatives to zero, the equilibrium moments are

$${\mu }_{10}={\displaystyle \frac{\nu \lambda \gamma }{\delta (\lambda +\gamma )(\delta +\lambda +\gamma )}},{\mu }_{11}={\displaystyle \frac{\nu \lambda (\delta +\lambda )}{\delta (\lambda +\gamma )(\delta +\lambda +\gamma )}}.$$

The total mean is

$${\mu }_1={\mu }_{10}+{\mu }_{11}={\displaystyle \frac{\nu \lambda }{\delta (\lambda +\gamma )}}.$$

Step 2: Time-dependent solution.

Solving (A1)–(A2) with initial condition ${\mu }_{10}\left(0\right)={\mu }_{11}\left(0\right)=0$ yields

$$\begin{array}{cc}\hfill {\mu }_{10}\left(t\right)& ={\mu }_{10}-A_1{e}^{-\delta t}-A_2{e}^{-(\lambda +\gamma )t}-A_3{e}^{-(\delta +\lambda +\gamma )t},\hfill \\ \hfill {\mu }_{11}\left(t\right)& ={\mu }_{11}-B_1{e}^{-\delta t}-B_2{e}^{-(\lambda +\gamma )t}+B_3{e}^{-(\delta +\lambda +\gamma )t},\hfill \end{array}$$

with constant coefficients $A_i,B_i$. The sum gives the time-dependent mean ${\mu }_1\left(t\right)$.

Appendix A.3. Second Moment (Variance)

Similarly, for $k=2$,

$$\begin{array}{cc}\hfill {\mu }_{20}^{\prime }\left(t\right)& =\delta {\mu }_{10}\left(t\right)+\gamma {\mu }_{21}\left(t\right)-(2\delta +\lambda ){\mu }_{20}\left(t\right),\hfill \\ \hfill {\mu }_{21}^{\prime }\left(t\right)& =\nu P_1\left(t\right)+2\nu {\mu }_{11}\left(t\right)+\delta {\mu }_{11}\left(t\right)+\lambda {\mu }_{20}\left(t\right)-(2\delta +\gamma ){\mu }_{21}\left(t\right).\hfill \end{array}$$

Step 1: Steady-state solution.

Setting derivatives to zero,

$$\begin{array}{cc}\hfill {\mu }_{20}& ={\displaystyle \frac{\nu \lambda \gamma \left[\nu \right(\delta +\lambda )+\delta (2\delta +\lambda +\gamma \left)\right]}{{\delta }^2(\lambda +\gamma )(\delta +\lambda +\gamma )(2\delta +\lambda +\gamma )}},\hfill \\ \hfill {\mu }_{21}& ={\displaystyle \frac{\nu \lambda (\delta +\lambda )\left[\nu \right(2\delta +\lambda )+\delta (2\delta +\lambda +\gamma \left)\right]}{{\delta }^2(\lambda +\gamma )(\delta +\lambda +\gamma )(2\delta +\lambda +\gamma )}}.\hfill \end{array}$$

The total second moment is ${\mu }_2={\mu }_{20}+{\mu }_{21}$, and the variance is

$${\sigma }^2={\mu }_2-{\mu }_1^2.$$

Appendix A.4. Third and Fourth Moments

Third moment:

$$\begin{array}{cc}\hfill {\mu }_{30}^{\prime }\left(t\right)& =\gamma {\mu }_{31}\left(t\right)-\delta {\mu }_{10}\left(t\right)+3\delta {\mu }_{20}\left(t\right)-(3\delta +\lambda ){\mu }_{30}\left(t\right),\hfill \\ \hfill {\mu }_{31}^{\prime }\left(t\right)& =\nu [P_1\left(t\right)+3{\mu }_{11}\left(t\right)+3{\mu }_{21}\left(t\right)]+\lambda {\mu }_{30}\left(t\right)-\delta {\mu }_{11}\left(t\right)+3\delta {\mu }_{21}\left(t\right)-(3\delta +\gamma ){\mu }_{31}\left(t\right).\hfill \end{array}$$

At steady state,

$${\mu }_3={\mu }_{30}+{\mu }_{31}.$$

Fourth moment:

$${\mu }_4^{\prime }\left(t\right)=\nu \left[P_1\left(t\right)+4{\mu }_{11}\left(t\right)+6{\mu }_{21}\left(t\right)+4{\mu }_{31}\left(t\right)\right]+\delta \left[{\mu }_1\left(t\right)-4{\mu }_2\left(t\right)+6{\mu }_3\left(t\right)-4{\mu }_4\left(t\right)\right].$$

Steady-state solution ${\mu }_4$ is obtained by setting derivative to zero and solving algebraically.

Appendix B.1. Fast-Switching Regime

In the fast-switching limit, the promoter switching rates are much larger than mRNA degradation,

$$\lambda +\gamma \gg 1.$$

The exact moments contain terms proportional to ${(\lambda +\gamma )}^{-1}$ and transient exponential terms $e^{-(\lambda +\gamma )t}$. Since these decay rapidly, we can neglect them at steady state. Expanding the exact expressions in powers of ${(\lambda +\gamma )}^{-1}$, we obtain

$$\begin{array}{cc}\hfill {\mu }_1& ={\displaystyle \frac{\nu \lambda }{\lambda +\gamma }},\hfill \\ \hfill {\mu }_2& ={\mu }_1^2+{\mu }_1+\mathcal{O}\left({(\lambda +\gamma )}^{-1}\right),\hfill \\ \hfill {\mu }_3& ={\mu }_1^3+3{\mu }_1^2+{\mu }_1+\mathcal{O}\left({(\lambda +\gamma )}^{-1}\right),\hfill \\ \hfill {\mu }_4& ={\mu }_1^4+6{\mu }_1^3+7{\mu }_1^2+{\mu }_1+\mathcal{O}\left({(\lambda +\gamma )}^{-1}\right).\hfill \end{array}$$

From the definitions of skewness and excess kurtosis,

$$\mathrm{Skew}={\displaystyle \frac{{\mu }_3-3{\mu }_1{\mu }_2+2{\mu }_1^3}{{\sigma }^3}},\mathrm{Kurt}={\displaystyle \frac{{\mu }_4-4{\mu }_1{\mu }_3+6{\mu }_1^2{\mu }_2-3{\mu }_1^4}{{\sigma }^4}},$$

we obtain to leading order

$$\mathrm{Skew}\sim {\displaystyle \frac{1}{\sqrt{{\mu }_1}}},\mathrm{Kurt}\sim {\displaystyle \frac{1}{{\mu }_1}}.$$

Thus, in the fast-switching regime, the distribution is close to Poisson, with weak asymmetry and light tails.

Appendix B.2. Slow-Switching Regime

In the slow-switching limit, mRNA degradation is fast relative to promoter switching,

$$\lambda +\gamma \ll 1.$$

In the exact expressions for ${\mu }_2,{\mu }_3,{\mu }_4$, the terms $1+\lambda +\gamma \approx 1$ and $2+\lambda +\gamma \approx 2$, etc. Retaining leading-order terms in $\lambda +\gamma$ yields

$$\begin{array}{cc}\hfill {\mu }_2& \sim {\mu }_1^2\left(1+{\displaystyle \frac{\gamma }{\lambda }}\right),\hfill \\ \hfill {\mu }_3& \sim {\mu }_1^3\left(1+3{\displaystyle \frac{\gamma }{\lambda }}+2{\displaystyle \frac{{\gamma }^2}{{\lambda }^2}}\right),\hfill \\ \hfill {\mu }_4& \sim {\mu }_1^4\left(1+6{\displaystyle \frac{\gamma }{\lambda }}+11{\displaystyle \frac{{\gamma }^2}{{\lambda }^2}}+6{\displaystyle \frac{{\gamma }^3}{{\lambda }^3}}\right).\hfill \end{array}$$

Plugging into the skewness and excess kurtosis formulas gives

$$\mathrm{Skew}\sim 2\sqrt{{\displaystyle \frac{\gamma }{\lambda }}},\mathrm{Kurt}\sim 6{\displaystyle \frac{\gamma }{\lambda }}.$$

Derivation step:

1. Start from the exact moments (see Appendix A), e.g.,

$${\mu }_2={\mu }_1\left[1+{\displaystyle \frac{\nu (1+\lambda )}{1+\lambda +\gamma }}\right]\approx {\mu }_1^2\left(1+{\displaystyle \frac{\gamma }{\lambda }}\right),$$

using $1\gg \lambda +\gamma$. 2. Similarly, expand ${\mu }_3$ and ${\mu }_4$ in powers of $\lambda +\gamma$, keeping the leading terms. 3. Compute central moments to obtain skewness and excess kurtosis.

Appendix B.3. Bursty Transcription Regime

The bursty limit corresponds to rare but intense transcription events,

$$\lambda \ll \gamma ,\nu \gg \delta ,$$

while keeping the burst size

$$B={\displaystyle \frac{\nu }{\gamma }}=\mathrm{finite}.$$

Step 1: Mean and variance

In this limit, the mean and variance reduce to

$${\mu }_1={\displaystyle \frac{\lambda B}{\delta }},{\sigma }^2={\mu }_1(1+B),$$

which coincide with a negative binomial distribution.

Step 2: Higher-order moments

Expanding the exact expressions of ${\mu }_3,{\mu }_4$ in powers of B and keeping leading-order terms,

$$\begin{array}{cc}\hfill {\mu }_3& \sim {\mu }_1(1+3B+2B^2),\hfill \\ \hfill {\mu }_4& \sim {\mu }_1^2(1+6B+11B^2+6B^3),\hfill \end{array}$$

which after standardization yield

$$\mathrm{Skew}\sim {\displaystyle \frac{2+B}{\sqrt{{\mu }_1}}},\mathrm{Kurt}\sim {\displaystyle \frac{6+6B+B^2}{{\mu }_1}}.$$

Step 3: Connection to negative binomial

The negative binomial PMF with parameters

$$r={\displaystyle \frac{{\mu }_1^2}{{\sigma }^2-{\mu }_1}},p={\displaystyle \frac{{\mu }_1}{{\sigma }^2}}$$

reproduces the above skewness and excess kurtosis, confirming that bursty transcription naturally leads to overdispersed transcript count distributions.

Taken together, the slow-switching and bursty limits illustrate how promoter kinetics and burst size control the higher-order statistics of mRNA distributions, including skewness and excess kurtosis.