Worked Example: Comparing Kinetic Models

Why Compare Kinetic Models?

RNA velocity is a model-based estimate of how expression is changing, not a direct measurement of a cell's future state. A steady-state model and a dynamical model can assign different velocities to the same observed unspliced and spliced abundances. The difference matters especially around a transcriptional switch: a cell may be moving away from equilibrium, so a steady-state relationship can be a poor description even when it is easy to fit.

This tutorial compares the models on a small, hypothetical dataset. The calculations use scaled abundances and an abstract time unit. They illustrate how to derive predictions, inspect the phase portrait, and compare residuals. They do not establish that the chosen kinetic model is biologically correct.

Definition: Let $u$ and $s$ denote scaled unspliced and spliced abundances. In the standard two-compartment kinetic model, transcription produces unspliced RNA at rate $\alpha$, unspliced RNA is spliced at rate $\beta u$, and spliced RNA is removed at rate $\gamma s$: $$ \frac{du}{dt}=\alpha-\beta u,\qquad \frac{ds}{dt}=\beta u-\gamma s. $$ The unspliced and spliced velocities are $v_u=du/dt$ and $v_s=ds/dt$, respectively. They have units of abundance per time unit.

Assumptions and Model Predictions

The steady-state model assumes that the cells being compared are sufficiently close to equilibrium for both derivatives to be approximately zero. Setting $ds/dt=0$ gives $\beta u=\gamma s$, or $u=(\gamma/\beta)s$. At this equilibrium, $v_s=0$. This relation is not a general identity between $u$ and $s$: it applies only when the equilibrium assumption is reasonable and the rates are stable.

A dynamical model instead uses the differential equations over time. In the worked calculations below, the rates are fixed at $\beta=\gamma=1$. Transcription is assumed to switch from $\alpha=4$ to $\alpha=0$ at time $t=2$, and both abundances are assumed to follow the stated induction dynamics up to that switch, with distinct trajectories. These are simplifying assumptions for the hypothetical example, not estimates from a biological experiment. The model also treats scaled counts as noisy observations of an underlying smooth trajectory; it does not claim that measured counts are exact.

Before the switch, assume $u(0)=s(0)=0$ and constant transcription $\alpha=4$. Solving the first differential equation gives $u(t)=4(1-e^{-t})$. Substituting this into $ds/dt=u-s$ and solving with $s(0)=0$ gives the spliced trajectory. After the switch, transcription is zero, so unspliced abundance decays exponentially. The spliced trajectory after the switch follows from solving $ds/dt=u-s$ with the values at $t=2$.

Derivation: Before the switch, $u(t)=4(1-e^{-t})$. The spliced equation is $ds/dt+s=u(t)$. Multiplying by the integrating factor $e^t$ gives $$ \frac{d}{dt}\{e^t s(t)\}=4(e^t-1). $$ Integrating from $0$ to $t$ and using $s(0)=0$ yields $s(t)=4\{1-(1+t)e^{-t}\}$. At the switch, $u(2)=4(1-e^{-2})$ and $s(2)=4(1-3e^{-2})$. If $\tau=t-2\geq 0$, then $$ u(2+\tau)=u(2)e^{-\tau},\qquad s(2+\tau)=e^{-\tau}\{s(2)+u(2)\tau\}. $$ The second expression follows by solving $ds/d\tau+s=u(2)e^{-\tau}$ with initial value $s(2)$.

This derivation also clarifies the phase portrait. During induction, $u$ rises toward its equilibrium level while $s$ follows with a lag. After transcription switches off, $u$ falls immediately, but $s$ can continue to rise briefly if the incoming splicing flux $\beta u$ still exceeds removal $\gamma s$. Thus a trajectory in the $(u,s)$ plane need not lie on the steady-state line. A loop or bend in that plane can be evidence of changing dynamics, but it is not by itself proof of a particular switch or its timing.

Worked Example: Calculate the Trajectory and Velocity After the Switch

At $\tau=0.5$, or $t=2.5$, calculate the dynamical predictions and the spliced velocity. Use the model and switch values derived above.

At the switch, $u(2)=4(1-e^{-2})\approx3.4587$ and $s(2)=4(1-3e^{-2})\approx2.3760$ abundance units. Therefore,

u(2.5)=3.4587e^{-0.5}\approx2.0978,\qquad s(2.5)=e^{-0.5}(2.3760+3.4587(0.5))\approx2.4900.

