Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Recovering Transcriptional Switching Times

RNA Velocity · Tutorial 394 of 700

Recovering Transcriptional Switching Times

Learn how a dynamical RNA velocity model uses unspliced and spliced counts to estimate when transcription switches and where a cell lies along the resulting gene-specific trajectory.

Advanced 11 min read

What You'll Learn

  • Write the two-compartment kinetic equations and identify the units of each rate.
  • Derive the unspliced and spliced abundance curves before and after transcription switches off.
  • Recover elapsed time within an on or off phase from model predictions.
  • Use both RNA compartments to check phase assignments and identify model misfit.
  • Distinguish a gene-specific switching time from a cell’s latent position and from clock time.

From RNA Abundance to a Switching Time

Unspliced and spliced RNA provide a dynamical record of transcription. When transcription starts, unspliced RNA tends to rise first; spliced RNA responds as molecules are processed. When transcription stops, unspliced RNA is depleted by splicing and degradation, while spliced RNA declines more slowly if its degradation rate is lower. The dynamical model turns these different responses into a trajectory through the two-dimensional space of unspliced and spliced abundance.

The aim here is to understand how a gene’s transcriptional switch is represented on that trajectory and how a cell’s observed pair of abundances can be related to the switch. A switching time is a time coordinate in the model, not automatically a known time in hours since a biological event. In a snapshot experiment, cells are observed at once; their positions along the modeled trajectory are inferred from their RNA measurements.

Definition: For one gene, let \(u(t)\) and \(s(t)\) denote its unspliced and spliced RNA abundances at model time \(t\). The two-compartment kinetic model is $$ \frac{du}{dt}=\alpha(t)-\beta u(t), \qquad \frac{ds}{dt}=\beta u(t)-\gamma s(t). $$ Here \(\alpha(t)\) is the transcription rate, \(\beta\) is the splicing rate, and \(\gamma\) is the spliced-RNA degradation rate. If abundance is measured in abundance units and time in hours, \(\alpha\) has units of abundance per hour, while \(\beta\) and \(\gamma\) have units of inverse hours. The switch-off time \(t_s\) is the model time when transcription changes from an on-state rate \(\alpha>0\) to an off-state rate of zero.

The equations assume constant \(\alpha,\beta,\gamma\) within each phase, nonnegative rates, and a single abrupt switch between the phases being modeled. The displayed formulas below assume \(\beta>0\) and \(\gamma>0\); equal-rate cases are interpreted by limits. They also treat abundance as continuous and do not explicitly represent measurement noise, cell-specific rate changes, or multiple transcriptional pulses. Observed counts are noisy measurements of abundance, so fitted positions need not lie exactly on the ideal curve.

Deriving the Two Phases of the Trajectory

For clarity, set the transcriptional onset to \(t=0\) and suppose \(u(0)=s(0)=0\). During the on phase, \(0\le t\le t_s\), transcription occurs at the constant rate \(\alpha\). Solving the unspliced equation gives \(u(t)=\alpha(1-e^{-\beta t})/\beta\). The spliced trajectory follows by integrating the splicing input \(\beta u(r)\), discounted by spliced-RNA degradation from time \(r\) to \(t\):

$$ s(t)=\int_0^t \beta u(r)e^{-\gamma(t-r)}\,dr =\frac{\alpha}{\gamma}(1-e^{-\gamma t}) -\frac{\alpha(e^{-\beta t}-e^{-\gamma t})}{\gamma-\beta}, \qquad \beta\ne\gamma. $$

At the switch, write \(u_s=u(t_s)\) and \(s_s=s(t_s)\). These are the abundances at the point where transcription turns off. After the switch, let \(\Delta=t-t_s\ge 0\). With \(\alpha(t)=0\), unspliced RNA decays exponentially. Spliced RNA receives input from the remaining unspliced pool while that pool is being depleted:

$$ u(t)=u_s e^{-\beta\Delta}, \qquad s(t)=s_s e^{-\gamma\Delta} +\frac{\beta u_s(e^{-\beta\Delta}-e^{-\gamma\Delta})}{\gamma-\beta}, \qquad \beta\ne\gamma. $$

These expressions meet continuously at the switch: at \(\Delta=0\), they give \(u=u_s\) and \(s=s_s\). If \(\beta=\gamma\), the expressions with denominator \(\gamma-\beta\) are replaced by their limits. In particular, on the off phase, \(s(t)=e^{-\beta\Delta}(s_s+\beta u_s\Delta)\); on the on phase, \(s(t)=\alpha[1-(1+\beta t)e^{-\beta t}]/\beta\).

