Likelihood-Based Kinetic Model Fitting

From Kinetic Trajectories to Count Likelihoods

A dynamical RNA velocity model describes how unspliced and spliced RNA change over time, while the data consist of counts measured in cells whose positions along the process are generally unknown. Likelihood-based fitting connects these two layers: it asks which kinetic rates and latent cell times make the observed counts most plausible under a specified count model.

For one gene, write \(u(t)\) and \(s(t)\) for the expected unspliced and spliced abundances at time \(t\). In the two-compartment model, transcription supplies unspliced RNA, splicing transfers it to the spliced compartment, and degradation removes spliced RNA:

Here \(\alpha(t)\) is the transcription rate, and \(\beta\) and \(\gamma\) are the splicing and degradation rates. The transcription rate may change at a switching time, giving different trajectory segments during induction and repression. Once the rates, switching time, and initial conditions are specified, these equations determine \(u(t)\) and \(s(t)\).

Let \(U_i\) and \(S_i\) be the observed unspliced and spliced counts for cell \(i\), and let \(\tau_i\) be that cell's latent time on the fitted trajectory. Counts need not be on the same scale as molecular abundances, so suppose known positive factors \(c_{Ui}\) and \(c_{Si}\) convert trajectory abundance to expected observed counts. These factors may represent fixed exposure or normalization offsets; estimating them jointly would add another source of ambiguity.

The Poisson model assumes counts are conditionally independent across cells and between the two compartments, and that each conditional variance equals its conditional mean. It also assumes the trajectory and offsets correctly specify the conditional means. These are modeling assumptions, not consequences of the kinetic equations. Real single-cell counts can be more variable or dependent than this model allows, so a high likelihood alone does not establish that the assumptions are adequate.

Constructing and Optimizing the Likelihood

Collect the kinetic quantities—such as \(\beta,\gamma\), transcription rates, switching time, and initial conditions—into \(\theta\). The latent times \(\tau_1,\ldots,\tau_n\) are also unknown. Under the stated observation model, the likelihood is a product of Poisson probabilities. Taking logarithms turns the product into a sum:

The maximum-likelihood fit chooses allowed parameter values and latent times that maximize \(\ell\). Equivalently, it minimizes the negative log-likelihood. The factorial terms depend on the observed counts but not on the parameters, so they can be omitted during optimization, though they belong in the full log-likelihood. Rates must be positive, means must be positive whenever a positive count is observed, and latent times must lie in the time domain of the trajectory. A practical optimizer may enforce positive rates by optimizing their logarithms and constrain times to a specified interval.

The latent times are not measured timestamps. They are fitted positions along the kinetic curve, estimated from the count pairs and shared kinetic parameters. Consequently, the optimization is joint: changing a rate changes the trajectory, and changing a cell's latent time changes the predicted means for that cell. If external ordering information is available, constraints can preserve that order; otherwise, an optimizer may assign cells to locations that improve the likelihood without respecting an intended biological ordering.

A Simple Rate Estimate and Its Derivation

The full model generally requires numerical optimization, but simple submodels show how likelihood optimization behaves. Suppose latent times \(t_i\) and a rate \(\beta\) are fixed, the initial unspliced abundance is zero, and transcription is constant at rate \(\alpha\). The unspliced solution is \(u(t_i)=\alpha f_i\), where \(f_i=(1-e^{-\beta t_i})/\beta\). With known count-scale factors \(c_i\), let the observed counts \(Y_i\) be conditionally independent Poisson with means \(\mu_i=c_i\alpha f_i\).

Ignoring terms that do not depend on \(\alpha\), the log-likelihood is \(\ell(\alpha)=\sum_i [Y_i\log(c_i\alpha f_i)-c_i\alpha f_i]\). Differentiating and setting the score equal to zero gives:

If \(\sum_iY_i>0\), require \(c_i f_i>0\) for every \(i\) with \(Y_i>0\), and \(\sum_i c_i f_i>0\); then the displayed estimate is the unique positive maximizer, since the second derivative is \(-\sum_iY_i/\alpha^2<0\) for \(\alpha>0\). If all counts are zero and \(\sum_i c_i f_i>0\), the likelihood has a supremum as \(\alpha\to0^+\) but no maximizer under the strict constraint \(\alpha>0\); if \(\sum_i c_i f_i=0\), it is constant in \(\alpha\). This derivation also illustrates why a fitted rate is tied to the assumed trajectory, times, and count scale.

Assessing Fit With Likelihood Contributions

A maximized likelihood is most useful when examined alongside the observations and fitted means. For a Poisson count \(y\) with fitted mean \(\widehat{\mu}\), its contribution to the Poisson deviance is \(d=2[y\log(y/\widehat{\mu})-(y-\widehat{\mu})]\) when \(y>0\), and \(d=2\widehat{\mu}\) when \(y=0\). The total deviance is the sum over counts. It compares the fitted model with a saturated model, which assigns each observed count its own fitted mean. A contribution is small when the fitted mean is close to the observed count and larger when they differ substantially.

The signed deviance residual is \(\operatorname{sign}(y-\widehat{\mu})\sqrt{d}\). Its sign records whether the observed count is above or below the fitted mean, while its magnitude summarizes the discrepancy on a deviance scale. These residuals are useful summaries, not a guarantee of a calibrated reference distribution: fitted latent times and kinetic parameters have been estimated from the same data.

Optimization Checks, Alternatives, and Limitations

Numerical fitting can fail even when the likelihood is correctly written. The objective may have multiple local optima, nearly flat directions, or optima near parameter boundaries. Latent time and kinetic rates can compensate for one another: for example, changing a rate may be partly offset by moving cells along the trajectory. The earlier material on identifiability explains why a good likelihood value does not necessarily imply precisely determined rates.

Useful checks include verifying that the optimizer converged from multiple initial values, confirming that rates and times satisfy their constraints, and plotting observed counts against fitted means or the fitted trajectory. Large, patterned residuals—such as systematic underprediction in one trajectory segment—suggest that the mean trajectory or observation model may be inadequate. Check also that the same cells and count definitions were used for the fitted likelihood and the displayed diagnostics.

The Poisson likelihood is a transparent baseline for count data. A Gaussian likelihood may be convenient for continuous measurements or transformed values, but it can assign probability to negative counts and needs a suitable variance model. A negative-binomial likelihood allows variance greater than the mean, often a relevant feature of count data, but requires a dispersion parameter and changes the optimization. These alternatives do not repair an incorrect kinetic trajectory automatically. Any likelihood-based approach remains conditional on its observation model, the chosen trajectory structure, and the treatment of cell-specific offsets and latent times.

A common mistake is to treat each cell's latent time as an observed value and then interpret the fitted kinetic rates as if timing were known exactly. The correction is to regard times and rates as jointly estimated quantities: uncertainty or weak identification in one can affect the other. Another mistake is to interpret the largest likelihood among a few trial times as the global optimum. A finite grid comparison only ranks the tried candidates; optimization or a sufficiently fine search over the allowed domain is needed to locate the maximum.