Gaussian Processes for Time Series with Lead-Lag Effects with application to biology data

The lead-lag effect, a widespread phenomenon where changes in one time series (the leading series) influence another time series (the influenced follower series) after a certain delay. This phenomenon is prevalent in a variety of fields, such as climate science (De Luca and Pizzolante, 2021), healthcare (Runge et al., 2019; Zhu and Gallego, 2021; Feng et al., 2021), economics (Wang et al., 2017; Skoura, 2019), and finance (Hoffmann et al., 2013; Bacry et al., 2013; Da Fonseca and Zaatour, 2017; Ito and Sakemoto, 2020). For example, Harzallah and Sadourny (1997) explores the lead-lag relationship between the Indian summer monsoon and various climate variables, such as sea surface temperature, snow cover, and geopotential height. However, this dynamic, rooted in causal regulatory relationships, is notably observed in biological processes, which serves as the primary inspiration for the proposed method. The lead-lag effect manifests itself in various processes, including development (Gerstein et al., 2010; Ding and Bar-Joseph, 2020; Strober et al., 2019), disease-associated genetic variation (Krijger and De Laat, 2016), and responses to various interventions and treatments (Faryabi et al., 2008; Lu et al., 2021).

To provide a clearer context for this concept, we consider temporal-omics datasets, particularly those used in the identification of potential enhancer-promoter pairs in various biological processes (Whalen et al., 2016; Fulco et al., 2019; Schoenfelder and Fraser, 2019; Moore et al., 2020). Enhancers are short cis-regulatory DNA sequences that are distal to the genes they regulate, as opposed to the promoter sequence that is adjacent to a particular gene. Since enhancers and promoters coordinate to generate time- and context-specific gene expression programs, that is, the amount of RNA produced at what time, we hypothesize that the dynamic patterns of enhancer regulation are typically followed by corresponding gene program patterns, albeit with a certain time lag. Approximately a million candidate human enhancers have been identified (Moore et al., 2020), though the intricate dynamics of how and which enhancers interact with promoters in specific processes is still an area of active and ongoing research.

Addressing this challenge, our study focuses on the multi-omics time series data presented by Reed et al. (2022). This dataset includes measurements of enhancer activity, and promoter transcription, and 3D chromatin structure at seven time points, from 0 to 6 hours (a later time point is excluded from this analysis). These data reveal complex dynamic patterns, as depicted in Figure 1. Our goal is to identify functional human enhancer-promoter pairs from over 20,000 candidate pairs, estimate the time lag, and quantify the extent of the lead-lag effect between enhancer activity and gene expression over time. The biological motivation for this goal is to uncover a relationship in genomics time series for a subset of enhancer-gene pairs, in particular those in which the enhancer state changes after cell activation and causes a subsequent change in expression of a target gene. Measuring the lag between enhancer activity and gene expression over time for these candidate enhancer-promoter pairs helps in understanding the sequence of regulatory events set off by activation, and how the time scale of gene regulation may vary across the genome.

In the context of our motivating application involving multi-omics time series data, and similar scenarios across various fields, there are key challenges that need careful consideration. Firstly, a common issue is the scarcity of time-point observations. This limitation is compounded by the irregular intervals between successive measurements, often due to the high cost of experimental design or sensor failures. Additionally, in cases where there is an abundance of time-point data, the focus may shift to pinpointing local time delays within smaller segments of the time series. These factors underscore the importance of accurately quantifying the uncertainty associated with the extent of the lead-lag relationship, particularly when dealing with limited temporal datasets.

Existing methods for aligning time-lagged time series to explore the magnitude of lead-lag effects are numerous and varied. This exploration can extend to tasks such as ranking leaders or followers among multiple time series (Wu et al., 2010) and performing network clustering of time series based on the pairwise lead-lag relationships (Bennett et al., 2022). For example, Wu et al. (2010) ranks the leaders of the streams to the sea surface temperature. Widely used models such as time-lagged cross-correlation (TLCC, Shen (2015)) by comparing the cross-correlation coefficients, dynamic time warping (DTW, Berndt and Clifford (1994)), and its differentiable approximation, soft-DTW (Cuturi and Blondel, 2017), and soft-DTW divergence (Blondel et al., 2021) have proven useful by comparing distance loss among multiple pairwise time series after alignment (De Luca and Pizzolante, 2021). TLCC, for instance, operates by sliding one signal over another and identifying the shift position with the highest correlation. However, these methods often require the same time intervals for each time series, which makes them less effective for irregular time series, potentially leading to distortions of the true underlying relationships.

