The chronODE framework for modelling multi-omic time series with ordinary differential equations and machine learning

Abstract

Many genome-wide studies capture isolated moments in cell differentiation or organismal development. Conversely, longitudinal studies provide a more direct way to study these kinetic processes. Here, we present an approach for modeling gene-expression and chromatin kinetics from such studies: chronODE, an interpretable framework based on ordinary differential equations. chronODE incorporates two parameters that capture biophysical constraints governing the initial cooperativity and later saturation in gene expression. These parameters group genes into three major kinetic patterns: accelerators, switchers, and decelerators. Applying chronODE to bulk and single-cell time-series data from mouse brain development reveals that most genes (~87%) follow simple logistic kinetics. Among them, genes with rapid acceleration and high saturation values are rare, highlighting biochemical limitations that prevent cells from attaining both simultaneously. Early- and late-emerging cell types display distinct kinetic patterns, with essential genes ramping up faster. Extending chronODE to chromatin, we find that genes regulated by both enhancer and silencer cis-regulatory elements are enriched in brain-specific functions. Finally, we develop a bidirectional recurrent neural network to predict changes in gene expression from corresponding chromatin changes, successfully capturing the cumulative effect of multiple regulatory elements. Overall, our framework allows investigation of the kinetics of gene regulation in diverse biological systems.

Fig. 1. Diagram illustrating how the logistic curve can be used to model time-series genomic signals.

Left panel: The curve represents a genomic signal (y axis) undergoing a logistic increase over time (x axis), modeled using the simplified form of the logistic ODE (Eq. 2) and with a positive growth rate constant (i.e., k* > 0). Note that the genomic signal y* represents a translated and normalized version of the original genomic signal z. The curve is constrained between a lower asymptote of 0 and an upper asymptote defined by the parameter b* (dashed blue line). t_start and t_end represent the initial and final time points, respectively, monitored by the experimental time course. The curve is composed of an acceleration phase followed by a deceleration phase, with the inflection point marking the transition between these phases at t_switch (solid red line). This curve can be used to model changes in genomic signals over time, such as gene expression (e.g., mRNA copy number measured via RNA-seq) or chromatin accessibility (e.g., number of nucleosome-free positions measured via ATAC-seq or DNase-seq). The arrows indicate how the number of mRNA copies and nucleosome-free positions increases during the acceleration and deceleration phases. Right panel: Analogous representation for a signal undergoing a logistic decrease (k* < 0). Elements of this figure were created in BioRender. Gerstein, M. (2025) https://BioRender.com/9ok1440.

Fig. 2. Kinetic characterization of the developing mouse brain transcriptome across three brain regions.

A Distribution of the kinetic parameters k (magnitude expressed in absolute value, y axis) and b (x axis) across all monotonic genes in the three brain regions. Genes are grouped into three quadrants (Q1–Q3) based on k-means clustering of their k & b combinations. No genes are found in Q4. B Lineplot showing, for each set of activated and repressed genes in Q1 through Q3 groups, the average expression (y axis) over time (x axis; PC post-conception). Note that the mean expression levels (expressed in ${\log }_2$(TPMs)) were rescaled to the range 0–100% to allow for comparison across Q1–Q3 groups. The error band corresponds to the standard deviation of the mean rescaled expression level (i.e., mean  ± SD). Q1 activated and repressed genes have an average expression profile that resembles the full logistic curve, whereas Q2 and Q3 genes align with the decelerator and accelerator parts of the curve, respectively. The arrow in Q3 genes indicates that these profiles remain far from reaching the saturation level b. C Examples of brain-related genes whose expression profile follows switcher, decelerator, and accelerator profiles, either increasing (upper panels) or decreasing (lower panels). Colored lines correspond to chronODE's fitted curve reconstructed using the kinetic parameters a, b, and k (see also Supplementary Fig. 2); gray dots correspond to raw data points. The red vertical line indicates the switching point (t_switch). D For each brain region (x axis), proportion (%) of genes (y axis) characterized by switcher (purple), decelerator (green), and accelerator (orange) profiles. E Number of genes (y axis) that show concordant (dashed black line) and discordant (continuous colored line) kinetic patterns across the three regions (x axis) (s switcher, d decelerator, a accelerator). In this case, we focused on the set of 4204 and 2489 genes that have consistently increasing (upper panel) or decreasing (lower panel) monotonic fits in all three regions, respectively.

Fig. 3. Cell type specificity of gene expression kinetics.

