This is an uncorrected proof.
Figures
Abstract
Single-cell transcriptomic data exhibit pervasive zero inflation, while traditional models either neglect this issue or fail to capture transcriptional burst-driven bimodality, hindering accurate gene regulatory studies. This study developed a zero-inflated telegraph model that integrates technical zero correction with the stochastic gene state-switching dynamics of the classical telegraph model. Systematic validation was conducted using synthetic data, human scRNA-seq data from lupus and breast cancer patients, and mouse embryonic stem cell scRNA-seq data. The model showed superior performance: it accurately fits mRNA distributions (including bimodal patterns), reliably estimates effective transcriptional burst parameters while preventing overfitting, thus enables correction of traditional models’ regulatory inference bias. It also outperforms conventional approaches in detecting differentially expressed genes, with notable advantages in small samples, and identifies unique disease-related genes (e.g., LDLR, GZMB for lupus, FAIM2, VDR for breast cancer). This biologically interpretable and robust tool advances single-cell transcriptomic analysis.
Author summary
Single-cell studies demonstrate that gene expression is an inherently stochastic process driven by transcriptional bursting. Genes randomly toggle between active and inactive states, causing variable mRNA abundance even among genetically identical cells. The widespread absence of transcripts in such data arises from both biological noise and sequencing dropouts, a phenomenon referred to as zero inflation. Many existing approaches struggle to simultaneously address zero inflation and burst-associated expression dynamics. To address this challenge, we present a zero-inflated telegraph model that integrates technical zero correction with stochastic gene regulatory dynamics. Tested on synthetic data, mouse and human single-cell transcriptomic datasets, our model produces reliable mRNA distribution fits, distinguishes biological signals from technical artifacts, and yields reasonable estimates of transcriptional burst parameters. It also performs well in differential gene detection, especially with limited samples. This framework may serve as a useful tool for single-cell transcriptomic analysis.
Citation: Yang C, Liao Y, Sheng Y, Jiao F (2026) Integrating zero-inflation correction and transcriptional kinetics for single-cell transcriptomic analysis. PLoS Comput Biol 22(9): e1014779. https://doi.org/10.1371/journal.pcbi.1014779
Editor: Chun-Chun Wang, Jiangnan University, CHINA
Received: March 9, 2026; Accepted: August 31, 2026; Published: September 25, 2026
Copyright: © 2026 Yang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The transcriptomic scRNA-seq datasets from both human patients and mouse samples used in this study were obtained from the cited references [11,17,31,32]: PBMC datasets from lupus patients are available via https://figshare.com/ndownloader/files/34464122 TIL datasets from breast cancer patients are available via https://github.com/cz-ye/DECENT-analysis/tree/master/savas Smart-seq2 data can be accessed via the GitHub link https://github.com/sandberg-lab/txburst.git while RamDA-seq data were obtained from the GEO database with the accession ID GSE132589.
Funding: This work was supported by grants from the Natural Science Foundation of China (No.\,12271118 to F.J.), Guangzhou Municipal Science and Technology Plan Basic and Applied Basic Research (No.\,2025A03J3089 to F.J.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Gene transcription is an inherently stochastic process, a phenomenon evidenced by the temporal fluctuations in messenger RNA (mRNA) copy numbers observed even within clonally homogeneous (isogenic) cell populations [1,2]. Quantifying mRNA molecules at the single-cell level generates large-scale datasets of mRNA copy number distributions (), which explicitly quantify the probability that a single cell contains exactly m mRNA molecules of the gene of interest [3,4]. When integrated with mechanistic mathematical models, these datasets provide a statistical framework for estimating key system parameters and elucidating the underlying mechanisms of stochastic gene regulation. Such analyses hold significant utility across diverse biological fields, including inferring the genetic circuits governing cell fate determination [5,6] and cellular responses to environmental perturbations [7,8], as well as identifying differentially expressed genes (DEGs) between control and pathogenic conditions to prioritize therapeutic targets [9,10].
To model the stochastic nature of transcription, the traditional constitutive model is insufficient. This is because the Poisson distribution it generates cannot describe the observed distribution, where the mean deviates significantly from the distribution peak [1]. In contrast, this phenomenon can be well captured by the classical telegraph model. It posits that a gene stochastically transitions between an active state and an inactive state; mRNA molecules are synthesized during the active state [1,3,11]. The telegraph model can effectively characterize the unimodal and bimodal distributions observed in bacteria [2], yeast [12], and mammalian cells [11], establishing a basic model framework for linking gene regulatory mechanisms with distribution data. Building on the telegraph model, transcription models with different gene state-switching mechanisms have been derived, such as those containing multiple gene states, multiple signaling pathways, and positive/negative feedback loops [13–15].
The telegraph model is widely utilized to fit mRNA expression distributions in single-cell RNA-seq (scRNA-seq) data, owing to its biologically interpretable structure [11,16–18]. For instance, it accurately models 90% of transcripts in scRNA-seq datasets and outperforms empirical bulk-focused gamma-Poisson models in capturing bimodal distributions [9]-a hallmark that directly characterizes heterogeneity in isogenic cell populations [1,12]. The model infers two core parameters: transcriptional burst frequency (the rate of switching from inactive to active transcription), and burst size (the average mRNA copies produced per active state). Notably, and estimated by the telegraph model are regarded as effective parameters. Even for genes regulated by more complex regulatory mechanisms, these effective metrics remain widely applicable and can partly reflect genuine changes in transcriptional burst frequency and size [11,17–20].
Accurate estimation of and is critical for deciphering the intensity of gene regulation and its contextual dynamics. For example, in mouse embryonic stem (ES) cells, transcription factors Oct4 and Nanog regulate post-replication [21]; in Escherichia coli (E. coli), over 20 promoters modulate in response to growth signals [8]; and in human T cells, HIV LTR promoters co-regulate both parameters under distinct chromatin environments [7]. Recent studies have extended the telegraph model from analyzing specific gene circuits to genome-wide transcription distributions in mouse and human cells. These extensions revealed that promoter elements (e.g., TATA box and initiator) synergistically enhance [11,22], whereas specific regulators and signaling pathways predominantly target either or [17,18]. For instance, MYC-bound genes exhibit larger , while AFF4-targeted genes show enriched [22].
The telegraph model has also gained growing attention for genome-wide identification of DEGs across distinct cell types or states [9,10,16,19]. This progress relies on efficient algorithms to compute its distribution-for example, the Finite-State Projection (FSP) algorithm [3,20] and Gauss-Jacobi quadrature method [9]. Compared to traditional DEGs analysis methods, it offers two key advantages. First, it deciphers DEGs’ intrinsic regulatory changes from a biophysical perspective by linking them to burst frequency and size . For example [19], showed that across mouse CAST and C57 alleles (fibroblasts vs. embryonic stem cells): predominantly regulates cell cycle-related DEGs, while -regulated DEGs enrich in signaling and energy metabolism. Second, it overcomes a key limitation of traditional methods: failure to detect DEGs with unaltered mean but changed expression distributions, making it ideal for capturing heterogeneity-linked differences. For instance [16], found that in mouse brain tissue, the gene Ndnf (encoding neuron-derived neurotrophic factor) showed significant differences in inferred burst size between the L6 corticothalamic neuronal subclass and others, despite no detectable change in mRNA mean level.
However, neither the telegraph model nor its extensions address a critical issue: pervasive zero inflation that exists in single-cell transcriptional datasets. Zero inflation stems from two distinct sources: biological zeros (true non-expression, arising from stochastic gene state switching inherent to single-cell heterogeneity) and technical zeros (well-documented artifacts in scRNA-seq, i.e., dropout events like inefficient mRNA capture or amplification failure) [23,24]. This conflates the zero probability P0 in mRNA distributions, which incorrectly aggregates biological and technical zeros. Consequently, P0 is artificially inflated in most mRNA distribution data, severely undermining the telegraph model’s ability to capture the true distributional characteristics [25]. For instance, the bimodal distribution of paternal Mbnl2 transcription in mouse fibroblasts is misfitted as unimodal by the telegraph model using standard maximum-likelihood estimation method (MLE) [11]. In contrast, our newly developed fitting approach-by ensuring precise fitting of P0 during distribution data fitting-enables the telegraph model to accurately recapitulate bimodal patterns [26]. Notably, this method also enhances fitting performance for unimodal mRNA distributions in mammalian cells [3,26].
While zero inflation has been integrated into stochastic gene transcription models, existing approaches primarily incorporate a zero-inflation parameter into the negative binomial distribution [23,27], a strategy widely adopted for single-cell RNA-seq data analysis [28–30]. Such models have been demonstrated to identify three classes of DEGs during embryonic development, with each class enriched in distinct biological processes [23]. However, large-scale single-cell datasets have confirmed that stochastic transcriptional regulation inherently generates bimodal distributions, whereas the negative binomial distribution alone exclusively yields unimodal patterns. Thus, such models generate bimodality solely through technical zero-inflation, and therefore fail to capture the transcriptional burst regulation that underlies bimodality. Therefore, it is imperative to develop a telegraph model integrated with zero inflation to investigate transcriptional regulatory mechanisms and DEGs.
In this study, we developed a zero-inflated telegraph model that integrates technical zero correction with the classical telegraph model’s ability to capture stochastic gene state-switching dynamics. We systematically validated this model using diverse datasets: synthetic data with dropout, extrinsic noise and varied sample sizes; human scRNA-seq data from lupus patient CD14+ monocytes [31] and breast cancer patient CD8+ T cells [32]; and mouse scRNA-seq data from embryonic stem cells [11,17]. We demonstrated that the model outperforms conventional models in fitting mRNA distributions, reliably estimates transcriptional burst parameters without overfitting, and enhances DEGs detection and their subsequent functional enrichment. Collectively, this model provides a robust, biologically interpretable tool for single-cell transcriptomic analysis.
Results
Zero-inflated telegraph model
In the classical telegraph model [1,3,4,8], as illustrated in the following diagram:
the gene is proposed to transition stochastically between an active (on) state and an inactive (off) state, governed by first-order rate constants: the activation rate kon and inactivation rate koff. During the active state, mRNA molecules are synthesized at a first-order production rate kb. For mRNA molecules, their degradation, which includes both decay and dilution during cell division, occurs at a first-order decay rate kd. Then the burst frequency and size are expressed as [1,3]
The telegraph model allows for deriving the steady-state exact formulas of mRNA distribution through solving master equations [3,15]:
where denotes the confluent hypergeometric function [33]. Note that the exact formula of is expressed in terms of the ratios and . To improve clarity, we standardize the degradation rate kd by setting kd = 1 [1,20]. Accordingly, all model parameters should be interpreted as the ratio of their actual values to kd. When applying the telegraph model to the fitting of genome-wide distribution data, directly employing its exact expression for computations yields unsatisfactory precision and efficiency [3,10]. Instead, we can treat of the telegraph model as a beta-Poisson distribution [4,9]
(1)which can be efficiently computed via the Gauss-Jacobi quadrature method [9].
The telegraph model has a well-known special form corresponding to the negative binomial distribution. This limiting case arises when the inactivation rate koff and mRNA synthesis rate kb are both sufficiently large, while the burst frequency kon and size stay within regular ranges. Under this scenario, genes remain predominantly in the inactive state; upon activation, they rapidly produce abundant transcripts, and the resulting expression counts follow a negative binomial distribution [3]:
(2)It has been shown that a substantial portion of genome-wide distribution data in mouse and human cells can be described by the negative binomial distribution [11,19]. However, it can be mathematically proven that this distribution only generates unimodal distributions, and thus fails to capture the bimodal distributions that characterize heterogeneity [1,2,11].
To use either the telegraph model or negative binomial distribution to describe read counts while accounting for excess zero inflation, we introduce the average technical zero rate for each cell, where technical zero refers to zeros caused by technical failures in mRNA detection. Let N denote the total number of cells in an isogenic cell population, and denote the number of cells generating exactly m mRNA copies (true mRNA counts). Let represent the observed distribution that incorporates these technical zeros, while denotes the true distribution generated by inherent gene regulatory mechanisms. Then we have
(3)When is expressed via Eq. (1), denotes the exact mRNA distribution formula for the zero-inflated telegraph model; when expressed via Eq. (2), this formulation corresponds to the zero-inflated negative binomial (ZINB) model [23], which is widely adopted in modern single-cell transcriptomic analyses [28–30]. Notably, although both the zero-inflated telegraph and ZINB models can produce bimodal mRNA distributions, their underlying mechanisms differ substantially. The telegraph model intrinsically generates bimodal patterns from random transitions between the gene’s active and inactive states during transcriptional bursting. By contrast, ZINB achieves bimodality by combining the non-zero peak derived from the negative binomial distribution and an extra zero peak introduced by the technical zero rate . We proved that Eq. (3) ensures a unique parameter set exists for every observed distribution (Methods).
Zero-inflated telegraph model robustly captures mRNA distribution with dropout
We first examined differences among mRNA distributions generated by the telegraph model, ZINB model, and zero-inflated telegraph model. It has been established that the telegraph model can generate only three distribution types: a bimodal distribution with both zero and non-zero peaks, and unimodal distributions with either a sole zero peak or a sole non-zero peak [4]. We generated synthetic datasets containing 4,000 cells using the stochastic simulation algorithm (SSA) under the zero-inflated telegraph model with a technical zero rate . We then fitted these synthetic datasets to the telegraph and ZINB models separately using MLE to test their fitting performance; see Fig 1A. When the underlying core telegraph process yields either a bimodal or a non-zero-peak unimodal distribution, the zero-inflated telegraph model can output bimodal data. While the ZINB model can visually produce bimodal curves when fitted to such data, the valley dividing the zero peak and the secondary peak is fixed rigidly at mRNA count m = 1. This artifact causes misalignment in both the inter-peak valley shape and the location of the secondary peak. By contrast, the telegraph model may misclassify these bimodal datasets as unimodal. For parameters where the core telegraph process generates zero-peak unimodal distributions, the zero-inflated telegraph model preserves this unimodal signature. Yet both the telegraph and ZINB models fit such data poorly; notably, the ZINB model may yield bimodal profiles instead of the true unimodal shape.
(A) Fitting performance of the telegraph and ZINB models against mRNA distribution data simulated from the zero-inflated telegraph model. All synthetic datasets are generated with 4000 cell samples and fixed kb = 30 and . Three distinct parameter regimes are displayed: left panel , middle panel and right panel . (B) Synthetic distribution data (with dropout events) were generated using the telegraph model across sample sizes (N = 500, 1000, 2000, 4000) and with/without 20% extrinsic noise. All parameters satisfied and . The zero-inflated telegraph model maintains correct selection rates of nearly 70% or higher, outperforming the other two models.
We further examined whether the distributions generated by the three models are distinguishable via model fitting and selection. To this end, we randomly sampled 104 parameter sets for the telegraph model, with . We then employed the SSA to generate synthetic distribution data for each parameter set across varying cell sample sizes (N = 500, 1000, 2000, 4000). Next, we introduced the dropout model [23,35] to add dropout events to the synthetic distribution data (S1 Appendix). Each of the three models then fits these dropout-containing synthetic datasets via MLE, with the corrected Akaike information criterion (AICc) subsequently applied to identify the optimal model for describing the data (S1 Appendix). AICc imposes a penalty based on parameter count: 4 for the zero-inflated telegraph model and 3 for the other two. As illustrated in Fig 1B, over 70% of the synthetic datasets were correctly identified as best described by the zero-inflated telegraph model. This indicates that the other two models’ distributions fail to accurately capture the synthetic profiles, and their pronounced deviations cannot be masked by the AICc penalty from the zero-inflated telegraph model’s additional parameter.
To assess the robustness of this finding, we further introduced extrinsic noise into the telegraph model’s parameters : each parameter was redefined as a log-normally distributed random variable, retaining its original mean while setting the standard deviation to 20% of the mean (i.e., a 20% noise level) [34,36]. Synthetic distribution data with extrinsic noise were first generated using SSA; dropout events were then added to these data following the same procedure as above. Notably, while the selection rate of the zero-inflated telegraph model decreased slightly compared to the noise-free scenario, it remained above 70% across most sample sizes. Even for the smallest sample size N = 500, the rate was only marginally lower at 68.5%; see Fig 1B. Collectively, these results highlight the zero-inflated telegraph model’s robust reliability in describing complex mRNA distributions with technical artifacts (e.g., dropout) and biological noise. In contrast, the other two models lack sufficient accuracy to match the target distributions.
Zero-inflated telegraph model reliably estimates effective transcriptional burst frequency and size
We first evaluated how well the telegraph model, ZINB model, and zero-inflated telegraph model estimate burst frequency and burst size , using the 104 groups of dropout-containing synthetic distributions described in Fig 1B. For each model, we plotted 104 points in Fig 2A, where each point compares the estimated (or ) against its true value. Ideally, if all points cluster tightly along the 1:1 reference line, it indicates consistent estimates of effective burst parameters (with a Pearson correlation coefficient, PCC = 1). In contrast, points straying far from this line signal unreliable estimation, typically accompanied by an extremely low PCC.
Synthetic distribution data (with dropout events) were generated via the telegraph model across sample sizes (N = 500, 1000, 2000, 4000) and with/without 20% extrinsic noise. All parameters were constrained to and . (A) Estimation of and (no extrinsic noise, N = 4000): Zero-inflated telegraph model estimates cluster along the 1:1 line (PCC = 0.96); telegraph model underestimates and overestimates , while ZINB model shows opposite bias. (B) Estimation with 20% extrinsic noise: All models have reduced precision, but zero-inflated telegraph model remains robust (PCC > 0.75) with consistent bias patterns of the three models. (C) Estimation across different sample sizes (with/without noise): Zero-inflated telegraph model outperforms other models (PCC > 0.7, except 0.67 at N = 500 with 20% extrinsic noise); telegraph model and ZINB model struggles with inferring and , respectively. (D) distribution of two zero-inflated models: is similar and stable across sample sizes and extrinsic noise, confirming robust dropout rate inference.
At the large sample size N = 4000 (Fig 2A), the and estimates from the telegraph model departed markedly from the 1:1 line (PCC < 0.5): points consistently lay below the line (indicating underestimation), while points hovered above it (indicating overestimation). The ZINB model showed the opposite bias: overestimating and underestimating . In striking contrast, the zero-inflated telegraph model’s estimates for both and clustered tightly along the 1:1 line, with high PCC of 0.96. This underscores the model’s ability to reliably estimate burst parameters. These trends persisted when fitting dropout-containing synthetic data with 20% extrinsic noise added to each system parameter (Fig 2B). Notably, estimation accuracy degraded across all models relative to the noise-free scenario, yet their patterns of overestimation or underestimation for and remained unchanged. Crucially, the zero-inflated telegraph model maintained its advantage on inferring and : despite the noise-induced decline in precision, it retained robust performance with a PCC above 0.75, reinforcing its superiority in reliably estimating burst parameters even under noisy conditions.
We further compared the estimation of and across the three models under varying sample sizes N = 500, 1000, 2000, 4000; see Fig 2C and Fig A in S1 Appendix. Several consistent principles emerge. First, the zero-inflated telegraph model consistently outperforms the others, reliably estimating both and with PCC > 0.7 across almost all sample sizes, whether with or without extrinsic noise. The only exception is a modest PCC of 0.67 under the smallest sample size N = 500 with extrinsic noise. Second, the telegraph model struggles more with (PCC < 0.25) than with , while the ZINB model shows the reverse: it fares better with but poorly with (PCC < 0.25). This means neither model can simultaneously estimate both and reliably. Taken together, these findings highlight the zero-inflated telegraph model’s unique robustness: it maintains reliable parameter estimation across diverse sample sizes and noise conditions, whereas the other two models are inherently limited by their inability to balance accurate estimation of both burst frequency and size.
We finally examined the estimation of the technical zero rate for the ZINB model and zero-inflated telegraph model. Notably, there is no “true” value, as dropout-containing synthetic distributions were generated by randomly introducing dropout events for each gene and cell (S1 Appendix). Fig 2D shows estimates for both zero-inflated models, with and without extrinsic noise. The two models yielded consistently similar values, with the ZINB model estimating slightly higher values. This indicates that the zero-inflated telegraph model’s superior performance in describing synthetic data does not stem from better capturing dropout strength, but rather from the telegraph model’s greater power to capture inherent gene regulatory patterns. Equally striking, estimates remained stable across varying sample sizes and regardless of extrinsic noise. This underscores the robustness of zero-inflated models in estimating the dropout rate even when experimental conditions change. We also explore the correlation between the estimated and the estimation reliability of and . Using the 104 synthetic datasets with dropout (sample sizes N = 500, 1000, 2000, 4000, see Fig 1B), we calculated the absolute relative errors of and . The results reveal a weak positive correlation between these errors and inferred , with PCC ranging between 0.1 and 0.4 (Fig B in S1 Appendix).
Zero-inflated telegraph model exhibits no over-fitting
We investigated how the zero-inflated telegraph model estimates burst parameters when fitting synthetic data from simpler models, focusing on sample sizes N = 500, 1000, 2000, 4000. Using SSA, we generated 104 synthetic distributions from the telegraph model with parameters and , alongside 104 dropout-containing distributions from the ZINB model. Notably, both sets share same burst frequency and burst size .
For telegraph model-generated data, we fitted three models and plotted and estimates against their true values Fig 3A and Fig C in S1 Appendix: the telegraph model itself yielded accurate estimates with a PCC near 1. In contrast, the ZINB model’s estimates strayed from the 1:1 line, with overestimated and underestimated , consistent with Fig 2A. Critically, the zero-inflated telegraph model’s estimates clustered tightly along the 1:1 line (PCC approaching 1). Its estimated technical zero rate was tiny across all sample sizes () and smaller than that from the ZINB model (Fig 3B), highlighting its superior ability to distinguish dropout events. For dropout-containing data from the ZINB model, the telegraph model underestimated and overestimated , while both zero-inflated models produced similarly reliable burst parameter estimates Fig 3C and Fig C in S1 Appendix. Their estimates also aligned (Fig 3D), showing that the zero-inflated telegraph model captures dropout rates well even when fitting simpler data. Together, these findings confirm the zero-inflated telegraph model’s robust ability to balance flexibility to capture dropouts and parsimony to avoid overparameterization.
(A) For telegraph model-generated data (N = 4000), the zero-inflated telegraph model’s and estimates cluster tightly along the 1:1 reference line (PCC approaching 1), while the ZINB model shows obvious and estimation bias. (B) For telegraph model-generated data, the zero-inflated telegraph model has tiny () across sample sizes, smaller than those of the ZINB model. (C) For ZINB model-generated dropout-containing data (N = 4000), the telegraph model underestimates and overestimates , while both zero-inflated models give reliable and estimates. (D) For ZINB model-generated dropout-containing data, the two zero-inflated models have aligned estimates, confirming the zero-inflated telegraph model’s ability to capture dropout rates. For both the telegraph model and ZINB model, all parameters were constrained to and , with and .
Generalization evaluation across different gene regulatory models
In previous analyses, synthetic datasets were generated based on the basic telegraph model and negative binomial distribution. Inspired by general approach of generating synthetic data based on complex gene regulatory networks for algorithm validation [37], we further tested the generalization capacity of our zero-inflated telegraph model using the three-state model, as shown in the following diagram [38]:
The three-state model is a biologically motivated extension of the classical telegraph model. To fit experimental observations, it expands the single inactive state into two sequential states controlled by separate activation rates kon,1 and kon,2. This modification resolves the telegraph model’s drawback of assuming exponentially distributed inactive durations, which fails to explain the non-exponential peak patterns widely detected across prokaryotic and eukaryotic genes [39,40]. Mechanistically, the three-state structure mimics multi-step promoter events such as chromatin remodeling, and also reproduces the refractory state of transcription [15].
All synthetic datasets were generated using SSA with 4000 cells. We created 104 datasets for the three-state model with parameters sampled within fixed intervals: and . We set to highlight the multi-step activation mechanism that distinguishes the three-state model from the telegraph model. Specifically, the three-state model degenerates into the telegraph model when . Meanwhile, prior evidence confirms that equal values of kon,1 and kon,2 yield the maximum deviation of transcriptional noise from the telegraph model [15]. For the three-state model, we prepared two types of synthetic data: raw datasets free of technical zero inflation, and datasets with dropout events to simulate sequencing-derived biases. Next, we benchmarked the performance of the classical telegraph model, ZINB model, and our newly proposed zero-inflated telegraph model.
Three major conclusions were drawn from the comparative analyses. First, we adopted the Hellinger distance (HD) to quantify distribution fitting performance (Fig 4A). On datasets free of technical zeros, the zero-inflated telegraph model and the classical telegraph model achieved the lowest HD values and optimal fitting results, whereas the ZINB model performed the worst. When technical zero inflation was introduced, the zero-inflated telegraph model still maintained the minimal HD, followed by the ZINB model. This performance difference demonstrates that integrating a dedicated zero-inflation correction is indispensable for analyzing scRNA-seq data contaminated by dropout artifacts.
Data were generated with parameters constrained to and . (A) Hellinger distance (HD) results show that the zero-inflated telegraph model outperforms the classic telegraph model and ZINB model in distribution fitting. (B) Estimates of the technical zero rate confirm that the zero-inflated telegraph model accurately distinguishes biological zeros from technical artifacts. (C) The zero-inflated telegraph model exhibits relatively higher accuracy for inferring burst frequency and burst size of the three-state model. (D) The zero-inflated telegraph model effectively tracks relative changes in overall transcriptional kinetics of the three-state regulatory system, whereas the telegraph model and ZINB model fail to do so. (E) Comparative fitting performance on synthetic genetic toggle switch datasets. In the absence of technical dropout noise, the zero-inflated telegraph model yields lower HD and smaller estimated technical zero rate than the ZINB model; with dropout noise introduced, it retains minimal HD while generating broadly consistent estimates relative to ZINB.
Second, we evaluated the estimation of the technical zero rate (Fig 4B). For synthetic data without technical zeros, our model yielded extremely small values, which were lower than those estimated by the ZINB model. In contrast, the two zero-inflated models produced comparable values when fitting datasets with dropout events. This result confirms that our model can accurately distinguish biological zeros from sequencing-induced technical zeros, and reliably judge whether input data contains zero-inflation artifacts.
Third, we assessed the inference accuracy of burst frequency and burst size . On datasets without zero inflation, the zero-inflated telegraph model showed comparable accuracy to the telegraph model, and both performed considerably better than the ZINB model (Fig D in S1 Appendix). This indicates that the negative binomial distribution in ZINB is unable to capture the mechanisms of the three-state gene regulatory model. When zero inflation was introduced, the zero-inflated telegraph model still achieved satisfactory parameter estimation, with PCC values above 0.8 for both and . Its performance was far better than that of the telegraph model and ZINB (Fig 4C), demonstrating that modeling zero inflation and intrinsic regulatory mechanisms is essential for accurate parameter inference.
Though our model inevitably produces quantitative deviations when fitting data generated from more complex three-state model, we mainly examined whether its estimates reflect true parameter changes. We used the root sum of squares (RSS), defined as , to characterize overall parameter variations. The RSS estimated by our model was strongly positively correlated with the true values, while neither the telegraph model nor ZINB model could reproduce such trends (Fig 4D). These results verify that our model reliably captures relative changes in transcriptional kinetics for the three-state gene regulatory system.
To further examine the generalization capacity of the zero-inflated telegraph model across distinct gene regulatory architectures, we generated synthetic datasets using the canonical genetic toggle switch circuit [37], which describes a two-gene system where each gene product represses the transcription of its counterpart. We randomly sampled 5,000 parameter groups with repression rates confined to (0,3) and gene product synthesis rates within (10,30), while fixing the degradation rate at unity. Synthetic single-cell expression profiles containing 4,000 cells per condition were simulated via the Python package gillespy2, yielding a total of 104 synthetic distribution datasets of two genes. We then fitted all synthetic data using three competing frameworks: the telegraph model, ZINB model, and our zero-inflated telegraph model, with quantitative comparisons summarized in Fig 4E.
For synthetic datasets free of sequencing dropout artifacts, the zero-inflated telegraph model and telegraph model yielded comparable HD values, both lower than those obtained from the ZINB model. Meanwhile, the technical zero rate estimated by our model remained smaller than that of the ZINB model. When sequencing dropout noise was incorporated into simulated profiles, the zero-inflated telegraph model produced the lowest HD, with its inferred values comparable with those output by the ZINB model. Collectively, these observations indicate that the zero-inflated telegraph model maintains robust performance in distribution fitting and avoids overfitting risks, even when applied to expression data generated from separate regulatory architectures.
Zero-inflated telegraph model exhibits superior performance in DEGs analysis
We investigated the performance of three models-the telegraph model, ZINB model, and zero-inflated telegraph model-in identifying DEGs from single-cell mRNA distribution data. For this analysis, we utilized distribution datasets derived from a scRNA-seq experiment of peripheral blood mononuclear cells isolated from lupus patients treated with interferon-beta (IFN-) [31]. Specifically, we focused on CD14+ monocytes (a key myeloid subset responsive to IFN-) and analyzed 15,701 genes retained after quality control. From these 15,701 genes, we randomly selected 10,000 genes for subsequent model fitting and DEGs analysis.
We fitted the zero-inflated telegraph model for the 10,000 randomly selected genes to estimate three parameters . These parameter estimates served as the control group. To simulate DEGs, we first randomly selected 2,000 parameter sets from the 10,000 parameter estimates of the control group, and artificially perturbed one parameter by 2-fold in each selected set. These 2,000 perturbed parameter sets, together with the remaining 8,000 unchanged parameter sets from the control group, constituted the 10,000 parameter sets for the DEGs group. Using the telegraph model, we then generated dropout-containing synthetic distributions for both control and DEGs groups under four cell sample sizes: N = 100, 500, 1000, 2000. We applied the three models to the synthetic datasets to distinguish DEGs from non-DEGs, with implementation details following pipelines adapted from prior protocols [9,23]. Details of the DEGs detection pipeline are in the S1 Appendix. For performance quantification, we constructed Receiver Operating Characteristic (ROC) curves and computed the Area Under the ROC Curve (AUC) for each model at all four sample sizes: AUC = 1 indicates perfect DEGs detection, while AUC = 0.5 indicates random performance.
As shown in Fig 5, the zero-inflated telegraph model achieved the highest AUC across all sample sizes, followed by the ZINB model, with the telegraph model exhibiting the lowest performance. This result underscores two considerations for DEG analysis: (1) Accounting for zero inflation is essential to avoid misclassifying non-DEGs that arise from spurious zeros, and (2) incorporating the telegraph-like on-off switching dynamics improves detection of DEGs driven by changes in gene regulation kinetics. Notably, the AUC values of the two zero-inflated models converged at larger sample sizes. However, at the smallest sample size N = 100, the zero-inflated telegraph model outperformed the ZINB model by a noticeable AUC margin. This observation aligns with prior findings that gene expression in mammalian cells is inherently bursty [3,7], and further demonstrates that the zero-inflated telegraph model enables more accurate identification of DEGs, particularly in contexts with limited sample sizes.
Synthetic datasets simulating DEGs (from telegraph model-generated dropout-containing data, N = 100, 500, 1000, 2000) were used to test three models. The zero-inflated telegraph model achieved the highest AUC across all sample sizes (notably superior at small N), followed by the ZINB model, with the telegraph model performing lowest.
Zero-inflated telegraph model enhances DEGs detection in human scRNA-seq data
We validated the zero-inflated telegraph model using two distinct datasets, each composed of two comparative subgroups, followed by model fitting, parameter estimation, and DEGs analysis to verify the model’s performance. The first dataset focuses on CD14+ monocytes isolated from peripheral blood mononuclear cells (PBMCs) of systemic lupus erythematosus (SLE) patients [31]. It comprises two groups: 2,931 cells from the unstimulated control group and 2,765 cells from the IFN- treated group, respectively, with 15,701 genes profiled across all cells. These cells were derived from multiplexed droplet scRNA-seq experiments, which supports investigating cell-type-specific transcriptional responses to IFN- stimulation. The second dataset centers on tumor-infiltrating lymphocytes (TILs) isolated from patients with triple-negative breast cancer (TNBC) [32]. For comparative analysis, two CD8+ T cell clusters were selected: CD8+ tissue-resident memory T (CD8+Trm) cells and CD8+ non-tissue-resident memory (CD8+non-Trm) cells, comprising 606 and 1,097 cells, respectively, with 10,061 genes analyzed. These subsets were initially identified through scRNA-seq of TILs, and their transcriptional heterogeneity is crucial for elucidating the functional differences and prognostic significance of CD8+ T cell populations in TNBC.
We applied three models to fit the distribution data from 4 groups (derived from PBMCs and TILs) and computed HD to quantify the fitting performance. As shown in Fig 6A, the ZINB model yields significantly lower HD than the telegraph model, suggesting that capturing zero inflation may be more critical than accounting for the inherent gene state-switching mechanism when fitting genome-wide scRNA-seq transcriptomic data in human cells. Furthermore, the zero-inflated telegraph model outperforms the ZINB model with even lower HD values, indicating that incorporating detailed gene regulatory mechanisms also plays a pivotal role in accurate data fitting. To further dissect the balance between capturing zero inflation and gene state switching, we plotted the estimated technical zero rate against mean expression levels for both PBMCs and TILs datasets (Fig 6B). The results reveal that decreases with increasing mRNA mean, implying that dropout events are rare for observed highly expressed genes, and thus these genes may be better described by the telegraph model than the ZINB model. Collectively, our findings show that capturing zero inflation is more critical for fitting scRNA-seq distribution data, while accounting for inherent gene regulatory mechanisms further improves it. By integrating both, the zero-inflated telegraph model offers superior fitting for genome-wide scRNA-seq data.
(A) The zero-inflated telegraph model shows lower HD (better fitting) for mRNA distributions than the telegraph model and ZINB model. (B) The technical zero rate decreases with increasing mRNA mean for both datasets. (C) The zero-inflated telegraph model detects more DEGs than Seurat, the telegraph model, and the ZINB model, with notable advantages in small TIL samples. (D) Among the top 20 GO terms of the zero-inflated telegraph model, its DEGs share 14-16 GO terms with other models (marked in black), alongside unique terms (marked in red).
We compared DEGs analysis results from the three models and the conventional Seurat algorithm across PBMC and TIL datasets. As shown in Fig 6C, Seurat identified the fewest DEGs, while the telegraph model detected more than Seurat, and the two zero-inflated models yielded the most. Notably, in PBMC datasets (2,765 and 2,931 cells), the zero-inflated telegraph model identified slightly more DEGs than the ZINB model; in TIL datasets (606 and 1,097 cells), it detected significantly more. This aligns with Fig 5: the two zero-inflated models perform similarly with larger samples, but the zero-inflated telegraph model outperforms in small samples. Among unique DEGs of the zero-inflated telegraph model: In PBMCs, LDLR participates in low-density lipoprotein (LDL) metabolism, and its functional impairment causes plasma LDL accumulation, further inducing atherosclerosis-a common SLE comorbidity [41]; GZMB, a CD8+ T cell marker, associates with cytotoxicity and interferon signatures, key in lupus immune dysregulation [42]. In TILs, FAIM2, downregulated in tumors, links to immune infiltration and prognosis in breast cancer [43]; VDR regulates immune modulation, critical for breast cancer microenvironment [44]; CLSPN, a DNA damage response protein, promotes tumorigenesis via genomic instability, central to breast cancer progression [45].
Finally, we performed GO term enrichment analysis on DEGs identified by the three models in the PBMC and TIL datasets (Fig 6D and Fig E in S1 Appendix). Notably, 14 of the top 20 GO terms from the zero-inflated telegraph model overlapped with the top 20 terms of at least one other model in PBMCs, while this overlap reached 16 in TILs (Fig 6D), confirming that the model retains its advantage of identifying more DEGs without compromising functional enrichment consistency. These overlapping terms align with the research focus: in PBMCs, shared terms include multiple proteasome-related functions, involved in antigen processing/presentation, cellular stress responses, and signal regulation, closely associated with SLE pathogenesis [46,47]. In TILs, conserved terms include multiple ribosome-related functions, a process closely associated with breast cancer progression [48]. Beyond these shared functions, the zero-inflated telegraph model also enriched for unique GO terms not detected by the other models. For example, in PBMCs, “U2-type precatalytic spliceosome” is linked to SLE, as the spliceosome autoantigen heterogeneous nuclear ribonucleoprotein A2 serves as a major T cell autoantigen in SLE patients [49]. In TILs, “endoplasmic reticulum chaperone complex” is associated with breast cancer via its involvement in protein folding-related genes [50], while “proteasome core complex” enhances the efficacy of proteasome inhibitors to activate CD8+ T cell-mediated anti-tumor immunity in breast cancer [51].
Zero-inflated telegraph model quantifies technical zeros in scRNA-seq data from mouse cells.
Here, we investigated whether the zero-inflated telegraph model can improve the description of scRNA-seq data. We selected three bimodal mRNA distributions: the C57 and CAST alleles of Mbln2 gene in mouse fibroblasts [11], and the Plac/ara promoter in E. coli [2]. Both the telegraph model and zero-inflated telegraph model were employed to fit these distributions. As shown in Fig 7A, the telegraph model failed to accurately capture the bimodality of all three distribution data; it incorrectly fitted the C57 and CAST alleles of Mbln2 as unimodal distributions with non-zero peaks. A consistent limitation of the telegraph model fits was that the zero count in the fitted distribution was lower than the actual zero count in the data, suggesting the presence of notable technical zeros in the datasets. We therefore applied the zero-inflated telegraph model to these distributions, and it exhibited good agreement with the bimodal patterns, accurately matching both peaks; see Fig 7A. Incorporating zero inflation into the model thus significantly improves the fitting of bimodal mRNA distribution data.
(A) For single-cell transcription data (mouse fibroblasts Mbln2 C57/CAST alleles [11] and E. coli Plac/ara promoter [2]), the telegraph model misfits bimodal mRNA distributions as unimodal, while the zero-inflated telegraph model accurately captures bimodality. (B) For mouse ES cells genome-wide scRNA-seq data (Smart-seq2 and RamDA-seq) [11,17], the zero-inflated telegraph model shows lower HD (better data fitting) than the telegraph model, with smaller HD differences in Smart-seq2 data. (C) The technical zero rate estimated by the zero-inflated telegraph model is lower in Smart-seq2 data (median = 0.14-0.16) than in RamDA-seq data (median = 0.3-0.5). (D) Relative to the zero-inflated telegraph model, the telegraph model underestimates and overestimates , while the ZINB model shows the opposite bias. Compared with RamDA-seq data, the zero-inflated telegraph model has smaller differences in the estimated and with the telegraph model in Smart-seq2, but larger differences in these parameters with the ZINB model.
To further investigate the impact of zero inflation on scRNA-seq data fitting, we analyzed genome-wide CAST allele mRNA distribution data from mouse embryonic stem (ES) cells using two distinct single-cell sequencing strategies: a Smart-seq2 dataset covering 9,895 genes across 188 cells [11], and a RamDA-seq dataset comprising 9,182 genes quantified over 419 cells [17]. For rigorous comparison, we focused on the 1,988 genes shared between the two datasets. We fitted both datasets using the telegraph model and its extended zero-inflated variant. As anticipated, Fig 7B demonstrates that the zero-inflated telegraph model yielded lower HD values (better fitting performance) for both datasets. Notably, the difference in HD values between the two models was smaller for Smart-seq2 data than for RamDA-seq data, suggesting that Smart-seq2 data fitting is less affected by dropout events. This aligns with previous reports demonstrating that Smart-seq outperforms RamDA-seq for low-input RNA samples [52]. To validate this observation, we estimated the technical zero rate from the zero-inflated telegraph model. As shown in Fig 7C, Smart-seq2 data exhibited substantially lower values: its median was at 0.16, compared to 0.3 for RamDA-seq data. We further validated the above trend using a mouse ES cell dataset containing C57 and 129 alleles. This dataset includes 958 shared genes profiled by Smart-seq2 for the C57 allele (188 cells) [11] and by RamDA-seq for the 129 allele (419 cells) [17]. The median was 0.14 for Smart-seq2 and 0.51 for RamDA-seq; see Fig 7C. Overall, parameter estimates from these two independent datasets demonstrate that the technical dropout rate is consistently lower for Smart-seq2 than for RamDA-seq.
Finally, we compared burst frequency and burst size across gene sets derived from mouse ES cells: 1,988 shared CAST allele genes, and 958 shared genes covering both C57 and 129 alleles. Parameters were inferred from Smart-seq2 and RamDA-seq data using the telegraph model, ZINB model and zero-inflated telegraph model. For each model and dataset type, we retained only genes with estimated and values within [0.01,100], thereby excluding extreme outliers without biological significance [11,53]. We subsequently analyzed the distribution of these retained parameters; see Fig 7D. We note that the inferred corresponds to the ratio of the true to the degradation rate kd. We thus multiplied the inferred by the experimentally measured kd from mouse cell data [54] to obtain the actual value of . For the CAST allele, compared with the standard telegraph model, the zero-inflated telegraph model yielded distinct parameter estimates: the standard telegraph model produced lower (below the medians of 0.09 for Smart-seq2 and 0.2 for RamDA-seq data) and higher (above the medians of 1.54 for Smart-seq2 and 8.91 for RamDA-seq data). In contrast, the ZINB model gave larger (above the medians of 0.95 for Smart-seq2 and 0.12 for RamDA-seq data) and smaller (below the medians of 2.65 for Smart-seq2 and 2.15 for RamDA-seq data). These results are consistent with the trends shown in Fig 2A. Comparable quantitative discrepancies were also observed for the C57 and 129 alleles (Fig 7D). Notably, differences in (and ) between the telegraph model and zero-inflated telegraph model were less pronounced in Smart-seq2 data than in RamDA-seq data, indicating that the higher technical zero rate in RamDA-seq compromises the parameter accuracy of the telegraph model. In contrast, such differences were more striking between the ZINB model and zero-inflated telegraph model in Smart-seq2 data, indicating that the lower technical zero rate in Smart-seq2 highlights the need for more detailed gene regulatory mechanisms in model fitting.
Conclusions and discussions
In this study, we addressed a critical gap in single-cell transcriptional data analysis by bridging zero inflation correction and biophysical transcriptional modeling. Conventional zero-inflated models correct for dropout but rely on unimodal distributions [23], thus failing to capture transcriptional bursting-driven bimodality, a key feature of stochastic gene regulation [1,12]. The classical telegraph model describes stochastic gene state-switching dynamics yet conflates biological and technical zeros. Notably, its more complex variants suffer from low computational efficiency [55,56], while being sensitive to small sample sizes, extrinsic noise, and overfitting-all of which introduce biases into model selection and parameter inference [20,57,58].
To this end, we developed the zero-inflated telegraph model, a robust framework that integrates technical zero correction with the classical telegraph model’s inherent ability to capture stochastic transcriptional bursting dynamics. Through systematic validation across diverse datasets, including synthetic data, human scRNA-seq data from lupus patient PBMCs [31] and breast cancer TILs [32], as well as mouse scRNA-seq data from ES cells [11,17], we establish the model’s superiority over conventional approaches. Notably, zero-inflated telegraph model stands out for its enhanced biological interpretability, reliable inference of transcriptional regulation parameters, comprehensive detection of DEGs, and high computational efficiency, making it a practical tool for single-cell transcriptomic analysis.
The zero-inflated telegraph model demonstrates enhanced biological interpretability, reflected in its ability to capture mRNA distribution patterns and stabilize model selection. First, it enables stable model selection: even with a small sample size (N = 500) and strong extrinsic noise (20%), its correct selection rate remains nearly 70% or higher (Fig 1B). This ratio is comparable to that of the classical telegraph model [20,58], indicating that incorporating technical zeros does not compromise the performance of this more complex model. Second, compared with conventional models, it shows an essential quantitative difference in mRNA distribution characterization (Fig 1A): it not only accurately fits bimodal distributions in E. coli and mouse mRNA distribution data (where traditional models even misfit bimodal distributions as unimodal ones; Fig 7A) but also performs better when fitting genome-wide transcriptional data from human (Fig 6A) and mouse (Fig 7B) scRNA-seq data. Notably, the ZINB model outperforms the classical telegraph model, collectively highlighting that accounting for zero inflation is more crucial for fitting genome-wide transcriptomic data, while incorporating the inherent mechanism of gene state-switching further improves fitting efficacy.
The zero-inflated telegraph model enables reliable inference of transcriptional bursting parameters: burst frequency (standardized by degradation rate) and burst size . For one thing, it maintains robust and estimation when analyzing data integrating gene state-switching and dropout, even with small sample sizes (e.g., N = 500) and high extrinsic noise (e.g., 20%); see Fig 2A and 2B and 2C. For another, it degrades to the corresponding simple model when fitting synthetic data without dropout or inherent gene state-switching regulation, thus preventing overfitting (Fig 3). Notably, the ZINB model overestimates and underestimates , while the telegraph model shows the opposite bias (underestimating and overestimating ), a discrepancy that may explain the excessively large relative to inferred by the telegraph model in bacteria, yeast, and mammalian cells [7,8,53]. For example, in contrast to the telegraph model, which estimates smaller and larger for the Plac promoter under high-growth stimulation [8], RNA imaging directly measures Plac-based promoter as having and under the similar conditions [2] (S1 Appendix). At the genome-wide scale, in mouse ES cells scRNA-seq data, the zero-inflated telegraph model corrects parameter inferences from the the telegraph model (median : 0.09-0.27; median : 1.5-9) and the ZINB model (median : 0.07-0.95; median : 1.5-2.7) to more accurate values (Fig 7D).
The zero-inflated telegraph model also enables reliable inference of the technical zero rate : inferred values are nearly independent of model type, sample size variations, or extrinsic noise, and decline to zero in synthetic data without dropout events (Fig 2D and Fig 3B and 3D). Furthermore, we performed a sensitivity analysis of system parameters based on the first three moments. The parameter showed the highest sensitivity among all parameters (Fig F in S1 Appendix). This result demonstrates that accounting for zero inflation is essential for model construction. Applying inference to genome-wide human scRNA-seq transcriptomic data reveals higher for lowly expressed genes and lower for highly expressed genes (Fig 6B). Notably, scRNA-seq typically exhibits low expression levels [17]. Comparing across two scRNA-seq datasets of mouse cells reveals lower technical zero rates for the Smart-seq2 platform relative to RamDA-seq. Specifically, the median absolute differences in between the two platforms reached 0.14-0.37 in ES cells (Fig 7C). This leads to a key observation (Fig 7B and 7D): the telegraph model performs better in data fitting and estimation of and for Smart-seq2 data than for RamDA-seq data. This underscores the need for caution when using the telegraph model for scRNA-seq data; ideally, the zero-inflated telegraph model should be used to eliminate the impact of technical zeros.
We validated the generalization ability of the zero-inflated telegraph model using synthetic data generated from biologically more realistic three-state transcription model and genetic toggle switch circuit [37,39]. Benchmarks on datasets with and without sequencing dropout artifacts show that our model outperforms the telegraph model and ZINB model in both distribution fitting accuracy and the estimation reliability of the technical zero rate , allowing effective discrimination between biological and technical zeros (Fig 4A and 4B and 4E). Furthermore, our model also yields substantially higher accuracy in estimating burst frequency and burst size (Fig 4C and 4D). These results suggest that the zero-inflated telegraph model can better characterize both technical zeros and the intrinsic gene regulatory mechanisms.
As expected, the zero-inflated telegraph model cannot perfectly infer and for complex regulatory systems. Nevertheless, its relative kinetic changes faithfully suggest the transcriptional mechanistic variations of the three-state model (Fig 4D). This finding is also supported by previous studies. The telegraph model has reliably captured genome-wide single-cell transcriptional bursting in mouse ES cells, fibroblasts and human cardiomyocytes [11,17–19]. For genes under complex regulatory circuits, relative changes in its estimated effective and reflect genuine regulatory alterations in native systems. For instance, burst kinetics derived from the telegraph model successfully identifies synthetic negative feedback regulation in human kidney cells [20]. Recent evidence further demonstrates that shifts in effective or correspond to major changes in the true parameters of complex transcriptional mechanisms [19]. Collectively, although the zero-inflated telegraph model fails to output absolute accurate values of bursting parameters, it remains informative for dissecting intrinsic regulatory mechanisms. These merits render it a robust tool to detect differentially regulated genes across cellular contexts.
The zero-inflated telegraph model exhibits superior performance in detecting DEGs, outperforming the classical telegraph model, ZINB model, and conventional algorithms like Seurat. When applied to synthetic datasets simulating DEGs, the model achieved the highest AUC across all sample sizes, with particularly more notable advantages at the smaller sample size (Fig 5). In validation with human datasets, including CD14+ monocytes from lupus patient PBMCs and CD8+ T cell subsets from breast cancer TILs, the zero-inflated telegraph model detected more DEGs than its counterparts (Fig 6C). Notably, it captured disease-relevant unique DEGs that were missed by other methods, e.g., LDLR (links to SLE-atherosclerosis) and GZMB (associates with cytotoxicity and interferon signatures) in lupus PBMCs; FAIM2 (links to immune infiltration) and VDR (regulates tumor microenvironment) in breast cancer TILs. Functional enrichment analysis confirmed the model’s DEGs aligned with disease mechanisms (Fig 6D): 14–16 of its top 20 GO terms overlapped with those of other models, while it also uncovered unique terms, e.g., “U2-type precatalytic spliceosome” links to SLE autoantigens in PBMCs and “endoplasmic reticulum chaperone complex” is critical for breast cancer protein folding in TILs, which further validate its biological relevance.
KEGG enrichment analysis further validated the biological functions of the DEGs identified by the zero-inflated telegraph model. LDLR was significantly enriched in the hsa05417 lipid and atherosclerosis pathway (). As atherosclerosis is a prevalent cardiovascular complication in patients with SLE [41], this result solidifies the association between aberrant LDLR expression and SLE-related vascular pathological lesions. In comparison, although FAIM2 was not enriched in canonical KEGG pathways, our enrichment analysis of breast cancer datasets revealed its significant enrichment in the apoptotic module of the hsa04210 cell growth and death pathway (). Mechanistically, FAIM2 acts as a key Fas apoptosis inhibitory molecule that specifically regulates Fas-mediated apoptotic signaling [43]. Consistent with this function, FAIM2 drives breast cancer progression by inhibiting tumor cell apoptosis. Its interacting and co-expressed genes are enriched in diverse apoptosis and immune signaling pathways, including the p53 and TNF signaling pathways, as well as natural killer cell-mediated cytotoxicity, highlighting the core function of FAIM2 in modulating breast cancer apoptosis and sustaining tumor microenvironment homeostasis.
Notably, we simplify technical-dropout modeling with a unified global , which enables our model to maintain the high computational efficiency and numerical stability of the original telegraph model. Several alternative methods handle variable technical zeros without manually introducing zero-inflation parameters. One approach attributes missing reads to cell-specific Beta-distributed capture rates, where stochastic binomial downsampling intrinsically yields zero values without artificial zero inflation [59,60]. A second method separates true biological zeros from technical sequencing dropouts by measuring spliced (mature) and unspliced (nascent) RNAs [16], though standard scRNA-seq datasets do not provide such paired RNA quantification. As previous mechanistic studies suggested [61,62], combining nascent and mature RNA measurements can further improve parameter estimation of the model. Joint analysis of nascent and mature RNA distinguishes transcriptional elongation from post-transcriptional decay, reducing systematic biases in and caused by fitting solely on mature mRNA data. Drawing on these strategies, we will incorporate gene- and cell-specific alongside paired nascent and mature RNA data in future model development.
A key limitation of this study is that the classical telegraph model fails to fully capture the quantitative features of more complex regulatory mechanisms. Accordingly, its inferred burst frequency and burst size are effective parameters rather than true microscopic biochemical rates, which requires cautious biological interpretation. Future work should address this by incorporating additional regulatory mechanisms (e.g., multiple gene states, crosstalk signaling pathways [13,15]) into computationally scalable frameworks, potentially integrating machine learning with stochastic modeling to improve parameter inference efficiency [56,63]. Another limitation is the exclusive focus on the steady-state data, so subsequent research ought to move beyond this by coupling transcriptional burst dynamics with time-resolved single-cell measurements [2,64] to understand cell fate determination and pathological processes. Additionally, the model standardizes mRNA degradation rate for computational simplicity, overlooking its variations across genes or cell states [11,17,18]. Future work should integrate gene-specific degradation kinetics to cover more comprehensive gene expression regulation layers and correct such biases.
Methods
In the following, we prove that the definition (3) guarantees that for any given observed distribution data , there exists a unique parameter set that generates the data. We proceed by contradiction. Suppose there exist two distinct parameter sets and that yield the identical observed distribution for all . Denote the distributions generated by the telegraph model with parameters and as and , respectively. We first claim that . Otherwise, if , then (3) yields for all . Since the parameters are uniquely determined by the first three moments of [34], we immediately obtain . This contradicts the assumption that the two parameter sets are distinct. Hence, . Given , the definition (3) yields the following relation between the two distributions:
Moreover, both distributions satisfy the model recurrence relations [4]:
For all , substitute , and into the above two recurrence equations and subtract the resulting equations. The common scaling factor a cancels out, leading to
Rearranging terms gives the ratio of successive distribution terms:
(4)We now show that . Suppose, for contradiction, that . Taking the limit of (4) as , we obtain
(5)However, from the explicit integral representation (1) of the distribution, we have
(6)As , the right-hand side tends to 0, which contradicts the existence of a finite non-zero limit (5). Therefore, must hold.
Substituting into the ratio formula (4) yields
Combining this result with the integral bound (6), we arrive at
This result is impossible, as the left-hand side varies with m, while the right-hand side does not. Consequently, the assumption that two distinct parameter sets exist is invalid, which completes the proof of parameter uniqueness.
Code accessibility
For the analyses presented in this study, the R code implementing the zero-inflated telegraph model, including mRNA distribution calculation, data fitting, parameter estimation and DEGs detection, is publicly available at https://github.com/Always-Stude/ZITM.git This repository also contains scripts to generate synthetic datasets for both zero-inflated and standard model configurations.
Supporting information
S1 Appendix.
Fig A: Estimation of and under varying sample sizes and extrinsic noise. Synthetic data with dropouts were generated at N = 500, 1000, 2000 (complementing the N = 4000 results in Fig 2A and 2B). Estimation performance was compared among the telegraph model, ZINB model and zero-inflated telegraph model. Fig B: Correlation between the estimated technical zero rate and absolute relative errors of and . A weak positive correlation was found between parameter errors and , with PCC values between 0.1 and 0.4. Fig C: Model performance on synthetic data from the telegraph model and ZINB model. The zero-inflated telegraph model provides robust parameter estimates for both types of datasets, while the other two models present distinct estimation biases. Fig D: Burst parameter inference on synthetic dropout-free data from the three-state transcription model. The telegraph model and zero-inflated telegraph model show similar accuracy, both outperforming the ZINB model. Fig E: GO enrichment analysis of DEGs from lupus and breast cancer datasets. Shared and unique enriched GO terms are marked in black and red, respectively. Fig F: Sensitivity analysis of the zero-inflated telegraph model parameters using distribution moments. The technical zero rate is the most sensitive parameter, highlighting the importance of zero-inflation correction.
https://doi.org/10.1371/journal.pcbi.1014779.s001
(PDF)
Acknowledgments
We are grateful to Prof. Xiaofei Zhang and Prof. Da Zhou for their valuable suggestions.
References
- 1. Munsky B, Neuert G, van Oudenaarden A. Using gene expression noise to understand gene regulation. Science. 2012;336(6078):183–7. pmid:22499939
- 2. Golding I, Paulsson J, Zawilski SM, Cox EC. Real-time kinetics of gene activity in individual bacteria. Cell. 2005;123(6):1025–36. pmid:16360033
- 3. Raj A, Peskin CS, Tranchina D, Vargas DY, Tyagi S. Stochastic mRNA synthesis in mammalian cells. PLoS Biol. 2006;4(10):e309. pmid:17048983
- 4. Jiao F, Sun Q, Tang M, Yu J, Zheng B. Distribution Modes and Their Corresponding Parameter Regions in Stochastic Gene Transcription. SIAM J Appl Math. 2015;75(6):2396–420.
- 5. Razooky BS, Pai A, Aull K, Rouzine IM, Weinberger LS. A hardwired HIV latency program. Cell. 2015;160(5):990–1001. pmid:25723172
- 6. Zong C, So L, Sepúlveda LA, Skinner SO, Golding I. Lysogen stability is determined by the frequency of activity bursts from the fate-determining gene. Mol Syst Biol. 2010;6:440. pmid:21119634
- 7. Dar RD, Razooky BS, Singh A, Trimeloni TV, McCollum JM, Cox CD, et al. Transcriptional burst frequency and burst size are equally modulated across the human genome. Proc Natl Acad Sci U S A. 2012;109(43):17454–9. pmid:23064634
- 8. So L-H, Ghosh A, Zong C, Sepúlveda LA, Segev R, Golding I. General properties of transcriptional time series in Escherichia coli. Nat Genet. 2011;43(6):554–60. pmid:21532574
- 9. Vu TN, Wills QF, Kalari KR, Niu N, Wang L, Rantalainen M, et al. Beta-Poisson model for single-cell RNA-seq data analyses. Bioinformatics. 2016;32(14):2128–35. pmid:27153638
- 10. Delmans M, Hemberg M. Discrete distributional differential expression (D3E)--a tool for gene expression analysis of single-cell RNA-seq data. BMC Bioinformatics. 2016;17:110. pmid:26927822
- 11. Larsson AJM, Johnsson P, Hagemann-Jensen M, Hartmanis L, Faridani OR, Reinius B, et al. Genomic encoding of transcriptional burst kinetics. Nature. 2019;565(7738):251–4. pmid:30602787
- 12. Blake WJ, Balázsi G, Kohanski MA, Isaacs FJ, Murphy KF, Kuang Y, et al. Phenotypic consequences of promoter-mediated transcriptional noise. Mol Cell. 2006;24(6):853–65. pmid:17189188
- 13. Jiao F, Zhu C. Regulation of Gene Activation by Competitive Cross Talking Pathways. Biophys J. 2020;119(6):1204–14. pmid:32861266
- 14. Sheng Y, Lin G, Jiao F, Jia C. Geometric Theory of Distribution Shapes for Autoregulatory Gene Circuits. SIAM J Appl Math. 2025;85(2):636–61.
- 15. Zhou T, Zhang J. Analytical Results for a Multistate Gene Model. SIAM J Appl Math. 2012;72(3):789–818.
- 16. Carilli M, Gorin G, Choi Y, Chari T, Pachter L. Biophysical modeling with variational autoencoders for bimodal, single-cell RNA sequencing data. Nat Methods. 2024;21(8):1466–9. pmid:39054391
- 17. Ochiai H, Hayashi T, Umeda M, Yoshimura M, Harada A, Shimizu Y, et al. Genome-wide kinetic properties of transcriptional bursting in mouse embryonic stem cells. Sci Adv. 2020;6(25):eaaz6699. pmid:32596448
- 18. Trzaskoma P, Jung S, Pękowska A, Bohrer CH, Wang X, Naz F, et al. 3D chromatin architecture, BRD4, and Mediator have distinct roles in regulating genome-wide transcriptional bursting and gene network. Sci Adv. 2024;10(32):eadl4893. pmid:39121214
- 19. Chen L, Wu Y, Yang C, Fang S, Liao Y, Wu Y, et al. Using the simple telegraph model to decipher transcriptional burst regulation across genome-wide data. iScience. 2026;29(8):117127. pmid:42620746
- 20. Jiao F, et al. What can we learn when fitting a complex gene expression model to a simple telegraph? PLoS Comput Biol. 2024;20:e1012118.
- 21. Skinner SO, Xu H, Nagarkar-Jaiswal S, Freire PR, Zwaka TP, Golding I. Single-cell analysis of transcription kinetics across the cell cycle. Elife. 2016;5:e12175. pmid:26824388
- 22. Mahat DB, Tippens ND, Martin-Rufino JD, Waterton SK, Fu J, Blatt SE, et al. Single-cell nascent RNA sequencing unveils coordinated global transcription. Nature. 2024;631(8019):216–23. pmid:38839954
- 23. Miao Z, Deng K, Wang X, Zhang X. DEsingle for detecting three types of differential expression in single-cell RNA-seq data. Bioinformatics. 2018;34(18):3223–4. pmid:29688277
- 24. Silverman JD, Roche K, Mukherjee S, David LA. Naught all zeros in sequence count data are the same. Comput Struct Biotechnol J. 2020;18:2789–98. pmid:33101615
- 25. Jia C. Kinetic Foundation of the Zero-Inflated Negative Binomial Model for Single-Cell RNA Sequencing Data. SIAM J Appl Math. 2020;80(3):1336–55.
- 26. Chen L, Zhu C, Jiao F. A generalized moment-based method for estimating parameters of stochastic gene transcription. Math Biosci. 2022;345:108780. pmid:35085545
- 27. Li HS, Ou-Yang L, Zhu Y, Yan H, Zhang XF. scDEA: differential expression analysis in single-cell RNA-sequencing data via ensemble learning. Brief Bioinform. 2022;23:1–11.
- 28. Risso D, Perraudeau F, Gribkova S, Dudoit S, Vert J-P. A general and flexible method for signal extraction from single-cell RNA-seq data. Nat Commun. 2018;9(1):284. pmid:29348443
- 29. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. pmid:25516281
- 30. Lopez R, Regier J, Cole MB, Jordan MI, Yosef N. Deep generative modeling for single-cell transcriptomics. Nat Methods. 2018;15(12):1053–8. pmid:30504886
- 31. Kang HM, Subramaniam M, Targ S, Nguyen M, Maliskova L, McCarthy E, et al. Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nat Biotechnol. 2018;36(1):89–94. pmid:29227470
- 32. Savas P, Virassamy B, Ye C, Salim A, Mintoff CP, Caramia F, et al. Single-cell profiling of breast cancer T cells reveals a tissue-resident memory subset associated with improved prognosis. Nat Med. 2018;24(7):986–93. pmid:29942092
- 33.
Olver FWJ, Lozier DW, Boisvert RF, Clark CW. NIST Handbook of Mathematical Functions. New York: Cambridge University Press; 2010.
- 34. Grima R, Esmenjaud P-M. Quantifying and correcting bias in transcriptional parameter inference from single-cell data. Biophys J. 2024;123(1):4–30. pmid:37885177
- 35. Pierson E, Yau C. ZIFA: Dimensionality reduction for zero-inflated single-cell gene expression analysis. Genome Biol. 2015;16:241. pmid:26527291
- 36. Fu X, Patel HP, Coppola S, Xu L, Cao Z, Lenstra TL, et al. Quantifying how post-transcriptional noise and gene copy number variation bias transcriptional parameter inference from mRNA distributions. Elife. 2022;11:e82493. pmid:36250630
- 37. Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 2018;19(1):15. pmid:29409532
- 38. Tang M. The mean and noise of stochastic gene transcription. J Theor Biol. 2008;253(2):271–80. pmid:18472111
- 39. Suter DM, Molina N, Gatfield D, Schneider K, Schibler U, Naef F. Mammalian genes are transcribed with widely different bursting kinetics. Science. 2011;332(6028):472–4. pmid:21415320
- 40. Zimmer C, Häkkinen A, Ribeiro AS. Estimation of kinetic parameters of transcription from temporal single-RNA measurements. Math Biosci. 2016;271:146–53. pmid:26522167
- 41. Rivas L, Zabaleta M, Toro F, Bianco NE, De Sanctis JB. Decreased transcription, expression and function of low-density lipoprotein receptor in leukocytes from patients with systemic lupus erythematosus. Autoimmunity. 2009;42(4):266–8. pmid:19811272
- 42. Yin H, Li L, Feng X, Wang Z, Zheng M, Zhao J, et al. 2D4, a humanized monoclonal antibody targeting CD132, is a promising treatment for systemic lupus erythematosus. Signal Transduct Target Ther. 2024;9(1):323. pmid:39551768
- 43. Cai J, Ye Z, Hu Y, Wang Y, Ye L, Gao L, et al. FAIM2 is a potential pan-cancer biomarker for prognosis and immune infiltration. Front Oncol. 2022;12:998336. pmid:36185230
- 44. Chen P, Li M, Gu X, Liu Y, Li X, Li C, et al. Higher blood 25(OH)D level may reduce the breast cancer risk: evidence from a Chinese population based case-control study and meta-analysis of the observational studies. PLoS One. 2013;8(1):e49312. pmid:23382798
- 45. Azenha D, Hernandez-Perez S, Martin Y, Viegas MS, Martins A, Lopes MC, et al. Implications of CLSPN Variants in Cellular Function and Susceptibility to Cancer. Cancers (Basel). 2020;12(9):2396. pmid:32847043
- 46. Ichikawa HT, Conley T, Muchamuel T, Jiang J, Lee S, Owen T, et al. Beneficial effect of novel proteasome inhibitors in murine lupus via dual inhibition of type I interferon and autoantibody-secreting cells. Arthritis Rheum. 2012;64(2):493–503. pmid:21905015
- 47. Yang M, Wang P, Liu T, Zou X, Xia Y, Li C, et al. High throughput sequencing revealed enhanced cell cycle signaling in SLE patients. Sci Rep. 2023;13(1):159. pmid:36599883
- 48. Luan T, Song D, Liu J, Wei Y, Feng R, Wang X, et al. A Ribosome-Related Prognostic Signature of Breast Cancer Subtypes Based on Changes in Breast Cancer Patients’ Immunological Activity. Medicina (Kaunas). 2023;59(3):424. pmid:36984424
- 49. Fritsch-Stork R, Müllegger D, Skriner K, Jahn-Schmid B, Smolen JS, Steiner G. The spliceosomal autoantigen heterogeneous nuclear ribonucleoprotein A2 (hnRNP-A2) is a major T cell autoantigen in patients with systemic lupus erythematosus. Arthritis Res Ther. 2006;8(4):R118. pmid:16859514
- 50. Sisinni L, Pietrafesa M, Lepore S, Maddalena F, Condelli V, Esposito F, et al. Endoplasmic Reticulum Stress and Unfolded Protein Response in Breast Cancer: The Balance between Apoptosis and Autophagy and Its Role in Drug Resistance. Int J Mol Sci. 2019;20(4):857. pmid:30781465
- 51. Tang D, Lin S, Zhou J, Lei JH, Shao F, Sun H, et al. Augment proteasome inhibitor efficacy activates CD8+ T cell-mediated antitumor immunity in breast cancer. Cell Rep Med. 2025;6(7):102211. pmid:40609539
- 52. Ura H, Niida Y. Comparison of RNA-Sequencing Methods for Degraded RNA. Int J Mol Sci. 2024;25(11):6143. pmid:38892331
- 53. Cao Z, Grima R. Analytical distributions for detailed models of stochastic gene expression in eukaryotic cells. Proc Natl Acad Sci U S A. 2020;117(9):4682–92. pmid:32071224
- 54. Herzog VA, Reichholf B, Neumann T, Rescheneder P, Bhat P, Burkard TR, et al. Thiol-linked alkylation of RNA to assess expression dynamics. Nat Methods. 2017;14(12):1198–204. pmid:28945705
- 55. Jia C, Grima R. Holimap: an accurate and efficient method for solving stochastic gene network dynamics. Nat Commun. 2024;15(1):6557. pmid:39095346
- 56. Fang Z, Gupta A, Kumar S, Khammash M. Advanced methods for gene network identification and noise decomposition from single-cell data. Nat Commun. 2024;15(1):4911. pmid:38851792
- 57. Nicoll AG, Szavits-Nossan J, Evans MR, Grima R. Transient power-law behaviour following induction distinguishes between competing models of stochastic gene expression. Nat Commun. 2025;16(1):2833. pmid:40121209
- 58. Zhu C, et al. Correlation and distinction between stochastic gene transcription models with and without polymerase dynamics. Phys Rev Res. 2025;7:023050.
- 59. Wang Y, Shu Z, Cao Z, Grima R. From noise to models to numbers: Evaluating negative binomial models and parameter estimations in single-cell RNA-seq. PLoS Comput Biol. 2026;22(3):e1014014. pmid:41838800
- 60. Tang W, Jørgensen ACS, Marguerat S, Thomas P, Shahrezaei V. Modelling capture efficiency of single-cell RNA-sequencing data improves inference of transcriptome-wide burst kinetics. Bioinformatics. 2023;39(7):btad395. pmid:37354494
- 61. Wen K, Liao Y, Wang J, Choubey S, Jiao F. A moments-based approach for inferring mechanisms of transcriptional regulation using nascent RNA data. Biophys J. 2026;125(5):1257–75. pmid:41572627
- 62. Sullivan DK, Hjörleifsson KE, Swarna NP, Oakes C, Holley G, Melsted P, et al. Accurate quantification of nascent and mature RNAs from single-cell and single-nucleus RNA-seq. Nucleic Acids Res. 2025;53(1):gkae1137. pmid:39657125
- 63. Huang Z, et al. Deep learning linking mechanistic models to single-cell transcriptomics data reveals transcriptional bursting in response to DNA damage. eLife. 2025.
- 64. Semrau S, Goldmann JE, Soumillon M, Mikkelsen TS, Jaenisch R, van Oudenaarden A. Dynamics of lineage commitment revealed by single-cell transcriptomics of differentiating embryonic stem cells. Nat Commun. 2017;8(1):1096. pmid:29061959