Time spent in the active state can be estimated from data in Table 1 of Ref. [29] since the mean active time is the inverse of ${\sigma }_b$.

This corresponds to a Markovian arrival process (MAP) in queueing theory [53]. Since not all MAPs are renewal processes [35], we focus on renewal MAPs where inter-arrival times follow a phase-type distribution [54].

In the 2-state model, this implies that burst size ($\rho /{\sigma }_{\mathrm{off}}$) scales with cell volume, consistent with experimental findings in eukaryotic cells [60, 61].

Our modeling approach does not explicitly capture fluctuations in mRNA counts due to the cell-cycle phase dependence of kinetic parameters, cell-cycle duration variability, doubling of gene copy number during replication, and cell division, which some other published models account for [29, 57]. A major reason for excluding such biological detail is that these models are time dependent, which makes the exact solution of the joint distribution of nuclear and cytoplasmic mRNA a difficult technical problem. Another reason is that these complex models have a considerably large number of parameters which makes it very challenging to reliably infer them using mRNA count data from a few thousand cells of one type (typical sample size).

Each cell follows the 2-state ENC model with transcription rate proportional to the normalized cell volume $\beta \sim \text{Gamma}(7.5,0.11)$. A gamma distribution or a similar unimodal distribution (such as a lognormal distribution) is a good approximation to the cell volume distribution when the cell-cycle duration is highly noisy; see Fig. 2A in Ref. [65]. Other papers have used such distributions of parameters to model extrinsic noise on the rates of the telegraph model [7, 66]. Delay SSA simulations [56] generated 15 batches of joint count data, each comprising 2000 cells with distinct kinetic parameters summarized in SM Table S2 [38].

The 2-state ENC model was fit using Eqs. (6)–(8) to infer the parameters $\{\overline{\rho },{\sigma }_{\mathrm{on}},{\sigma }_{\mathrm{off}},\tau \}$. For the telegraph model [12], transcription was assumed volume-dependent with four reactions: $G\stackrel{\overline{\rho }\beta }{\to }G+M$, $G\underset{{\sigma }_{\mathrm{on}}}{\overset{{\sigma }_{\mathrm{off}}}{\leftrightharpoons }}{G}^{*}$, and $M\stackrel{d=1}{\to }\varnothing$, where $M$ denotes total mRNA. The empirical PGF for total counts was defined as $\overline{G}(z)={\sum }_i{z}^{n_i+m_i}/n_c$, and the loss function adapted to $J(\varphi )={\int }_0^1|\mathcal{G}(z,\varphi )-\overline{G}(z){|}^2\mathrm{d}z$, with $\mathcal{G}(z,\varphi )$ computed by integrating the square of $G=M({\sigma }_{\mathrm{on}}/d,({\sigma }_{\mathrm{on}}+{\sigma }_{\mathrm{off}})/d,\overline{\rho }\beta (z-1)/d)$ over $\widehat{p}(\beta )$. The modified telegraph model requires inference of three parameters $\{\overline{\rho },{\sigma }_{\mathrm{on}},{\sigma }_{\mathrm{off}}\}$. The relative errors are averaged over four parameters for the 2-state ENC model and three for the telegraph model, with the mean and standard error computed across 100 replicate datasets.

The other fixed parameters are $\overline{\rho }=17$, $\tau =0.5$, $d=1$, and $n_c={10}^4$.

We excluded the 3-state model, as simulations (SM Fig. S4, Table S4 [38]) showed that, with sample sizes of a few thousand cells, it is difficult to distinguish it from the 2-state model.

Here $\mathrm{cor}(x,y)$ denotes the correlation coefficient between score vectors $x$ and $y$.