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

chronODE: modeling time-resolved genomic signals with the logistic function

Given the cooperative and saturating nature of genomic signals (Fig. 1), we propose modeling transcriptomic and epigenomic kinetics over time using the logistic function, previously applied in fields like bacterial growth^31,32. In this context, the signal of a genomic locus (e.g., gene expression, chromatin accessibility, or histone modification) is defined as a time-dependent positive variable z, constrained within the interval [a, b]. The rate of change of z over time t can be studied using the following formula, which corresponds to the generalized form of the logistic ODE (see Supplementary Note 1, Propositions 1 and 5, and Supplementary Fig. 2):

where a and b correspond to the lower and upper asymptotes, respectively, of the logistic curve in the real range of the data.

However, fitting this general form directly can lead to instability, particularly when the number of experimental time points available is limited³³. To mitigate this and provide a simpler and easier-to-understand perspective, we fit the data using a simplified form of the ODE (after appropriate translation and normalization of the data; see Supplementary Note 1, Propositions 6–9, and Supplementary Fig. 2):

where t_start represents the initial time point of the experimental time course, and \({y}_{{{\rm{start}}}}^{ * }\) represents the normalized gene expression level or chromatin signal of the genomic locus at t_start.

We first solve Eq. 2 analytically to obtain an explicit formulation for y* (see Supplementary Note 1, Proposition 1, for the derivation of the analytical formula):

Then, we numerically estimate the two kinetic parameters k* and b* that best fit the data. Specifically, our proposed chronODE package optimizes this solution by repurposing a generic initial value problem solver and constraining the numerical solution to the logistic space, so that parameters k* and b* (and their corresponding k and b values back-transformed to the original range) can be easily interpreted. In this context, b* is a positive number and corresponds to the saturation level of the normalized logistic curve, i.e., the horizontal asymptote that approximates the maximum predicted level of y*, while k* represents the growth/decay rate constant, which dictates how fast the signal ramps up or slows down over time. In the case of increasing (i.e., activated) profiles, k* > 0, whereas decreasing (i.e., repressed) profiles have k* < 0 (Fig. 1).

We describe the logistic curve as consisting of two biologically meaningful phases. First, the acceleration phase (where \(| {y}^{ *\prime }| \) progressively increases over time) models the cooperative aspect of epigenetic and transcriptional processes, which tend to occur more frequently and easily after the initial state^{10,11,12,13,14,16}. Second, the deceleration phase (where \(| {y}^{ *\prime }|\) progressively decreases) models the saturating aspect of these processes (Fig. 1). We define t_switch as the time point that marks the transition between these two phases as follows:

Equation (1) not only models the kinetics of genomic signals while accounting for the underlying biophysical and biochemical constraints, but it also provides the flexibility to capture diverse patterns beyond a standard logistic curve. In particular, it allows for variability in the behavior of genes or regulatory elements, accommodating patterns such as exponential, linear, or logarithmic-like trajectories, which correspond to the accelerating, switching (inflection point), and decelerating phases of the logistic curve, respectively. In essence, it does not assume that all signals strictly follow a full logistic profile.

We note that this approach is well-suited for modeling genomic signals that exhibit monotonic increases or decreases. While these monotonic fits describe most of the observed transcriptional programs, a small fraction of genes may display peak-like expression profiles over time, such as those involved in circadian clocks, cell cycle regulation, or stress responses^34,35,36. Thus, we further extended our methodology by introducing a piecewise fitting approach that can accommodate such profiles by combining two logistic curves (or two portions of them). Although alternative functions (e.g., Gaussian curves) could be employed for these cases, our choice to use logistic functions maintains consistency within the framework and allows for comparable kinetic parameters across both monotonic and piecewise fits. Both approaches are implemented as standalone pipelines, providing a scalable solution for genome-wide applications at both bulk and sc levels.

Kinetic profiling of the developing mouse brain transcriptome reveals a balance between rate and saturation level