The model-based spliced velocity is $v_s=u-s$, because $\beta=\gamma=1$. Substitution gives $v_s\approx2.0978-2.4900=-0.3922$ abundance units per time unit. The negative sign means that the modeled spliced abundance is decreasing at this point. It does not mean that spliced abundance is negative, nor does it establish that every cell with similar measured counts is decreasing.

Compare Predictions With Hypothetical Counts

The table gives six hypothetical observations. The first four are assigned to the induction segment and the last two to the post-switch segment. The count at $t=3$ is deliberately discordant with the assumed dynamical trajectory: it helps illustrate why a visible change in unspliced counts is not sufficient grounds for trusting a model.

Time	Segment	Observed $u$	Observed $s$
0.5	Induction	1.55	0.35
1.0	Induction	2.55	1.05
1.5	Induction	3.05	1.78
2.0	Switch time	3.50	2.38
2.5	After switch	2.10	2.50
3.0	After switch	1.25	3.1853

For a transparent descriptive comparison, define a residual as observed abundance minus model-predicted abundance. For $m$ observations, the root mean squared error is $\mathrm{RMSE}=\sqrt{\sum_{i=1}^{m}r_i^2/m}$. It summarizes prediction discrepancies in the same units as the abundance. Here it is a descriptive measure, not a formal test or a guarantee that the model assumptions hold.

Worked Example: Check Dynamical Predictions and Residuals

Evaluate the dynamical model at each listed time and calculate the RMSE separately for unspliced and spliced abundance. The predicted $u$ values, rounded to four significant figures, are $1.5739,\ 2.5285,\ 3.1075,\ 3.4587,\ 2.0978,\ 1.2724$. The first four follow $4(1-e^{-t})$; the last two follow $u(2)e^{-\tau}$.

The predicted $s$ values are $0.3608,\ 1.0570,\ 1.7687,\ 2.3760,\ 2.4900,\ 2.1464$, using the induction and post-switch formulas. For unspliced abundance, the residuals are approximately $-0.0239,\ 0.0215,\ -0.0575,\ 0.0413,\ 0.0022,\ -0.0224$. Their squared sum is approximately $0.0065514$, so

\mathrm{RMSE}_u=\sqrt{\frac{0.0065514}{6}}\approx0.03304

abundance units. For spliced abundance, the residuals are approximately $-0.0108,\ -0.0070,\ 0.0113,\ 0.0040,\ 0.0100,\ 1.0389$. The squared residuals sum to approximately $1.0796$, giving

\mathrm{RMSE}_s=\sqrt{\frac{1.0796}{6}}\approx0.4242

abundance units. Combining all twelve coordinate residuals gives $\sqrt{(0.0065514+1.0796)/12}\approx0.3009$ abundance units. The large final spliced residual dominates this combined summary. It signals that the chosen trajectory does not describe both coordinates well, even though its unspliced predictions are close to the six observations.

What the Steady-State Comparison Says

For these illustrative rates, the steady-state relation is $u=s$, and the steady-state spliced velocity is zero. Comparing observed $u$ with $s$ gives residuals $1.20,\ 1.50,\ 1.27,\ 1.12,\ -0.40,\ -1.9353$ abundance units. The RMSE for this fixed equal-rate relation is approximately $\sqrt{10.4627/6}=1.321$ abundance units. It does not capture the strong change in unspliced abundance across the switch.

This is a comparison with a specified steady-state relation, not a claim that every steady-state fitting procedure must have that RMSE. A fitted relationship may change the numerical residuals, but it cannot make the equilibrium assumption true when the observations represent a transition. Conversely, a dynamical model has more temporal structure, but its extra structure is useful only if its assumptions and predictions are compatible with the data.

Here the unspliced trajectory is consistent with a rise followed by a fall, while the last spliced observation continues rising far above the dynamical prediction. Thus the dynamical model is a more plausible description of the apparent unspliced switch than the fixed steady-state relation, but it is not a good joint explanation of both coordinates. A defensible report would label the $-0.3922$ estimate at $t=2.5$ as conditional on the stated dynamical model and avoid presenting it as an established biological direction.

Worked Example: Choose a Velocity Estimate With Appropriate Caution