Beyond quantifying the lead-lag effect between pairwise time series, there is often a desire to estimate time lag considering the irregularities in time series data, one approach is to first impute values at regular time points, then apply methods like TLCC to calculate the time lag. Common techniques for this separate modeling of each time series include using splines (Rehfeld et al., 2011) or separate Gaussian processes (SGP, Li et al. (2021)). However, such methods do not leverage shared information across time series, potentially leading to sub-optimal results. Hierarchical GP models such as Magma (Leroy et al., 2022, 2023) integrate shared information across multiple time series but are primarily designed for prediction and clustering individual time series. This is in contrast to our primary objective, which is to cluster pairwise time series by a measure of similarity. More direct approaches for irregular time series include the use of the negative and positive lead-lag estimator (NAPLES, Ito and Nakagawa (2020)) and Lead-Lag method (Hoffmann et al., 2013), which utilize the Hayashi–Yoshida covariance estimator to estimate a fixed time-lag. Alternatively, multi-task GPs through latent variable models introduce additional variables to account for temporal discrepancies such as GPLVM (Dürichen et al., 2014), which may limit flexibility due to the presumed independence between tasks and temporal sequences (also known as separable kernels). Moreover, more complex misalignment techniques (Kazlauskaite et al., 2019; Mikheeva et al., 2022) might not perform well with limited time observations for biology data where simpler linear mismatch assumptions suffice, and they can face operational challenges like Cholesky decomposition failures. Table 1 compares the aforementioned methods, highlighting their suitability for limited or irregular time points, their ability to quantify the magnitude of lead-lag effects, and estimate time lags:

In response to these challenges and limitations of existing methods, we propose the Gaussian Process for lead-lag time series (GPlag) model. GPlag introduces a novel class of GP kernels for time series with lead-lag effects, featuring interpretable parameters that estimate both the time lag and the degree of lead-lag effects. These parameters are crucial for tasks such as identifying associated factors, ranking or clustering time series based on lead-lag relationships. Our model also provides theoretical justification for the identifiability and interpretability of these key parameters, and integrates a Bayesian framework for greater flexibility in handling time series with constraints like limited observations or irregular time gaps.

To support our interpretation to $A\coloneqq(a_{ll^{\prime}})\in\mathbb{R}^{L\times L}$ and $S\coloneqq(s_{l})\in\mathbb{R}^{L}$ in Section 2, we discuss the identifiability, i.e., different parameters lead to different models. Note that not all parameters of commonly used GP models are identifiable, which means that different parameter values can lead to the same model. For example, in the classical Matérn covariance function, the spatial variance $\sigma^{2}$ and lengthscale parameter $b$ are not identifiable, meaning that it is not possible to distinguish between different values of these parameters based on the data alone (Stein, 1999; Zhang, 2004; Kaufman and Shaby, 2013; Li et al., 2023; Li, 2022). This can be a limitation of GP models and it is important to consider when interpreting the results. In the case of GPlag, if the key parameters $A$ and $S$ are not identifiable, i.e., the GPlag determined by $(A,S)$ is equivalent to the one determined by $(\widetilde{A},\widetilde{S})$, then it is unconvincing to interpret $A$ and $S$ as times series dissimilarity and time lag. Without identifiability, $A$ and $S$ do not admit any clear interpretation and the results of the model would be hard to interpret.

In this section, we prove that the parameters $A$ and $S$ in GPlag model are identifiable, allowing them to be used to as measures of the lead lag effect between time series (measured by $A$) and the time lag (measured by $S$). Here, we assume all parameters are positive, finite, real numbers, and the domain is fixed (also known as infill domain). The fixed domain assumption is reasonable for multi-omics time series, as the dynamics of gene regulation in such experiments occur during a fixed window of time. Investigators will decide, based on available budget, the time points needed to capture the dynamics during, say activation, or differentiation. If allocated more budget, they will then fill in the gaps between existing observations, to have better resolution over the same time period. We begin by introducing the necessary definitions.