A Distribution of the kinetic parameters k (y axis) and b (x axis) for activated genes (k > 0) across 41 brain cell types. Cell types have been grouped based on their day of appearance (between 8 and 9, 9 and 10, 10 and 11, 12 and 13, and 14and 15 embryonic [E] days). B For each cell type (x axis) within these five groups, proportion of genes (%; y axis) that are characterized by switcher (purple), decelerator (green), and accelerator (orange) profiles. C Distribution of the kinetic parameter k (y axis) for essential (red) and non-essential (gray) genes (x axis). Box plots present the median as the center, 25% and 75% percentiles as box limits, and whiskers extending to the largest and smallest values within the 1.5 × interquartile range (IQR) of the box limits. We applied the Wilcoxon Rank-Sum test (two-sided, unpaired) to compute significant differences (***: p value  < 0.001 [p value = 4.2E-203]; n = 54,212; we note these are unreplicated single-cell data from Qiu et al.). D Proportion (%; y axis) of essential (Ess.) and non-essential (Non-Ess.) genes (x axis) belonging to the three kinetic classes (color-coded as in B).

Fig. 4. Genes linked to poly-pattern cCREs are more dynamic during brain development and are enriched in brain-specific functions.

A Schematic illustrating the distinction between mono-pattern and poly-pattern genes. Mono-pattern genes (left panel) can be linked to one or more cCREs, but if associated with multiple cCREs, all of these cCREs exhibit the same regulatory profile—either all activated (k > 0, red, upper panel) or all repressed (k < 0, blue, lower panel). Genes linked to a single cCRE represent the simplest case of mono-pattern regulation. In contrast, poly-pattern genes (right panel) are associated with multiple cCREs, which may be predominantly activated (upper panel) or repressed (lower panel). B Distributions of gene expression fold-change (${\log }_2$(TPM); y axis) for mono-pattern (mono-p.) and poly-pattern (poly-p.) genes (x axis) across the three brain regions. Box plots present the median as the center, 25% and 75% percentiles as box limits, and whiskers extending to the largest and smallest values within the 1.5  × IQR of the box limits. We applied the Wilcoxon Rank-Sum test (two-sided, unpaired) to compute significant differences (***: p value  < 0.001). Forebrain: p value = 7.6E-33 (n = 7992; 2 biological replicates); Midbrain: p value = 1.5E-5 (n = 5999; 2 biological replicates); Hindbrain: p value = 9.5E-16 (n = 6772; 2 biological replicates). C Proportion of distal cCREs (y axis) linked to mono- and poly-pattern genes (x axis). Distal cCREs are defined as those located  >2 Kb from the nearest annotated transcription start site of the gene. We report mean proportion and standard deviation across the three brain regions (n = 3). Overlayed black dots indicate the proportion of distal cCREs in each of the three regions. D Gene Ontology biological process terms (y axis) enriched among mono- and poly-pattern genes. Significance of GO terms' enrichment was determined with a hypergeometric test (two-sided). Only GO terms reporting a Benjamini-Hochberg adjusted p value < 0.01 were considered. Among these, we show in the panel the 15 most significant GO terms with the corresponding −log₁₀ p value (x axis).

Fig. 5. Predicting time-series gene expression from chromatin signals of associated cCREs.

A Architecture of the biRNN-based model. At each time point, the input for a specific multicCRE-gene pair i is a vector c of dimension 1 × m, where m represents the number of cCRE signals associated with the gene. The model generates an output vector, representing the gene expression levels g across all time points. B Density plot showing cross-gene expression correlation for each of the four regulatory mechanisms. Each scatterplot shows true (x axis) versus predicted (y axis) expression values (${\log }_2$-transformed TPMs) across all genes in the test set, alongside their degree of correlation as measured by Pearson’s correlation coefficient (r) and the corresponding Benjamini–Hochberg adjusted p value. C Representative examples of gene expression predictions over time demonstrating our model’s capability of capturing the three kinetic patterns (switchers, decelerators, and accelerators). True and predicted expression values are color-coded. Expression values correspond to ${\log }_2$-transformed TPMs. D Barplot showing, for each of the four regulatory mechanisms, the relative feature importance (computed using the SHAP method; y axis), for different cCREs associated with genes in the test set (x axis). The x axis represents the proximity rank of each cCRE to its associated gene, with rank 1 indicating the closest cCRE. SHAP importance values are shown for each of the eight time points (color-coded). The analysis highlights that cCREs closer to the gene have higher predictive importance at the initial and final time points, while more distant cCREs show greater importance at intermediate time points.