Suppose the immediate question is whether spliced abundance is predicted to rise or fall at $t=2.5$. The steady-state model gives $v_s=0$ by its equilibrium assumption. The dynamical model gives $v_s=u-s\approx-0.3922$ abundance units per time unit.

The observations at $t=2.5$ show $u=2.10$ and $s=2.50$, consistent in sign with a negative spliced velocity under the chosen rates. The unspliced RMSE is also small, at $0.03304$ abundance units. These facts favor the dynamical estimate over a zero-velocity steady-state answer for describing this particular transient, provided the fixed-rate switch model is a reasonable approximation.

However, the spliced RMSE of $0.4242$ abundance units and the late residual of about $1.0389$ show substantial conflict with the assumed full trajectory. Therefore the appropriate conclusion is conditional: the dynamical model estimates a decline at $t=2.5$, but the small dataset does not justify treating that estimate as robust. If the question requires a well-supported velocity rather than a model-conditional calculation, neither estimate should be reported without qualification.

Diagnostics, Alternatives, and Limitations

A phase portrait plots observed or predicted $u$ against $s$, while a residual plot places prediction errors against time or fitted abundance. Both views can reveal information hidden by one summary statistic. In this example, the trajectory bends after the switch, and the large late spliced residual stands out. A low unspliced RMSE alone would conceal that problem.

Useful checks include whether residuals show a systematic run of positive or negative signs, whether one observation dominates the squared-error total, and whether the predicted trajectory has the expected direction in each segment. Also check that the time ordering and switch point are plausible; a sequence of cells is not automatically a measured time course. Sparse counts, scaling choices, unmodeled changes in degradation or splicing rates, and biological heterogeneity can all make a smooth, constant-rate trajectory inadequate.

The main alternative here is the steady-state relation: it is simpler and can be useful when the cells are near equilibrium, but it cannot represent a transient switch as a trajectory. The dynamical model explicitly represents time and can distinguish induction from repression, but depends on the assumed kinetic form, rates, and switch. More elaborate fitting does not automatically resolve poor fit or identify the true biological mechanism. Residuals and phase portraits diagnose incompatibility; they do not prove why it occurred.

Common pitfall: “The unspliced counts go down after the switch, so the dynamical velocity is reliable.” That conclusion ignores the spliced observations and the model's other assumptions. In this example, the unspliced fit is close but the spliced residual at $t=3$ is large. The correction is to inspect both coordinates, report which model and time point define the estimate, and qualify the result when the trajectory fails to explain the data jointly.

A Practical Comparison

State the candidate assumptions.
Use a steady-state relation only when near-equilibrium behavior is plausible. For a dynamical fit, specify the kinetic equations, rates, and any assumed switch.

Derive or evaluate the predictions.
Calculate $u$, $s$, and the corresponding velocities from the stated model before interpreting their signs.

Inspect residuals in both coordinates.
Use plots and summaries such as RMSE, and check whether a single point or a systematic pattern drives the mismatch.

Match the conclusion to the evidence.
Choose a model only to the extent that its assumptions and residuals support it. If neither is adequate, report the model-conditional estimate and its limitations rather than implying certainty.

Check Your Understanding

Answer each question using the definitions and hypothetical calculations above.

Why does the steady-state model give $v_s=0$, and what assumption is required?
At $\tau=0.5$, use the displayed predictions to calculate $v_s$. State its units and direction.
Which coordinate contributes most to the combined dynamical RMSE in the example, and why?
If $u$ falls after a switch but $s$ continues to rise briefly, must the dynamical model be inconsistent? Explain.
What conclusion is justified when unspliced residuals are small but a late spliced residual is large?

Solutions. (1) At equilibrium $ds/dt=0$, so $\beta u=\gamma s$, which makes $v_s=\beta u-\gamma s=0$; the assumption is that the state is sufficiently close to equilibrium. (2) $v_s\approx2.0978-2.4900=-0.3922$ abundance units per time unit, so the model predicts a decline. (3) Spliced abundance contributes most: its squared residual sum is about $1.0796$, far above the unspliced sum of $0.0065514$, largely because of the late observation. (4) No. Spliced abundance can still rise while $\beta u>\gamma s$, even as $u$ is falling; the equations, not the direction of $u$ alone, determine $ds/dt$. (5) The model describes the unspliced observations well but fails to explain both coordinates jointly. Any velocity estimate should be presented as conditional on the model, with the mismatch made explicit.