If $K\equiv\widetilde{K}$, then $K$ can never be correctly distinguished from $\widetilde{K}$ regardless how dense the observed time points are (Zhang, 2004). Focusing on a parametric family of covariance functions $K_{\theta}$, if there exists $\theta_{1}\neq\theta_{2}$ such that $K_{\theta_{1}}\equiv K_{\theta_{2}}$, then any estimator of $\theta$ based on $n$ observations $\{t_{i},y_{i}\}_{i=1}^{n}$, denoted by $\widehat{\theta}_{n}$, cannot be weakly consistent (Dudley, 1989).

As a result, it suffices to show the following necessary and sufficient conditions for two GPlags being equivalent, which suggests identifiable parameters for further interpretation:

The proof is given in Appendix A.2. Although $\sigma^{2}$ and $b$ in GPlag are not identifiable, the proposed interpretation of the model does not rely on these two parameters. Instead, the parameter $A$ is interpreted as a measure of dissimilarity between time series under a time lag represented by the parameter $S$.

It is important to note that GPs are commonly used for prediction and regression tasks. For this purpose, even though some kernel parameters such as $\sigma^{2}$ and $b$ are not identifiable in these models, the prediction performance remains asymptotically optimal when the kernel parameters are mis-specified (Kaufman and Shaby, 2013). This is because regression is a distinct problem from parameter inference in GPs.

It is also worth noting that the identifiability of kernel parameters in GPs depends on the dimension of the domain. For time series, the domain is a 1-dimensional interval in $\mathbb{R}^{1}$. However, our theorem holds for 2-dimensional and 3-dimensional domains. In dimensions greater than $3$, it becomes a challenging problem and we treat it as a future work since our focus of this work is one-dimensional time series.

The main focus of our experiments was to illustrate that GPlag could accurately estimate both the time lag and the extent of lead-lag effects. This capability aided in the ranking or clustering of multiple time series that existed in a lead-lag relationship with a target time series. We demonstrated the theoretical properties of GPlag using synthetic data, as these properties were best illustrated under a fully controlled environment. All code and additional implementation details were available in the Section B.1.

Baseline models 1) For estimating the time lag, we compared with Lead-Lag (Hoffmann et al., 2013). 2) By utilizing the GP framework, we were able to predict values for two time series at unobserved time points. Then we evaluated the prediction error against three baselines: natural cubic splines with five breakpoints, SGP and Magma. 3) To evaluate the effectiveness of the estimate of lead-lag effects, we employed four baseline methods: TLCC, DTW, soft-DTW, and soft-DTW divergence. TLCC and DTW could output aligned time series, allowing us to calculate correlation coefficients from the aligned data. On the other hand, soft-DTW and soft-DTW divergence measured the amount of warping needed to align two time series and provided a distance loss metric. We used both correlation coefficients and distance loss as measures of dissimilarity between the time series.

Parameter estimation and inference. To validate the identifiability of our model, we simulated a pair of time series from GPlag with three kernels in Equations 2, 3 and 4, with 100 replicates. We set the mean to zero for both series and specified parameters as $b=0.3,\leavevmode\nobreak\ a=1,\leavevmode\nobreak\ s=2,\leavevmode\nobreak\ \sigma^{2}=4$ for each kernel. For LMat, we set the smoothness parameter as $\nu=3/2$. The number of measurements within each time series varied across $\{20,50,100,200\}$ and the time points were generated following $t_{i}=i+\epsilon$ and $t^{\prime}_{i}=t_{i}+s$, where $i\in\{1,2,\dots,n\}$ and $\epsilon\sim\text{Unif}(-\frac{1}{4},\frac{1}{4})$ represents noise.

