Fitting Steady-State Phase Portraits

Fitting a Gene-Wise Steady-State Trend

A phase portrait places a gene’s unspliced abundance against its spliced abundance across cells. If the cells are near steady state and the gene’s kinetic rates are shared across them, their points should follow an approximately linear equilibrium relation. Fitting that relation is a key step in the steady-state velocity model: the fitted slope estimates a rate ratio, and a cell’s vertical deviation from the line can be converted into a spliced-RNA velocity if the splicing rate is available.

The fit is not simply a line-drawing exercise. Candidate cells must be plausibly near steady state, the chosen regression must match the equilibrium model, and the expression measurements must be on comparable scales. A line can be numerically well fitted yet biologically misleading if these conditions are not credible.

Definition: For one gene, let \(u_i\) and \(s_i\) denote the measured unspliced and spliced abundances, respectively, for cell \(i\). Under the two-compartment kinetic model, \(\frac{du}{dt}=\alpha-\beta u\) and \(\frac{ds}{dt}=\beta u-\gamma s\), where \(\alpha\) is the transcription rate, \(\beta\) is the splicing rate, and \(\gamma\) is the spliced-RNA degradation rate. At steady state, the expected phase-portrait relation is \(u_i=a s_i\), where \(a=\gamma/\beta\). A vertical residual from a fitted origin-constrained line is \(r_i=u_i-\hat a s_i\).

The equilibrium relation follows by setting both time derivatives to zero. The first equation gives \(\alpha=\beta u\), and the second gives \(\beta u=\gamma s\). Therefore \(u=(\gamma/\beta)s\). If transcription varies across cells while \(\beta\) and \(\gamma\) remain fixed, steady-state cells may occupy different locations along this same line. A single phase portrait generally identifies the ratio \(\gamma/\beta\), not the two rates separately.

Selecting Candidate Cells and Choosing a Fit

For each gene, define a candidate set \(C_g\) of cells whose measurements may reasonably represent the steady-state relation. Candidate selection is a scientific judgment supported by the data, not a guarantee that every selected cell is at equilibrium. Depending on the analysis, useful considerations include whether the gene has enough nonzero measurements to define a trend, whether the relevant cells span a usable range of abundance, and whether independent cell-state or time-course information makes equilibrium plausible.

Avoid selecting cells only because they lie close to the line you intend to fit. That practice makes the fit look better by construction and can remove precisely the cells whose deviations are important for estimating velocity. Equally, a cluster of near-zero measurements does not provide a reliable slope merely because it contains many cells. Record the selection rule and check whether the fitted relation is stable across reasonable changes to the candidate set.

The equilibrium equation has no intercept. This gives a mechanistic reason to fit through the origin when measured abundances are appropriately scaled and there is no relevant additive background offset. Ordinary least squares (OLS) then chooses the slope that minimizes the sum of squared vertical deviations. An unconstrained intercept can be useful as a diagnostic for a possible offset, but its fitted line is not the stated equilibrium equation. Robust regression can reduce the influence of extreme points, although it does not make an incorrectly selected candidate set or a non-steady-state model valid.