We used chronODE to investigate the kinetics of gene expression regulation during mouse brain development. We first analyzed time-series RNA-seq maps generated for three regions of the mouse brain (forebrain, midbrain, and hindbrain) across eight developmental time points, from post-conception (PC) day 10.5 to the first postnatal day (PN) (Supplementary Table 1)²⁷. On average across the three regions, we identified  >12,000 genes differentially expressed (DE) over time (Supplementary Fig. 3A). We applied chronODE on the set of DE genes (defining, for a given gene, t_start = 10.5 and \({y}_{{{\rm{start}}}}^{ * }\) = normalized expression level of the gene at time point 10.5; see “Methods”), and identified an average of 87% monotonic fits, suggesting that the logistic ODE is appropriate for modeling the majority of time-series patterns in our dataset (Supplementary Fig. 3B). Conversely, a small fraction of genes (average of 8%) displaying peak-like patterns were captured using piecewise fits (Supplementary Fig. 3C–F). The remaining  ~5% of genes exhibited complex or variable patterns that were not well-suited for fitting with the current models.

When analyzing the distribution of k and b values across genes and brain regions, we found that the monotonic genes are clustered into three major quadrants based on their k & b combinations: high k & low b (Q1), low k & low b (Q2), and low k & high b (Q3) (Fig. 2A). This distribution exhibits an L-shaped pattern (Pearson’s r correlation = −0.44, p value < 2.2e-16), where high values of both k and b are not observed together (i.e., Q4 is absent). A closer examination revealed that Q1 genes (27%) change expression faster, with expression profiles following the full logistic curve (Fig. 2B). In contrast, Q2 genes (59%) and Q3 genes (14%) align with only portions of the logistic curve and are characterized by slower kinetics. Specifically, Q2 genes approach saturation, while Q3 genes remain far from saturation throughout the observed time frame.

Within these slowly saturating Q3 genes, those repressed during brain development are associated with nucleotide biosynthesis and amino acid metabolism—metabolic processes essential for early cell proliferation but downregulated later on during lineage commitment^37,38 (Supplementary Fig. 4A). This result is consistent with these genes being already far from saturation by day 10.5 (Fig. 2B). In contrast, Q3 activated genes are enriched for functions related to oxidative phosphorylation and mitochondrial ATP synthesis (Supplementary Fig. 4B). These key components of aerobic respiration become increasingly prominent at later developmental stages, particularly during the peri- and post-natal periods³⁹, aligning with the slower activation and delayed saturation of Q3 activated genes (Fig. 2B).

We performed a similar analysis for genes with peak-like expression profiles, and found that, also in this case, only a small fraction (~8%) corresponds to Q4 genes (Supplementary Fig. 4C–F). Altogether, this suggests that genes cannot undergo rapid expression changes while simultaneously achieving high saturation. Consequently, the expression level of rapidly saturating genes is fundamentally constrained (Q1 genes).

Although most genes follow monotonic expression patterns, our results indicate that not all of them transition from full inactivation to full activation (or vice versa) at the same pace (Fig. 2B), and that these differences may correlate with the timing of when their functions are required. To systematically capture the temporal progression of gene expression kinetics, we developed a classification scheme that directly ties the gene’s kinetic profile, as dictated by the k and b parameters, to the specific time window defined by t_start and t_end (in our case, days 10.5 and PN, respectively). This classification identifies the switching time (t_switch) of the gene—the point at which its expression shifts from accelerating to decelerating (see Eq. 4)—and determines whether t_switch occurs before (decelerator), within (switcher), or after (accelerator) the time window monitored by the time-series study (Fig. 2C, Supplementary Fig. 5A, and Supplementary Data 1). This allows us to describe the full signal curve in the context of the time window [t_start, t_end] (Supplementary Fig. 5B, C). Switcher profiles were the most abundant among both activated and repressed genes, while only a small fraction of genes were classified as accelerators (Fig. 2D). Moreover, although gene activation and repression were largely conserved, with nearly 100% of the genes showing either always an increasing or always a decreasing trend across the three brain regions, their kinetic patterns varied significantly. In fact, 47% and 40% of activated and repressed genes, respectively, were assigned to a different kinetic category in at least two regions (Fig. 2E).

Overall, chronODE allowed us to decipher the kinetics of gene expression during brain development at the resolution of individual genes, highlighting how changes in gene expression result from a balance between the rate of expression and the extent of saturation.