Our results showed that the MLE of $a$ and $s$ were very close to their truth among all three kernels when $n=100$. This was consistent with the theoretical analysis in Theorem 3 (Figure 2 A, B) that GPlag could accurately recover the true lead-lag effects. Compared to Lead-Lag, GPlag exhibited similar or smaller variance, specifically when the kernel was Matérn with $\nu=3/2$. Additionally, we performed simulations with a varying number of time points and found that the variance of the MLEs of $a$ and $s$ decreased as the sample size increased, which empirically confirmed the consistency of the MLEs (Figure 2 C). Plots for other parameters of secondary interest, including $\sigma^{2}$, $b$, and $\sigma^{2}b^{2\nu}$, are provided in Figure S1. We also performed simulations with different settings of $b=0.3,5,10$, and found that the consistency of the MLEs for $a$ and $s$ was maintained across these different values of $b$ (See Figure S2 for more details).

Prediction. We assessed the accuracy of GPlag’s predictions by generating data from two different processes. In the first, we used the kernel LExp in Equation 3 with parameters $b=1,\leavevmode\nobreak\ a=0.3,\leavevmode\nobreak\ s=2$, and $\sigma^{2}=4$. In the second, we generated pairwise time series from linear functions with non-Gaussian noise (t-distributed): $y=2t+3+5\epsilon_{1}$ and $y=2(t-20)+3+5\epsilon_{2}$ where $t=0,1,2,\ldots,100$. The noise term, $\epsilon_{1},\epsilon_{2}$, were drawn from a t-distribution with 5 degrees of freedom. We then randomly selected $50\%$ of the data as training data for model fitting and used the remaining $50\%$ as testing data for prediction in both settings. We applied the best linear unbiased prediction (BLUP, see Section B.1 for more details) and used the mean squared error (MSE) on the testing set as the metric to evaluate the prediction accuracy.

Among all methods, GPlag achieved the highest log-likelihood and smallest MSE (Figure 3 A-B, Figure S3). The advantage of GPlag over SGP and splines was that it modeled two time series together which allowed for sharing of information between time series through $a$. When $a\rightarrow\infty$, it was equivalent to SGP. While a small estimated $a$ indicated that two time series were similar, and the covariance function then allowed for more information to be shared between them. Although Magma adopted a hierarchical GP to specifically model the mean function for irregular data, it did not consider time lag information across time series.

Clustering or ranking time series based on dissimilarity to target time series.

We evaluated the relationship between the magnitude of parameter $a$ and the degree of lead-lag effects (the level of synchrony with a designated target signal) to assess the effectiveness of GPlag. We simulated nine time series from the function $f_{k}(t)=\frac{\arctan(k(t+s))}{\arctan(k)}$ for every pair of $k$ and $s$, where $k\in\{0.01,1,10\}$ and $s\in\{0,0.5,1\}$. Among them, the time series simulated from the function of $k=0.01,\leavevmode\nobreak\ s=0$ was treated as the target time series, which was close to a linear function (Figure 3 C, ts1). Each time series had $50$ times points ranging from $[-2,2]$. Here, $k$ determined the degree of distortion, where a smaller $k$ (e.g. 0.001) corresponded to a curve close to the identity transformation, and a larger $k$ represented the cases with more upward/downward distortion. By comparing different columns in Figure 3 C, we could see that the degree of dissimilarity between the target time series and other time series was controlled by $k$, whereas $s$ controlled the time lag along the $t$-axis.

Nine pairs of time series were fitted by GPlag with the LExp kernel. We observed that three pairs of time series in each column with the same $k$ were grouped together (Figure 3 D), and $a$ increased as $k$ increased. This implied that GPlag could accurately rank the degree of lead-lag effects among multiple pairs of time series based on the parameter $a$, through precise estimation of the time lag parameter $s$. A larger estimate of $a$ indicated a stronger dissimilarity between pairs of time series. To further quantify the clustering performance of lead-lag effects among the nine pairs of time series, we used the Adjusted Rand Index (ARI, Hubert and Arabie (1985)) and Normalized Mutual Information (NMI, Strehl and Ghosh (2002)) as evaluation metrics. To define a cluster from the pairwise dissimilarity parameter (correlation coefficients for TLCC and distance loss for DTWs), we used K-means (Shannon, 1948) with 3 clusters. As the focus of this work was on determining pairwise dissimilarity rather than the clustering algorithm, we opted for a simple clustering method. GPlag achieved the highest ARI and NMI in clustering the lead-lag effects among 9 pairwise time series (Table 1, column 2-3). We also tested GPlag using the LRBF and LMat kernels, and found that ARI and NMI scores remained at 1, consistent with the LExp results. The inferred $a$ values are shown in Figure S4, further indicating that the results are not sensitive to the choice of kernel.