Derivation: For candidate cells \(i\in C_g\), the origin-constrained OLS slope minimizes \(Q(a)=\sum_{i\in C_g}(u_i-a s_i)^2\). Differentiating gives \(Q'(a)=-2\sum_{i\in C_g}s_i(u_i-a s_i)\). Setting this derivative to zero yields \(\hat a=\frac{\sum_{i\in C_g}s_i u_i}{\sum_{i\in C_g}s_i^2}\), provided \(\sum_{i\in C_g}s_i^2>0\). The second derivative is \(Q''(a)=2\sum_{i\in C_g}s_i^2>0\), so this solution is the unique minimum.

This calculation requires at least one nonzero spliced abundance, equivalently \(\sum_{i\in C_g}s_i^2>0\); a useful range of abundance may improve stability or informativeness but is not required for the algebraic solution. Interpreting OLS as estimating the conditional mean line also relies on a suitable linear relation and errors with mean zero conditional on \(s_i\); the usual classical uncertainty calculations additionally assume independent errors with constant variance. Those assumptions can be questionable for single-cell abundance measurements, whose variability often changes with expression. A descriptive least-squares slope can still be computed, but the formula alone does not establish that its uncertainty or kinetic interpretation is reliable.

Worked Example: Fitting an Origin-Constrained Trend

Suppose three hypothetical candidate cells have spliced abundances \(s=(1,2,3)\) and unspliced abundances \(u=(1.1,1.9,3.2)\), all in the same normalized abundance units. Assume for this calculation that the cells are plausible steady-state candidates and that an origin-constrained fit is appropriate.

The slope is \(\hat a=\frac{\sum s_i u_i}{\sum s_i^2}=\frac{1(1.1)+2(1.9)+3(3.2)}{1^2+2^2+3^2}=\frac{14.5}{14}=1.035714\). Thus the fitted relation is \(\hat u=1.035714s\). The slope is a ratio of unspliced to spliced abundance, so it has no units when both measurements use the same abundance scale.

Using the unrounded slope \(14.5/14\), the vertical residuals are \(1.1-(14.5/14)(1)=0.9/14=0.06429\), \(1.9-(14.5/14)(2)=-2.4/14=-0.1714\), and \(3.2-(14.5/14)(3)=1.3/14=0.09286\) abundance units. The residuals sum to zero after weighting by \(s_i\), as expected for an origin-constrained least-squares fit.

If the splicing rate were known to be \(\beta=0.40\) per hour, the corresponding model-based spliced velocities would be \(\hat v_{s,i}=\beta r_i\): \(0.40(0.9/14)\approx0.02571\), \(0.40(-2.4/14)\approx-0.06857\), and \(0.40(1.3/14)\approx0.03714\) abundance units per hour. These are conditional on the assumed rate and fitted equilibrium relation; the phase-portrait slope alone does not supply the value of \(\beta\).

From a Fitted Line to Velocity

For a cell with observed abundances \(u_i\) and \(s_i\), the model’s instantaneous spliced-RNA velocity is \(v_{s,i}=\frac{ds_i}{dt}=\beta u_i-\gamma s_i\). At equilibrium, the fitted slope estimates \(a=\gamma/\beta\), so \(\gamma\) is represented by \(\beta a\). Substitution gives the practical plug-in estimate \(\hat v_{s,i}=\beta(u_i-\hat a s_i)=\beta r_i\). The residual is in abundance units; multiplying by a splicing rate in inverse time gives velocity in abundance units per time.

This conversion has a crucial qualification: \(\hat a\) estimates a rate ratio, whereas \(\beta\) sets the time scale. If \(\beta\) is not known or estimated separately, \(r_i\) can be reported as a scaled velocity proxy, but not as an absolute velocity in abundance units per hour or another time unit. The conversion also presumes that the fitted relation represents the equilibrium line for the gene and that the measurements are comparable across cells.

A fitted intercept changes the interpretation. For a line \(\hat u=c+\hat a s\), the residual \(u_i-(c+\hat a s_i)\) is not, in general, the kinetic expression \(u_i-(\gamma/\beta)s_i\). An intercept may signal an additive measurement offset or other mismatch, but should not automatically be folded into the steady-state velocity equation. Investigate the reason for it before using residuals as velocities.

Worked Example: Comparing Origin and Intercept Fits

Consider four hypothetical candidate cells with \(s=(1,2,3,4)\) and \(u=(1.4,2.3,3.2,4.1)\) normalized abundance units. Assume the cells are plausible candidates and inspect both regression choices to see whether an additive offset matters.

For OLS with an intercept, \(\bar s=2.5\) and \(\bar u=2.75\). The centered cross-product sum is \(\sum(s_i-\bar s)(u_i-\bar u)=4.5\), and the centered squared-spliced sum is \(\sum(s_i-\bar s)^2=5\). Thus the slope is \(4.5/5=0.9\), and the intercept is \(2.75-0.9(2.5)=0.5\). The fitted line is \(\hat u=0.5+0.9s\).

For the fit constrained through the origin, the numerator is \(\sum s_i u_i=1.4+4.6+9.6+16.4=32\), and the denominator is \(\sum s_i^2=1+4+9+16=30\). The slope is \(32/30=1.0667\). The two fits differ appreciably because the data follow a line with a positive offset. Although that offset is informative diagnostically, the intercept fit is not the no-offset steady-state equation. It would be inappropriate to use its ordinary residuals as \(\beta(u_i-\hat a s_i)\) without resolving the mismatch.

Diagnostics, Alternatives, and Limitations

After fitting, inspect the points and residuals rather than relying only on the slope. A convincing fit should not be determined by one extreme cell, should cover a meaningful range of spliced abundance, and should not show a systematic curve or changing spread that contradicts the proposed linear trend. Compare the origin-constrained fit with an intercept fit as a diagnostic, and examine whether a small number of high-abundance cells dominate the sum of squares. If a narrow range or a few influential points control the estimate, the slope may be unstable even when a software routine returns a precise-looking number.

The alternatives answer different questions. Origin-constrained OLS respects the equilibrium model but can be biased by an uncorrected additive offset. OLS with an intercept can reveal such an offset, but its line does not itself define the model’s equilibrium velocity relation. Robust regression can limit the impact of isolated extreme measurements, but its result depends on the chosen robust loss and still requires appropriate candidate cells. None of these choices repairs poor normalization, inadequate expression range, or systematic departures from steady state.

A practical fitting sequence is:

Define candidates.
State why the selected cells are plausible steady-state observations for this gene; do not define them solely by proximity to the desired fitted line.

Check the measurements.
Confirm that unspliced and spliced abundances are on comparable scales, and that the candidate cells provide a useful range rather than only near-zero values.

Fit and inspect.
Estimate the origin-constrained slope when its assumptions are credible; compare alternatives and inspect points, residual patterns, and sensitivity to candidates.

Convert cautiously.
Multiply the vertical residual by a justified splicing rate for velocity in abundance per time, or label the residual as a scaled proxy if the rate is unavailable.

The steady-state fit is limited when candidate cells are not near equilibrium, when kinetic rates vary across cells, or when observed values are dominated by measurement noise or technical offsets. A visually straight cloud does not prove equilibrium, and a slope estimate does not identify absolute kinetic rates. These limitations matter because an apparent velocity can reflect a poor equilibrium reference rather than a reliable kinetic deviation.

Worked Example: Assessing an Influential Candidate

Suppose four hypothetical cells with \(s=(1,2,3,4)\) and \(u=(1.1,2.0,2.9,4.0)\) appear to follow a steady-state trend. A fifth cell has \(s=5\) and \(u=8\). Assume independent cell-state information raises a concern that the fifth cell may not be at steady state; do not classify it as an outlier solely because it is far from the fitted line.

Using only the first four cells, the origin-constrained slope is \(\frac{1.1+4+8.7+16}{1+4+9+16}=\frac{29.8}{30}=0.9933\). If the fifth cell is added, the slope becomes \(\frac{29.8+5(8)}{30+25}=\frac{69.8}{55}=1.269\). This large change shows that the estimate is sensitive to the candidate set. The fifth point’s residual under the five-cell fit is \(8-1.269(5)\approx1.655\) abundance units.

The calculation does not decide whether that cell should be excluded. The appropriate response is to investigate its status using information that is not simply the residual itself, and to report how the estimate changes under defensible candidate definitions. If there is no basis for distinguishing the fifth cell, the sensitivity is a limitation of the estimate, not permission to discard inconvenient data.

Common Mistakes

Assuming every observed cell is at steady state. The equilibrium line describes a model condition, not an automatic property of all cells in a dataset.
Using an intercept without changing the interpretation. A fitted offset can be diagnostically useful, but an intercept residual is not automatically the kinetic residual \(u-(\gamma/\beta)s\).
Calling a residual an absolute velocity. A vertical residual has abundance units. An absolute velocity requires multiplication by a splicing rate with known time units.
Removing points just because they worsen the fit. Candidate selection needs a defensible rationale; selecting by fit residual alone can conceal real departures and understate uncertainty.
Assuming a good-looking line proves the rates are identified. The steady-state slope estimates \(\gamma/\beta\), not \(\gamma\) and \(\beta\) separately.

Key takeaway: Fit the origin-constrained relation \(u=\hat a s\) only when the candidate cells and measurement scale make that equilibrium model plausible. Its vertical residual becomes a model-based spliced velocity only after multiplication by a justified splicing rate.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

Why does the steady-state phase-portrait relation have zero intercept under the two-compartment model?
A gene has a fitted slope \(\hat a=1.2\), and a cell has \(s=2.0\) and \(u=2.7\) abundance units. Calculate its vertical residual.
If the splicing rate in Question 2 is \(0.5\) per day, calculate the plug-in spliced velocity and give its units.
Why might an intercept fit be useful even though its residual is not automatically a steady-state kinetic residual?
What does a large change in slope when one candidate cell is added indicate, and what does it not establish?

Solutions. 1. Setting the derivatives to zero gives \(\beta u=\gamma s\), hence \(u=(\gamma/\beta)s\), which passes through the origin. 2. The residual is \(r=u-\hat a s=2.7-1.2(2.0)=0.30\) abundance units. 3. The velocity is \(\hat v_s=\beta r=0.5(0.30)=0.15\) abundance units per day. 4. An intercept fit can reveal an additive offset or a mismatch with the origin constraint, but its residual includes that fitted offset and need not equal \(u-(\gamma/\beta)s\). 5. It indicates sensitivity to candidate selection or influence; it does not, by itself, prove that the added cell is invalid or justify excluding it.