Cell type emergence and gene essentiality shape gene expression kinetics

After mathematically quantifying the kinetics of gene expression during brain development, we next sought to uncover their biological meaning. Our classification indicates that the kinetic patterns of certain genes differ substantially between the forebrain, midbrain, and hindbrain (Fig. 2E). We hypothesized that these variations could be influenced by the cellular composition of each brain region and the changes in this composition over time. In fact, each region is composed of specific cell types, with specific characteristics, that emerge at different developmental stages, potentially shaping the kinetics of gene expression⁴⁰. To test this hypothesis, we extended our analysis from tissue-level data to sc resolution. Specifically, we applied our chronODE framework to a recently published time-series scRNA-seq atlas of the developing mouse embryo²⁹. This dataset includes 41 brain-related cell types that emerge at various time points between embryonic days E8.5 and E14.25. For our subsequent analyses we focused specifically on activated genes, since they were detected across all cell types, in contrast to repressed genes which were mostly restricted to cell types appearing between E8.5 and E9 (Supplementary Fig. 6A).

Genes with distinct kinetic behaviors, as defined by the three quadrants in Fig. 2A, were observed also across most individual cell types, consistent with the diversity of gene expression kinetics identified in the bulk tissue analysis (Supplementary Fig. 6B). However, the L-shaped distribution of these three kinetic subpopulations changed markedly over time as new cell types emerged. Stratifying cell types into five temporal groups revealed that genes become progressively constrained to a smaller range of k–b combinations at later developmental time points (Fig. 3A). In other words, late cell types, such as amacrine precursors and cholinergic amacrine cells, show more constrained kinetic patterns of gene activation, compared to early cell types (e.g., neural progenitors and motor neurons) (Supplementary Fig. 6B). Moreover, when taking a closer look at their temporal progression, we found that genes activated in early-appearing cell types predominantly followed decelerating profiles and were already close to their saturation level b by day E8.5 (Fig. 3B). In contrast, genes activated in late-emerging cell types displayed mostly switching profiles and, to a lesser extent, accelerating profiles. These findings confirmed our initial hypothesis, suggesting that the timing of cell type appearance seems to play a critical role in shaping gene expression kinetics.

Still, among genes activated in the early cell types, we observed a fraction following slower, accelerating profiles (average of 11% in E8-E9 cell types; Fig. 3B). This suggests that even though these cells are present from early stages of development, they selectively express certain genes at slower rates. This observation prompted us to investigate additional factors that might contribute to the regulation of gene expression kinetics. In particular, we hypothesized that a cell may tend to prioritize the expression of genes that are essential for its survival. To this end, we integrated in our sc analyses a recently published catalog of essential genes⁴¹. We observed that, especially in the earliest cell types (E8-E9 appearance days), essential genes tend to ramp up faster (higher k values) and are closer to their saturation levels (i.e., higher proportion of decelerator and switcher profiles) compared to non-essential genes (Fig. 3C, D). Thus, our kinetic analysis provides an alternative perspective: the essentiality of a gene for cell viability likely influences its expression kinetics, determining how rapidly it is expressed.

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

Besides estimating gene expression kinetics, chronODE is also suitable for studying the kinetics of chromatin activity at regulatory elements in the genome, providing insights into gene regulatory mechanisms. To this end, we also analyzed the kinetics of chromatin accessibility at candidate cis-regulatory elements (cCREs) from the ENCODE project⁴². Specifically, we applied chronODE to time-series DNase-seq and ATAC-seq maps²⁸ available for the same time points and brain regions as the bulk RNA-seq data (Supplementary Tables 2 and 3). We identified cCREs with both increasing (k > 0) and decreasing (k < 0) trends, and classified their kinetics in each of the three regions (Supplementary Data 1). We then investigated how activated and repressed cCREs relate to the expression of target genes over time. To achieve this, we employed a simple proximity criterion to pair mouse cCREs with target genes by linking each cCRE to its nearest gene based on linear distance^43,44. Although this is a simplified approach, we note that our framework could also accommodate more complex linking methods, such as those recently applied by ENCODE in the context of the human genome⁴⁵. Our analysis revealed that genes are not only linked to multiple cCREs, as previously reported⁴⁴, but that these cCREs also exhibit distinct regulatory roles. Specifically, on average 50% of genes are regulated by both activated (k > 0) and repressed (k < 0) cCREs (Supplementary Fig. 7A). We refer to these genes, which display mixed cCRE trends, as poly-pattern (poly-p.) genes (Fig. 4A). In contrast, mono-pattern (mono-p.) genes are exclusively associated with a single type of cCRE trend—either all activated or all repressed.