Extension to three time series. To demonstrate the applicability of GPlag to more than two time series, we conducted simulation studies using three time series generated from the LExp kernel. In this scenario, all constraints were linear, allowing us to utilize optimizers that support linear constraints. Similar to the pairwise simulation setting, we set the mean to zero for all series, and the time points were generated according to the equations $t_{i}=i+\epsilon$ and $t^{\prime}_{i}=t_{i}+s_{i}$, where $i\in\{1,2,\dots,50\}$ and $\epsilon\sim\text{Unif}(-\frac{1}{4},\frac{1}{4})$ represented noise. The specified parameter values were $b=0.3,a_{12}=a_{13}=a_{23}=1,s_{2}=2,s_{3}=4,\sigma^{2}=4$. The maximum likelihood estimate (MLE) of the parameters closely approximated the true values (represented by the red dashed line), providing further validation of the theoretical properties. Additionally, the proposed method achieved the lowest mean squared error (MSE) in this particular setup (see Figure 4).

In this section, we applied GPlag to a time-resolved multi-omics dataset involving enhancers and genes in human cells activated by interferon gamma from a resting state (Reed et al., 2022). Originally, the dataset comprised over 20,000 pairs. Each pair consisted of one time series for each enhancer and likewise one time series for each gene, at seven irregularly-spaced time points (0, 0.5, 1, 1.5, 2, 4, and 6 hours after a perturbation). The time series captured gene expression measurements for genes and histone tail acetylation measurements (H3K27ac), a measure of activity for the enhancers, both approximately variance stabilized. It is widely believed that the alterations in acetylation occur prior to the changes in expression of the regulated genes (Reed et al., 2022). We focused only on pairs that became functional after stimulation, by applying a filter on the total range of stabilized values across the timescale, yielding 4,776 dynamic pairs of enhancer-gene time series. Subsequently, the objective of our analysis, after estimating the model parameters of GPlag for these 4,776 pairs, was to identify the subset of all possible enhancer-gene pairs that represent functional enhancer-gene pairs, i.e. the enhancer responds to cell activation and regulates the expression of the gene, after some time lag. Many types of time profiles can be captured by such a time series experiment, e.g. increasing across the range, decreasing, or increasing then decreasing, etc. Particular functional enhancer-gene pairs will have subtle differences in the dynamics of the coordinated change in activity and expression, compared to other functional pairs, which we hope to capture with our GPlag model.

To assess the performance of GPlag in ranking enhancer-gene connections, we split the pairs into two groups according to the parameter $a$. The top 20% pairs were designated as the candidate enhancer-gene group, with the remaining pairs forming the negative group. To compare against baselines, we split all pairs into two groups using their respective metric of similarity and a 20% cutoff. Subsequently, we assessed the performance from three perspectives.

Firstly, we compared the Hi-C (High-throughput Chromosome Conformation Capture, Lieberman-Aiden et al. (2009)) counts between the candidate enhancer-gene group and the negative group. Since Hi-C counts measured the frequency of physical interactions between enhancers and genes in three-dimensional chromatin structure, we expected the candidate enhancer-gene group to have higher Hi-C counts, while the Hi-C information was not used in our ranking of pairs. Our results showed candidate enhancer-gene pairs exhibited significantly higher Hi-C counts compared to the other pairs (Table 2, column 4-5), with GPlag producing the smallest $p$-value (Wilcoxon rank-sum test, $p=0.019$).

