henhaobio
Sign in

bioRxivSystems BiologyPreprint

From time-course expression to gene regulation: direct linear ODE inference without finite-difference approximation

Xiaoqing HuangiD, Andersen Man Shun AngiD, Aatman Pushkarkumar VasoyaiD, Kun HuangiD, Teresa M. PrzytyckaiD, Yijie WangiD

PDFDOICode
In-depth full-text digest available · Jump to digest ↓

Abstract

Inferring gene regulation from time-course expression profiles is essential for understanding how cells transition between states during development, differentiation, and disease progression. Existing approaches often model expression dynamics with ordinary differential equations (ODEs). However, due to the computational complexity of directly solving these ODE models, most methods rely on finite-difference approximations of temporal derivatives, which can amplify measurement noise, introduce discretization bias, and lead to unstable or biased parameter estimates. To fill this gap, we develop the first computational method to directly learn a linear ODE model for gene regulation inference without relying on finite-difference approximations. We first formulate an optimization problem that directly exploits the closed-form solution of the linear ODE system. We then solve this problem via gradient descent, deriving analytical gradients with respect to the model parameters; these gradients involve matrix exponentials and integrals, which are challenging to directly compute. To make the computation efficient, we further use high-order Taylor approximations of the gradients whose truncation error is on the order of machine precision. In addition, we establish theoretical results demonstrating an inherent, non-vanishing gap between our exact solution and solutions derived from finite-difference approximations, which underscores the theoretical advantages of our approach. Finally, we demonstrate that our method consistently outperforms competing approaches on both simulated data and real-world scRNA-seq datasets in terms of AUROC. Our source codes can be accessed here: https://github.com/EJIUB/ExactLinearODE

In-depth summary

Claude Opus 5.5 · xhigh effort · 8 Nov 2026 · from the full text of version 1 (20 May 2026) · CC BY 4.0

Fitting a linear ODE with an equilibrium term via its matrix-exponential solution, not finite differences, raises TF-network AUROC over competing methods on simulated and scRNA-seq data.

ODE-based gene regulatory network inference from time-course or pseudotime expression usually substitutes finite differences for derivatives, which amplifies noise, adds discretization bias under large or irregular time gaps, and destabilizes parameter estimates. This work optimizes the closed-form solution directly, keeps a nonzero steady state that SCODE and GRISLI omit, and proves a non-vanishing finite-difference bias in the homogeneous case.

Key results

  1. ShownOn simulated linear-ODE networks with irregular sampling, the exact-solution method recovered nonzero interaction entries better than every dynamic baseline at all noise levels. It also improved when more time points were available.

    n = 30 TFs; SD1 (m = 30) and SD2 (m = 90), α from 0 to 0.10, 10 datasets per α (60 per setting). Beat SCODE, GRISILI, dynGENIE3, PROB and SWING in every setting; baselines hovered near AUROC 0.5 regardless of noise; SD2 consistently above SD1 at each α. Method AUROC means appear only in Fig. 2. · Data simulation procedure; Simulation data; Benchmarking results

  2. ShownOn scRNA-seq of mouse fibroblasts reprogrammed into induced neurons, ordered by pseudotime, the method ranked TF-TF edges against a TFRN gold standard with the top AUROC among both dynamic and steady-state methods.

    405 cells sampled at days 0, 2, 5 and 22, Monocle pseudotime, 100 TFs. Highest overall AUROC, by a clear margin per the authors (values shown only in Fig. 3(a)), against dynamic ODE/Granger and steady-state baselines. Per-TF AUROC in mouse: Maf 0.78, Atf3 0.76, Hmga2 0.70, Zfp238 0.73. · Dataset; Benchmarking results

  3. ShownThe same AUROC advantage held on two embryonic stem cell differentiation time courses, human ES to definitive endoderm and mouse ES to primitive endoderm. Each was scored as a 100-TF network against TFRN.

    Human ES to DE: 758 cells, 0 to 96 h; per-TF AUROC BNC2 0.72, NAGOG 0.74, NFE2L3 0.70, PITX2 0.69. Mouse ES to PrE: 456 cells, 0 to 72 h; BHLHE40 0.80, CREB3 0.85, RUNX1 0.76, SNAI1 0.68. Highest overall AUROC in both (Figs. 4(a), 5(a); values not tabulated); random level 0.5. · Dataset; Benchmarking results; Dataset; Benchmarking results

  4. ShownFor the homogeneous linear ODE fitted with squared loss, the optimal first-order finite-difference estimator of the interaction matrix is generally biased relative to the exact-solution estimator, and the gap does not vanish.

    Theorem 1 (proof in supplement): the bias of the finite-difference estimate grows linearly with the largest sampling step Δ and with the spectral condition number κ(XX^T). Wide time gaps or ill-conditioned data therefore enlarge the discrepancy from the true A. · 4 Finite–difference approximation vs. exact ODE…

  5. SuggestedInteraction matrices learned by the method regenerate expression trajectories that sit closer to the observations than those from SCODE, the only baseline that estimates numeric values of A.

    SD1 at α = 0.04: trajectories regenerated from the learned A matched observations more closely than SCODE's (Fig. 2). On the mouse MEF-to-iN data, fits followed transient (Maf), gradual (Atf3), declining (Hmga2) and sigmoidal (Zfp238) patterns, while SCODE gave an early sharp spike then near-constant levels; both ES datasets showed similarly better fits. Assessed visually only. · Benchmarking results; Benchmarking results; Benchmarking results; Benchmarking results