Across all three brain regions, genes with poly-pattern regulation were linked to more distal cCREs and exhibited larger expression changes over time compared to mono-pattern genes (Fig. 4B, C). This suggests that poly-pattern genes play a more significant role in brain development, as they are enriched in neural and brain-specific functions, including neurogenesis, trans-synaptic signaling, and sensory system development (Fig. 4D). In contrast, mono-pattern genes were predominantly associated with housekeeping functions such as gene expression and compound biosynthesis. A central question is what drives the highly dynamic behavior and specialized functions of poly-pattern genes: is it the number of regulatory elements or the diversity of their regulatory trends? To address this, we disentangled the contributions of cCRE count and trend diversity (mono- vs. poly-pattern regulation). Our analysis demonstrated that both the number of cCREs and the diversity of their trends significantly correlate with a gene’s dynamic expression profile and its involvement in brain-specific functions (Supplementary Fig. 7B, C). This finding expands on our conclusion that gene expression is fundamentally a constrained process by having a balanced regulation involving both activation and repression mechanisms.

Overall, these results suggest that the precise expression of genes important for brain development may be governed by a more complex network involving regulatory elements with mixed trends. In contrast, those genes associated with non-brain-specific functions may rely on a simpler regulome, potentially reflecting a differential degree of control based on the biological significance and impact of these genes.

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

The insights gained from our previous analyses improve our understanding of how gene expression evolves over time, particularly in terms of expression kinetics and the impact of epigenetic changes at regulatory elements. This knowledge can be applied to develop temporal models that predict gene expression at future time points, accounting for the additive effect of multiple cCREs associated with a given gene⁴⁴. Beyond prediction, incorporating chromatin regulation into these models can also provide insights into the complex interplay between transcriptome and epigenome.

Towards this goal, the first step in our analysis was to examine the correlation, over time, between gene expression and chromatin accessibility at associated cCREs. We observed a bimodal distribution, with 66% of the gene-cCRE pairs showing positive correlation and 34% showing negative correlation (Supplementary Fig. 8A). In this context, we make the simple assumption that when a cCRE is open, TFs can access the gene, facilitating expression. In contrast, a closed cCRE restricts TF access and potentially reduces gene activity. Based on this assumption, the positively correlated cCRE-gene pairs were interpreted as cases where the cCRE acts as an enhancer of gene expression, since the accessibility of the cCRE aligns with the direction of gene expression. In contrast, the negatively correlated cCRE-gene pairs were interpreted as cases where the cCRE exerts a silencer effect on gene expression. This information is crucial for determining the direction of gene expression changes (activation or repression).

On the other hand, we found that the kinetics of gene expression were correlated with the diversity of cCRE regulation associated with the gene (mono-p. vs. poly-p.; Supplementary Fig. 8B). Specifically, activated genes with decelerator profiles were more likely to undergo poly-pattern regulation compared to switchers and accelerators, while the opposite was observed for repressed genes. This suggests that the type of epigenetic regulation a gene undergoes plays a key role in shaping the kinetics of its expression profile, making it important to incorporate this information into the predictive model.