Secondly, we analyzed the enrichment of looped enhancer-gene pairs between the two groups. A looped enhancer-gene pair is one where the enhancer element is physically brought into close proximity with the regulated gene, enabling the enhancer’s activity to influence the expression of the associated gene. Looping in this dataset was determined using Hi-C counts, but the determination of loops uses spatial pattern detection and statistical models to determine if the signal rises to the level of a “loop”. We observed that candidate enhancer-gene pairs identified by GPlag had a significantly higher likelihood of being “looped” enhancer-gene pairs compared to the matched pairs from the other group, selected based on genomic distance (Davis et al., 2022) ($\chi^{2}$ test, $p=0.05$). However, no significant difference in enrichment was found between the two groups of pairs identified using baseline methods.

The lack of statistical significance observed in DTW-based methods may be due to the limited information available for alignment because of the small number of observations. Additionally, with irregularly spaced time points, TLCC might introduce bias when estimating time lag and coefficients. Given the constraint of TLCC, which can only generate time lag estimates from the observed time points, its performance was compromised when applied to this dataset with irregular time gaps (Figure 5). However, GPlag, by estimating a continuous time lag, was able to identify lags that restored strong similarity between the pairwise time series.

Thirdly, to evaluate if our candidate enhancer-promoter pairs were detected as having a functional relationship in an orthogonal experimental data type, we downloaded the genomic locations of response expression quantitative trait loci (response eQTL) found in human macrophages exposed to IFN_γ (Alasoo et al., 2018) (that is, using one of the two stimuli used in the Reed et al. (2022) experiment, and in the same cell type). Because the number of eQTLs was small compared to the number of enhancers, we used a relatively stringent threshold to label a putative functional enhancer-promoter pair: we split the candidate enhancer-promoter pairs into groups according to posterior mean correlation ($>0.98$). There were 74 H3K27ac peaks overlapping with macrophage response eQTL, among which 12 pairs were in our defined candidate enhancer-promoter pair sets. Therefore, there was a significant enrichment of our candidate enhancer-promoter pairs to eQTLs ($\chi^{2}$ test $p=0.18\times 10^{-3}$) futher suggesting the relevance of the enhancer-promoter pairs highly ranked by GPlag.

GPlag is a flexible and interpretable model specifically designed for irregular time series with limited time points that exhibit lead-lag effects, a scenario often encountered in biology, which inspired its development. These irregularities can be caused by missing values or study design restrictions. We believe that our method’s adaptability renders it valuable for time series applications in various areas such as longitudinal gene expression, epigenetics, proteomics, cytometry, and microbiome studies. It can be applied to a variety of problems including but not limited to ranking or clustering time series based on the degree of lead-lag effects, and interpolation or forecasting. The core feature of our methodology is the transform of the multi-output GP into a single-output GP, which facilitates the identification of parameters. This approach is grounded in theoretical support, affirming its validity as a covariance kernel for GP. This transform leverages the shared information among time series and enables simultaneous modeling of the time lag parameter $s$ and dissimilarity parameters $a$, thereby achieving desirable interpretability and expanding the application of the model. Additionally, our model can be extended to accommodate multiple time series, beyond just pairwise time series, enhancing its applicability and versatility.

The work here suggests some exciting future directions. First, when the sample size $n$ is small, the choice of the kernel is crucial in GP models. The extension to a non-stationary kernel allows for more flexibility when modelling non-stationary time series with more heterogeneity. The two main challenges are to find a valid non-stationary kernel with interpretable parameters and perform efficient and effective statistical inference. On the other hand, for time series data with large $n$, i.e., dense time points, the vanilla GP is computationally expensive with $\mathrm{O}(n^{3})$ complexity (Rasmussen, 2004). In this situation, we can apply scalable GP methods such as nearest neighbor Gaussian processes (NNGP, Datta et al. (2016)) or variational inference GP (Tran et al., 2015) to reduce the complexity to $\mathrm{O}(n\log n)$ for large datasets. In addition, incorporating additional parameters, such as varying the variance $\sigma^{2}$ and lengthscale $b$ across different time series, could potentially enhance model performance. Finally, the question of whether the parameters in Theorem 3 can be consistently estimated remains open. Even for simpler kernels such as the Matérn, consistent estimators for identifiable parameters (e.g., $\sigma^{2}b^{2\nu}$ and $\nu$) were only recently developed (Loh et al., 2021; Loh and Sun, 2023). For more complex kernels, like the ones proposed in this paper, constructing consistent estimators is an open and challenging problem.