Methods, data and reproducibility

TF expression follows a linear ODE around a learnable equilibrium x* with interaction matrix A. A Huber loss on the matrix-exponential closed-form solution is minimized by gradient descent from A = 0, with gradients from Fréchet derivatives evaluated via high-order truncated Taylor series. Edges are ranked by |aij|.

  • Simulated irregularly sampled linear-ODE time courses, setting SD1 (m = 30 time points), noise α from 0 to 0.10 · n = 60 datasets (10 per α), 30 TFs each
  • Simulated irregularly sampled linear-ODE time courses, setting SD2 (m = 90 time points), noise α from 0 to 0.10 · n = 60 datasets (10 per α), 30 TFs each
  • scRNA-seq of mouse embryonic fibroblast reprogramming into induced neurons (days 0, 2, 5, 22) with Monocle pseudotime; 100-TF network · mouse · n = 405 cells
  • scRNA-seq of human ES cell differentiation to definitive endoderm (0, 12, 24, 36, 72, 96 h) with provided pseudotime; 100-TF network · human · n = 758 cells
  • scRNA-seq of mouse ES cell differentiation to primitive endoderm (0, 12, 24, 48, 72 h) with provided pseudotime; 100-TF network · mouse · n = 456 cells
  • Transcription Factor Regulatory Network (TFRN) database, source of gold-standard TF edges for AUROC
Methods & tools
ExactLinearODE (this work), SCODE, GRISLI (written GRISILI in benchmarks), PROB (Bayesian Lasso linear ODE), dynGENIE3, SWING (Granger causality), GENIE3, JUMP3, IDA (Bayesian network inference), MI, Cor, LM (lm) and Lasso (msgps) steady-state baselines, Monocle (pseudotime, from data providers), Al-Mohy and Higham Fréchet-derivative approach, with the series approximated via a cited Algorithm 6.4
Experimental conditions
About 1 minute on simulated data, around 6 minutes on the scRNA-seq datasets of Sections 5.2 and 5.4, around 10 minutes on that of Section 5.3; gradient cost O(n^3). No hardware stated.
Reproducibility gaps
The text gives no accessions or download links for the three scRNA-seq datasets or the TFRN gold standard. Proofs, Taylor error bounds and implementation details sit in supplementary materials that are not included here.
Limitations (5)
  • The Huber-loss objective is non-convex and the authors note that initialization is crucial; every experiment initializes A at the zero matrix. · stated by the authors
  • Theorem 1 covers only the homogeneous system (x* = 0) fitted with squared error. It corresponds to the full Huber model only when δ is very large and x* = 0. · stated by the authors
  • The Taylor approximation of the matrix exponential assumes a small ‖A‖2; larger norms need a scaling-and-squaring fallback, and the error estimates are left to the supplement. · stated by the authors
  • The simulated data come from the same linear ODE family the method fits (Hurwitz, diagonally dominant A), so the simulation favors it; baselines sitting near AUROC 0.5 may reflect model mismatch. · noted in this summary
  • Overall real-data AUROC values appear only in figures, with no variance or significance tests. Each dataset is scored as a single 100-TF network against TFRN, and the highlighted TFs were picked where the method scored high. · noted in this summary

Read this if you infer gene regulatory networks from time-course or pseudotime-ordered expression with ODE models and want an alternative to finite-difference derivative estimates.

Citations
0

Reproducibility ledger

Nobody has reported trying to reproduce this paper yet. Did you replicate an experiment or run the code? Your report saves the next lab a week.

Sign in to file a reproducibility report.

Related papers

Similarity Map
SelectedSame fieldCross-field
THISBeheler-Ama '25Zamanighomi '15Fehr '19Zhao '25Zhang '24Matsumoto '17
61% semantic matchSystems Biology
Maggie Beheler-Amass, Chrisopher A Jackson, David Gresham+1 more·bioRxiv·Sep 19, 2025· 3 citations