Justification: During the off phase, multiplying \(ds/dt+\gamma s=\beta u_s e^{-\beta\Delta}\) by \(e^{\gamma\Delta}\) gives $$ \frac{d}{d\Delta}\bigl(e^{\gamma\Delta}s\bigr) =\beta u_s e^{(\gamma-\beta)\Delta}. $$ Integrating from \(0\) to \(\Delta\) yields $$ e^{\gamma\Delta}s-s_s =\frac{\beta u_s(e^{(\gamma-\beta)\Delta}-1)}{\gamma-\beta}. $$ Multiplication by \(e^{-\gamma\Delta}\) gives the off-phase expression above. Thus the second term in \(s(t)\) represents RNA newly spliced after transcription has stopped; simply letting \(s\) decay from \(s_s\) would omit that contribution.

If onset occurs at a nonzero time \(t_{\mathrm{on}}\), replace \(t\) in the on-phase formulas by the elapsed time \(t-t_{\mathrm{on}}\). The switch time is then \(t_s=t_{\mathrm{on}}+\tau_s\), where \(\tau_s\) is the duration of the on phase. More complicated initial conditions change the trajectory; the zero-abundance onset is a modeling assumption, not a fact guaranteed by the data.

Worked Example: Calculating the Abundances at Switch-Off

Suppose a hypothetical gene has \(\alpha=12\) abundance units per hour, \(\beta=1\) per hour, and \(\gamma=0.5\) per hour. Assume transcription begins from zero abundance and switches off after \(t_s=2\) hours. Calculate the switch-point abundances.

The on-phase unspliced abundance is \(u_s=12(1-e^{-2})=12(1-0.1353)=10.38\) abundance units.

For spliced RNA, substitute \(t_s=2\): \(s_s=12(1-e^{-1})/0.5-12(e^{-2}-e^{-1})/(0.5-1)\). Evaluating the expression with unrounded exponential values gives \(s_s\approx9.590\) abundance units.

The switch point is approximately \((u_s,s_s)=(10.38,9.590)\). This is a point on the gene’s trajectory, not a claim that an experimentally observed cell must have exactly those abundances.

Recovering a Cell’s Position Relative to the Switch

If the kinetic rates and phase are known, abundance can be translated into elapsed model time. In the on phase, the unspliced equation can be inverted:

$$ t=-\frac{1}{\beta}\log\left(1-\frac{\beta u}{\alpha}\right), \qquad 0\le u<\frac{\alpha}{\beta}. $$

The range condition matters: the on-phase unspliced abundance approaches \(\alpha/\beta\) asymptotically, so a value at or above that limit cannot yield a finite on-phase time under this model. The spliced abundance provides a second check: the observed pair \((u,s)\) should be compatible with the same time on the on-phase curve.

On the off phase, if the switch-point abundance \(u_s\) is known, unspliced abundance gives elapsed time since switch-off: \(\Delta=-\log(u/u_s)/\beta\), for \(0<u\le u_s\). The predicted spliced abundance at that same \(\Delta\) can then be compared with the observed \(s\). In practice, both compartments help constrain the phase and position; a single abundance alone may be ambiguous when rates, the switch point, or phase are uncertain.

Worked Example: Inferring Time in the On Phase

Consider a hypothetical gene with \(\alpha=8\) abundance units per hour and \(\beta=2\) per hour. A cell has \(u=3.20\) abundance units. Assume its unspliced abundance follows the zero-start on-phase model and that its phase is known to be on.

The limiting on-phase abundance is \(\alpha/\beta=8/2=4.00\) units, so \(u=3.20\) is within the allowable range. The inferred elapsed time since onset is \(t=-\log(1-2(3.20)/8)/2=-\log(0.20)/2=0.8047\) hours.

As a two-compartment check, let \(\gamma=1\) per hour. The predicted spliced abundance at \(t=0.8047\) hours is \(s=8(1-e^{-0.8047})-8(e^{-1.6094}-e^{-0.8047})/(1-2)\). Using the rounded exponential values, the second contribution is \(-1.978\), so the displayed rounded terms give \(4.422-1.978=2.444\); evaluating the unrounded expression gives \(s\approx2.445\) abundance units. A measured spliced abundance near 2.445 would agree with this idealized on-phase position; a substantial discrepancy would prompt a phase or model-fit check.

Worked Example: Inferring Time Since Switch-Off

Use the hypothetical switch point from the first example: \(u_s=10.38\) and \(s_s=9.590\), with \(\beta=1\) per hour and \(\gamma=0.5\) per hour. A cell has \(u=3.817\) abundance units, and the model assigns it to the off phase. Estimate its time since switch-off and check its spliced abundance.