Building on these findings and utilizing our many-to-one linking schema—where multiple cCREs are linked to a given gene over time—we designed a model architecture (Fig. 5A) capturing four types of gene regulatory mechanisms that differ in the direction (enhancer/silencer) and diversity (mono-p./poly-p.) of cCRE regulation as follows: (1) mono-pattern genes subject to enhancer regulation (enhancer mono-p.), where all cCRE signals linked to a gene are positively correlated with its expression; (2) mono-pattern genes subject to silencer regulation (silencer mono-p.), where all cCRE signals linked to the gene are negatively correlated with its expression; (3) poly-pattern genes subject to enhancer regulation, where the majority of cCRE signals are positively correlated with gene expression, while a minor fraction exhibit silencer-like effects (enhancer poly-p.); and (4) poly-pattern genes subject to silencer regulation, where most cCRE signals are negatively correlated with gene expression, with a smaller proportion displaying enhancer-like effects (silencer poly-p.). This approach enabled us to train our model towards predicting both the direction and kinetic patterns of gene expression over time. Each subset of genes was used to train the model independently, with an 80:20 split for training and testing, respectively.

Our modeling approach predicts gene expression at each time point through regression, using a bidirectional RNN-based architecture, which is well-suited for time-series data. Given a multicCRE-gene pair i, we define the complete cCRE time-series input signal associated with a given gene as the matrix C_i, such that:

At each time point t (t ∈ {1, 2, …, 8}) we define c_i,t as the vector of m cCRE signals:

C_i is the input used to predict g_i, i.e., the gene expression profile over time:

We use a biRNN to model gene expression g_i,t at a particular time point t, as a function of the cCRE chromatin accessibility signal vector c_i,t. In fact, capturing the interplay between chromatin and gene expression is not a trivial task, since multiple scenarios are possible. For instance, changes in chromatin accessibility at cCREs upstream of the promoter can potentially influence gene expression, for example, by enabling TF binding or facilitating Pol II release. Conversely, transcriptional activity—through Pol II elongation—can also induce changes in chromatin accessibility at both upstream and downstream cCREs. The biRNN allows us to cover this wide range of scenarios by leveraging chromatin signals at both past and future time points to make predictions of gene expression. Specifically, we define the following relation:

where \({{\mathcal{D}}}\) represents a dense layer that takes as input the output of the biRNN.

This modeling approach effectively handles the highly dimensional nature of the omics data and enables accurate regression of gene expression levels exhibiting diverse kinetic patterns over time. We note that although here we focused exclusively on chromatin accessibility as input data, our modeling architecture is versatile and can potentially incorporate additional cCRE chromatin features (such as histone post-translational modifications) at each time point t in the form of a tensor.

To evaluate the performance of the model, we first computed cross-gene correlations between true and predicted expression values. Our predictions on the test set showed a strong positive correlation (r > 0.85) in all four types of regulatory mechanisms (Fig. 5B), higher compared to previous reports of chromatin-gene expression predictions⁴⁶. We also evaluated the model performance by computing the distribution of mean squared error (MSE) and cross-timepoint correlation between true and predicted gene expression values (Supplementary Fig. 8C). In particular, we achieved mean cross-timepoint correlations ranging between 0.57 and 0.90, higher than equivalent cross-cell type correlations reported previously⁴⁶. To further showcase the predictive capability of our model, we also included some examples of true vs. predicted gene expression time-series profiles (Fig. 5C). We note that the usage of the biRNN is especially effective for predicting the expression of poly-pattern genes, which are particularly challenging to model due to the mixed contributions of enhancer and silencer cCRE signals (Supplementary Fig. 8D). The biRNN leverages the more comprehensive learning context provided by both past and future time points to address this complexity. Additionally, incorporating a dense layer enhances the model’s capacity to learn complex kinetic relationships between cCREs and their associated genes.

Finally, we investigated how the distance between cCREs and their target genes influences gene expression predictions using the SHapley Additive exPlanations (SHAP) method⁴⁷. Specifically, we computed the average importance of each cCRE across all the samples in the test set. Our analysis revealed that closer cCREs have a greater impact on predicting gene expression, underscoring their biological importance in gene regulation, consistent with previous findings⁴³. Notably, these proximal cCREs play a dominant role in predictions at the first and last time points, whereas middle-to-distant cCREs contribute more to the predictions at intermediate time points, thus directly influencing the kinetic patterns of the gene expression signals (Fig. 5D). This finding indicates that multiple cCREs provide the flexibility needed to capture complex kinetic patterns of gene expression. Altogether, this suggests that, during brain development, gene expression kinetics are shaped by the cumulative effects of numerous cCREs, especially the more distal ones.