From unspliced RNA, \(\Delta=-\log(3.817/10.38)/1\approx1.000\) hour. At this elapsed time, the off-phase prediction is \(s=9.590e^{-0.5} + 10.38(e^{-1}-e^{-0.5})/(0.5-1)\). Using the unrounded exponential values gives terms of about 5.8166 and 4.9544, respectively, so \(s\approx10.771\), or 10.77 abundance units. A measured value near 10.77 supports the proposed off-phase location under these parameters. It does not establish that exactly one clock hour has passed in the biological sample; it is one hour on this model’s time scale.

What the Switching Time Means—and What It Does Not

For a gene, \(t_s\) is the model time at which its transcription rate changes. Its switch-point abundances are determined by the full history up to that time, not by setting either RNA compartment to zero. A cell’s inferred position may be before or after that switch. If onset is set to zero, \(t_s\) is the duration of the modeled on phase; a cell on the off phase has model time \(t_s+\Delta\).

The switching time is gene-specific, whereas a cell’s position is its location on that gene’s trajectory. Across genes, different switches can occur at different positions along a broader cell-state progression. The earlier material on inferred pseudotime concerns positions along a cell-state path; a kinetic time coordinate for one gene is not automatically the same as that path’s pseudotime. Nor does either coordinate, without external calibration, identify elapsed clock time or prove that one cell biologically became another.

The main contrast is with the steady-state velocity model. That model describes deviations from a steady-state relationship between compartments; the dynamical model instead permits a transient trajectory with distinct on and off phases. This makes switching times meaningful in the dynamical model, but it also makes conclusions dependent on the assumed kinetic rates, initial conditions, and switch structure. The likelihood-based fitting material addresses parameter estimation; here the emphasis is how a given fitted trajectory encodes switch and cell times.

Diagnostics, Failure Modes, and Limitations

A time estimate is meaningful only if its assumptions and phase assignment are plausible. Useful checks include whether a proposed on-phase \(u\) is below \(\alpha/\beta\), whether off-phase \(u\) is no greater than \(u_s\), and whether the predicted \(s\) at the inferred time is compatible with the measured \(s\). A point that cannot be placed sensibly on either branch may reflect noise, an incorrect phase assignment, inappropriate rates, or a trajectory that is not well represented by a single on-to-off event.

  • Forgetting residual splicing after switch-off. Setting \(s=s_s e^{-\gamma\Delta}\) assumes no additional input into the spliced compartment. The correct off-phase solution includes the contribution from unspliced molecules still being spliced.
  • Confusing time scales. A fitted coordinate of \(0.8\) hours is not necessarily 0.8 hours of measured biological time. It is in hours only to the extent that the rate scale and model calibration justify that interpretation.
  • Using one compartment as decisive evidence. Similar unspliced values can occur on different portions of a trajectory. Check the paired spliced abundance and the model’s predicted phase curve.
  • Ignoring rate or initial-condition mismatch. Variable rates, nonzero starting abundance, multiple transcriptional pulses, or technical distortion can make a simple switch time misleading.

The dynamical approach is useful when cells span transient activation and repression, while a steady-state approach may be more suitable when the data support a stable relationship and transient timing is not the aim. Neither approach repairs violations of its assumptions. A trajectory that fits poorly should not be forced into a precise switching-time interpretation; model adequacy and uncertainty require separate assessment.

Key takeaway: The switch-off time identifies the point where the modeled transcription rate changes, while the paired unspliced and spliced abundances locate a cell before or after that point. The inferred time is a kinetic coordinate conditional on the model, rates, and initial conditions—not automatically a biological clock time.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. If \(\alpha=10\) abundance units per hour and \(\beta=2\) per hour, what is the limiting on-phase unspliced abundance? Can \(u=5.2\) yield a finite on-phase time under the zero-start model?
  2. For \(\beta=1\) per hour, \(u_s=8\), and off-phase \(u=4\), find the elapsed time since switch-off.
  3. Why is \(s=s_s e^{-\gamma\Delta}\) generally incorrect after transcription turns off?
  4. Does an inferred gene-specific switching time establish the actual clock time of a cell transition? Explain.

Solutions. (1) The limit is \(\alpha/\beta=10/2=5\) abundance units. Since \(5.2>5\), the on-phase inversion has no finite solution. (2) \(\Delta=-\log(4/8)/1=\log 2=0.6931\) hours. (3) Unspliced RNA present after switch-off continues to be spliced, adding input to \(s\); the off-phase formula includes this contribution. (4) No. It is a model time conditional on the kinetic assumptions and calibration, and is not automatically an observed clock time or proof of a biological transition.

Clinical Biostats