Inferring directed spectral information flow between mixed-frequency time series

TL;DR

This paper introduces a new method to analyze how information flows between time series with mixed frequencies, outperforming traditional approaches in accuracy and efficiency.

Contribution

The paper proposes a nonparametric method (MF-TFCCA) for assessing directed spectral information flow in mixed-frequency time series with nonlinear interactions.

Findings

MF-TFCCA outperforms traditional parametric models in computational efficiency and detection accuracy.

The method successfully recovers dominant driving frequencies in mixed-frequency time series.

MF-TFCCA is validated on real-world data from finance, climate, and neuroscience.

Abstract

Identifying directed spectral information flow between multivariate time series is important for many applications in finance, climate, geophysics and neuroscience. Spectral Granger causality (SGC) is a prediction-based measure characterizing directed information flow at specific oscillatory frequencies. However, traditional vector autoregressive (VAR) approaches are insufficient to assess SGC when time series have mixed frequencies (MF) or are coupled by nonlinearity. Here we propose a time-frequency canonical correlation analysis approach (“MF-TFCCA”) to assess the strength and driving frequency of spectral information flow. We validate the approach with extensive computer simulations on MF time series under various interaction conditions and further assess statistical significance of the estimate with surrogate data. In various benchmark comparisons, MF-TFCCA consistently outperforms the traditional parametric MF-VAR model in both computational efficiency and detection accuracy, and recovers the dominant driving frequencies. We further apply MF-TFCCA to real-life finance, climate and neuroscience data. Our analysis framework provides an exploratory and computationally efficient nonparametric approach to quantify directed information flow between MF time series in the presence of complex and nonlinear interactions.

Figures

b](#F1)): 1) inferring directionality of information flow; 2) identifying the driving frequency; and 3) quantifying the relative information flow metric (normalized by self-predictability in analogy to autocorrelation) and assessing it statistical significance (related to surrogate data). For the first subtask, we first use MF-TFCCA to compute the CC between a time series and a time-frequency process, where the time-frequency representation (TFR) can be computed by a moving window short-time Fourier transform (STFT). We further determine the directionality based on the asymmetric shape of lag-

a](#F2). We then produced a modulated signal $[eqn]$ that was coupled with the phase of $[eqn]$ , yielding a nonlinear coupling between $[eqn]$ and $[eqn]$ (Methods and [Fig. 3a](#F3)). Upon down-sampling of $[eqn]$ , we found that xcorr ([Fig. 3b](#F3)) or MF-TFCCA correlating complex-valued STFT of $[eqn]$ with $[eqn]$ ([Fig. 3c](#F3)) failed to discover the nonlinear causal relationship. In contrast, MF-TFCCA correlating the STFT amplitude of $[eqn]$ with $[eqn]$ could uncover the nonlinear GC relationship ([Fig. 3d](#F3)). The reason for such result discrepancy is that the causal informati

a](#F4)). More specifically, $[eqn]$ Granger causes $[eqn]$ at two driving frequencies (5 Hz and 80 Hz), and $[eqn]$ drives $[eqn]$ at one frequency (50 Hz), with the MF setting as $[eqn]$ and $[eqn]$ ([Fig. 4b](#F4)). Direct xcorr analysis suggested an asymmetric information flow between $[eqn]$ and $[eqn]$ ([Fig. 4c](#F4)). Provided that we only had access to $[eqn]$ and $[eqn]$ from a partially observed system, the inferred SGC would suggest a direct causality between $[eqn]$ and $[eqn]$ because of the latent variable $[eqn]$ ; however, the SGC inferred from the fully observed system would

c](#F5)). MF-TFCCA was more computationally efficient than MF-VAR (by about 1–2 order of magnitude) at various settings ([Fig. 5d](#F5)). From the lag-CC curve ([Fig. 5e](#F5)), we inferred a bidirectional information flow. The CCA coefficients (absolute values) at both directions revealed two peaks, suggesting two driving frequencies (one at theta and another at gamma band). To further recover the dominant driving frequency, we systematically conducted the “filter-one-frequency-out” analysis and found a greater CC contribution at the theta-band driving frequency from its relative GC gains ([F

a](#F6)). The CPI and OIL shares the same sampling rate, and standard VAR method detects GC from OIL to CPI ([Fig. 6b](#F6)). We then used TF-CCA to estimate the relationship between OIL and GDP and the relationship between CPI and GDP. The two lag-CC profiles revealed the causality from OIL to GDP ([Fig. 6c](#F6)), and a weaker causality from CPI to GDP ([Fig. 6d](#F6)). Next, we used partial CCA to assess the effect of CPI on the relationship between OIL and GDP, and found that causality relationship between OIL and GDP diminished ([Fig. 6e](#F6)). In conclusion, MF-TFCCA and GC methods reve

a](#F7)). The goal here is to infer the causal relationship between monthly precipitation and monthly temperature, monthly precipitation and yearly temperature, and monthly temperature and yearly precipitation. Interestingly, there was only weak or non-significant GC between rainfall and temperature after removing the linear trend in the time series ([Fig. 7b](#F7)). By examining the MF time series, however, MF-TFCCA discovered a directed information flow between monthly rainfall and yearly temperature time series ([Fig. 7c](#F7)), but not between monthly temperature and yearly rainfall time s

b](#F8)). Based on the previously detected hippocampal ripple events during non-rapid-eye-movement (NREM) sleep (n=753), we conducted STFT on the filtered hippocampal ripple band (100–250 Hz) activity using a 40-ms moving window with 30-ms overlapping window. Comparing the lag-CC profile and peak timing before, during and after the detected hippocampal ripple events ([Fig. 8c](#F8), [i](#F8)–[iii](#F8)), we detected a bidirectional information flow between the RSC and hippocampus (HPC), with a slightly stronger strength in the RSC→HPC direction before and during hippocampal ripples. Additional

Tables

Step 1: We convert univariate X(t) in the time domain into a multivariate TFRX~(t,ω) in the time-frequency domain (such as STFT) to match the sampling rate of Y(t). The window size Wl of STFT is chosen based on the sampling rate of X and desired frequency resolution: the longer the window, the higher frequency resolution. The overlap of moving window is (Wl-1/Fy). For calibration, we apply z-score normalization on each frequency of the complex-valued time-frequency matrix X~(t,ω).

Note: the z-score normalization on each frequency implies prewhitening of X~ in the frequency-axis but not in the time-axis. The degree of autocorrelation in the time-axis depends on the size of overlapping window: a larger overlap implies a higher degree of autocorrelation.

Step 2: We use regularized CCA to assess CC where the linear projection of X~ is maximally correlated with Y, from which we compute the frequency vs. CCA coefficient (absolute value) and identify a canonical frequency (say f0) that is associated with the peak CCA coefficient. The canonical frequency often corresponds to the dominant oscillatory frequency fx in the source signal X. The value of CC depends on both temporal and frequency resolution of X~ (i.e., the autocorrelation of the source signal) as well as the temporal resolution of Y. To calibrate the CC value, we may normalize CC with an autocorrelation-like CC factor, which is derived from MF-TFCCA between X~ and a down-sampled version of X(t).

Note 1: CC will be computed between real-valued Y and complex-valued X~=STFT(X). We may consider two options for assessing CC and CCA coefficients. Option 1: first applying complex-valued CCA and then taking the magnitude of complex-valued CCA coefficients (CC is always real-valued). Option 2: first taking the absolute value abs(X~)=abs(STFT(X)) and then applying standard CCA to obtain real-valued CCA coefficients and CC.

Note 2: As a sanity check, if we filter out the dominant oscillatory frequency fx in the source signal X (i.e., “filter-one-frequency-out”) and repeat the same procedure, the resultant CC shall be significantly reduced (i.e., a negative CC gain). Therefore, the spectral contributions of individual frequencies can be assessed by comparing their relative CC gains (compared to the original unfiltered setting).

Step 3: We shift the time series by a lag parameter and repeat the above procedure to construct a lag-CC curve. The lag-CC curve is usually asymmetric with respect to 0. A large CC value at negative lag indicates that X→Y in information flow, whereas a large CC value at positive lag indicates an opposite direction Y→X. If the lag-CC curve is approximately symmetric, an additional step (such as “filter-one-frequency-out”) is required to identify whether causality is bidirectional.

Step 4: To assess the statistical significance of lag-CC, we create the surrogate data of both X and Y, and generate a null distribution of CC. If the original CC value is significantly greater than the surrogate data statistics, a significance outcome will be reported (t-test’s P-value for multiple trials, or two-sample Kolmogorov-Smirnov test’s P-value). By comparing the lag-CC profiles derived from the original and surrogate data, we may determine the directed information flow between X and Y.

Full text

INTRODUCTION

Quantifying directed information flow between two or more simultaneously measured time series is a common problem in science and engineering, including economics, finance, meteorology, transportation, geophysics, and neuroscience^1–4^ (Granger, 1969; Surgihara et al., 2012; Silva et al., 2021; Yuan and Shou, 2022). These time series are sometimes collected at different temporal resolution. For instance, this may occur in finance by comparing a daily stock price with a weekly trade transaction; or may occur in weather forecast by linking monthly rainfall amount at the city level to yearly rainfall at the national level; or may occur in neuroscience when concurrent neural recordings of different modalities (e.g., calcium imaging and electrophysiology; concurrent EEG-fMRI) are simultaneously collected (Fig. 1a).

Statistical dependency of two time series may be characterized in the time domain (such as lagged cross-correlation, or “xcorr”) or in the frequency domain (such as coherence), which bears many other names in different contexts, such as functional connectivity^5^ (Bastos and Schoffelen, 2016), partial directed coherence (PDC)^6^ (Baccala and Sameshima, 2001), Granger causality (GC) or GC-like measures^7–11^ (Granger, 1980; Geweke, 1984; Ding et al., 2006; Barnett and Seth, 2014; Shojaie and Fox, 2022), and transfer entropy^12^ (Schreiber, 2000). These measures are either undirected or directed, where the direction implies an asymmetry in information flow between the “source” and “sin” (or “target”). Two generalized scenarios are noteworthy. First, multiple (greater than two) time series are observed; in this case, multivariate statistical methods have been developed to assess their conditional dependence (e.g., conditional GC, partial correlation, partial directed coherence). Second, two time series have distinct temporal resolution; in this case, different analysis strategies are required^11^ (Shojaie and Fox, 2022). One naive strategy is to apply subsampling or temporal averaging to the signal with higher sampling frequency to match the one with lower sampling frequency, and then to characterize their interactions based on standard metrics. Another strategy originated from the econometrics literature is GC testing with mixed frequency (MF) data using a MF-VAR model^13–15^ (Ghysels et al., 2016; Gotz and Hecq, 2019; Anderson et al., 2016). The MF-VAR method, however, only works well when the discrepancy of MF is small or the lag is relatively small, and cannot account for nonlinear interactions. Several lines of work have recently studied the identifiability of a structural VAR(1) model under arbitrary subsampling^16^ (Gong et al., 2015) and MF^17^ (Tank et al., 2019) settings. The nonlinear interactions have also been studied in the framework of non-Gaussian structural VAR model (e.g., REF^18^ Hyvarinen et al., 2010). To our best knowledge, however, these approaches rely on parametric model estimation (such as REF^16,17^ Gong et al., 2015; Tank et al., 2019), and cannot recover the driving frequency in the SGC method.

In this paper, we propose a computationally efficient, nonparametric approach to address the above-mentioned challenges, identify the prominent driving frequency (in SGC sense), and assess directed spectral information flow between MF time series. Our approach is motivated from canonical correlation analysis (CCA) in the time-frequency (TF) domain but adapted for MF signals, which we coin “MF-TFCCA” to distinguish from the traditional canonical GC^19,20^ (Sato et al., 2010; Ashrafulla et al., 2013). To validate our approach, we conduct systematic computer simulations on bivariate or trivariate VAR processes (with known SGC ground truth) to validate our analysis pipeline, and further apply MF-TFCCA to three real-life datasets. In the case of neuroscience data, our analysis results provide insights into the directed neural interactions between the mouse hippocampus and downstream cortical structures, and further reveal directed hippocampal-neocortical information flow before, during and after hippocampal sharp wave ripples (SWRs) in sleep.

A methodological overview of GC and directed information flow

Correlation is the most statistical measure that characterizes linear dependence of two random variables. In the frequency domain, the complex-valued coherency corresponds to the Fourier transform of cross-correlation, and the modulus of the coherency is coherence. While correlation and causation can coexist, correlation does not imply causation. Additionally, two variables can be statistically uncorrelated even if one is caused by the other based on nonlinear transformation (e.g., $[eqn]$ ). Therefore, the relationship between causality and correlation depends on multiple factors, such as linearity and Gaussianity.

Identifying causality in a complex physical or physiological system is an important task. In such systems, the recorded signals may have single or multiple dominant oscillatory components in the frequency domain. The notion of GC is different from causality in the physical sense. GC or GC-like methods provide a framework that uses temporal predictability as opposed to correlation to identify causation between time-series variables. Specifically, a variable $[eqn]$ is said to “Granger cause” a variable $[eqn]$ (notation: as $[eqn]$ ) if the predictability of the target variable $[eqn]$ declines when the source variable $[eqn]$ is removed from the universe of all possible causative variables^9^ (Ding et al., 2006). GC can include both lagged and instantaneous effects between time series. This assumes that causes can be separated from effects, so that a variable is identified as causative if predictability declines when that variable is removed. Here, we restrict the GC effect on the history information only and exclude the instantaneous causal effect. Although GC in the original definition is in a general sense^7^ (Granger, 1980), its practical use is often based on linear predictability assumption and Gaussian assumption, and has been facilitated by well-documented software^10^ (Barnett and Seth, 2014). Within the VAR framework, directed GC $[eqn]$ in the time domain can be also expressed equivalently in terms of SGC in the frequency domain, which displays as a high amplitude at a single driving frequency (or the most prominent driving frequency among several modes). The driving frequency is usually linked to the inherent oscillatory frequencies of the “source” signal. In a general setting of bidirectional GC (denoted as $[eqn]$ ), two driving frequencies in $[eqn]$ and $[eqn]$ directions are different, and $[eqn]$ and $[eqn]$ can play both the roles of source and target but at different timescales. GC is a special case of transfer entropy under the linear Gaussian assumption, which defines a direct information transfer between jointly dependent random processes^12,21,22^ (Schreiber, 2000; Barnett et al., 2009; Stramaglia et al., 2021).

Since GC is a prediction-based measure, it is generally impossible to determine GC from the asymmetric shape of “xcorr” between two random variables in the time domain^23^ (Borysov and Balatsky, 2014). Lagged cross-correlation is neither sufficient nor necessary condition for GC (see Supplementary Fig. 1 for simple examples). However, the lead-lag relationship may provide a hint to the directionality of information flow (but not the driving frequency). Notably, cross-correlation function is influenced not only by the relationship between the signals but also by the autocorrelation of the signals. In other words, autocorrelation within the target variable can also help prediction, leading to a spurious GC from the regressor to the target variable. Therefore, signal prewhitening may override the limitation of revealing the causal relationship between random variables^24^ (El-Gohary and McNames, 2007). CCA has also been used to determine linear Granger causal relationship^19,20,25^ (Sato et al., 2010; Ashrafulla et al., 2013; Zhuang et al., 2020). A naïve canonical GC approach is to maximize the canonical correlation (CC) between eigen-time series at different time lags while removing the influence of autocorrelation within the target time series^19^ (Sato et al., 2010). Notably, there is a conceptual difference between CCA-based GC analysis and VAR-based GC analysis in that the former approach is based on two sets of variables, and the latter approach allows a linear combination of more than one variable in the regressor.

Finally, it is uncommon to observe that time series have different temporal resolution induced by subsampling or due to data acquisition discrepancy. Traditional GC approaches do not work for signals with dramatically different temporal resolution. It is also well known that subsampling may produce a bias for the GC estimate^26–28^ (Solo, 2007; Barnett and Seith, 2017; Anderson et al., 2019). Several methods have been developed to account for nonlinear and non-Gaussian observations and allow for subsampled and MF time series^11^ (Shojaie and Fox, 2022), but remains unclear how various MF and nonlinear coupling conditions may affect the bidirectional SGC estimate.

Conceptual ideas of our approach

To gain intuition, let’s start our reasoning with a simple example. Let’s consider $[eqn]$ and $[eqn]$ as two univariate random variables from respective time series $[eqn]$ and $[eqn]$ , with the known ground truth: $[eqn]$ Granger causes $[eqn]$ at a driving oscillatory frequency $[eqn]$ , denoted as $[eqn]$ . From a linear system perspective, a directed information flow $[eqn]$ suggests that the history of $[eqn]$ improves linear predictability of $[eqn]$ , or $[eqn]$ contains linearly filtered information from $[eqn]$ . However, if the interaction component acts like a high-pass filter, the causal information transfer $[eqn]$ will be largely lost after down-sampling of $[eqn]$ because the autocovariance structure of $[eqn]$ does not contain all information required to estimate the original VAR model. Therefore, the temporal causal relationship in the high-frequency range may not be recovered.

When there is a unidirectional GC relationship $[eqn]$ and two variables $[eqn]$ and $[eqn]$ have MF, we consider two possibilities: the sampling rate of $[eqn]$ , denoted by $[eqn]$ , is greater than the sampling rate of $[eqn]$ , denoted by $[eqn]$ ; or the opposite, namely $[eqn]$ has a higher sampling rate than $[eqn]$ such that $[eqn]$ . In both cases, if $[eqn]$ , then it is impossible to recover the causal relationship in light of Nyquist’s sampling theorem. Only if $[eqn]$ , inferring the driving frequency $[eqn]$ becomes possible. When there is a bidirectional GC relationship $[eqn]$ and two directions have distinct driving frequencies, the outcome may become more complex and depend on multiple factors, which will be illustrated in the subsequent Results section.

The problem of identifying directed spectral information flow consists of three subtasks (Fig. 1b): 1) inferring directionality of information flow; 2) identifying the driving frequency; and 3) quantifying the relative information flow metric (normalized by self-predictability in analogy to autocorrelation) and assessing it statistical significance (related to surrogate data). For the first subtask, we first use MF-TFCCA to compute the CC between a time series and a time-frequency process, where the time-frequency representation (TFR) can be computed by a moving window short-time Fourier transform (STFT). We further determine the directionality based on the asymmetric shape of lag-CC derived from MF-TFCCA by moving the time lag between two time series. If the peak of lag-CC curve deviates from zero (at either positive or negative lag), it implies that the history of one variable can improve the predictability of the other variable. For the second subtask, the driven frequency (or frequencies) may be recovered by the dominant peak (or peaks) of CCA coefficients at the frequency axis. Furthermore, the contribution of individual driving frequencies can be assessed by the relative CC gain based on the “filter-one-frequency-out” strategy. For the third subtask, we resort on permutation statistics using surrogate data^4,10,29^ (Yuan and Shou, 2022; Barnett and Seth, 2014; Lancaster et al., 2018). The purpose of the surrogate data is to destroy the signal’s temporal regularity by introducing randomness in phase while maintaining the magnitude in the Fourier domain.

RESULTS

Analysis pipeline

We first consider the condition where $[eqn]$ and $[eqn]$ are two univariate time series, and then relax the assumption to more than two time series. Without the loss of generality, unless specified otherwise, we assume that $[eqn]$ has a higher sampling rate than $[eqn]$ , and call $[eqn]$ as a high-frequency (HF) process, and $[eqn]$ as a low-frequency (LF) process. The MF-TFCCA analysis pipeline (Fig. 1c) is described in Box 1.

**: **

In the following sections, we will validate MF-TFCCA using computer simulations and real-world datasets. A summary of computer simulation results is given in Supplementary Table 1.

Computer Simulations

Bivariate bidirectional systems with distinct driving frequencies.

We started with a bivariate and bidirectional GC system $[eqn]$ and $[eqn]$ . The system has two variables that mutually excite each other at different time lags at two different frequencies: $[eqn]$ and $[eqn]$ (Fig. 2a). The original time series were generated from a bivariate 4th-order AR process with a sampling rate of 200 Hz (Methods), with a total of 100 trials and 20 s per trial (Fig. 2b). After down-sampling of $[eqn]$ by a factor of 5, two time series had distinct temporal resolution: $[eqn]$ and $[eqn]$ . At the first examination of lag cross-correlation (“xcorr”), there was two noticeable positive peaks: one peak at negative lag of ~150 ms and the other peak at positive lag of ~100 ms (Fig. 2c). To compute the STFT of $[eqn]$ , the window length was chosen to be 200 ms, and the overlap for the moving window was chosen as 200–25=175 ms. Consequently, the time-frequency matrix $[eqn]$ had the same length in the time-axis as the time series vector $[eqn]$ . In this case, MF-TFCCA recovered similar lag-CC shapes as in xcorr: the peak at the negative lag indicates $[eqn]$ and the peak at the positive lag indicates $[eqn]$ (Fig. 2d). Furthermore, examining the absolute CCA coefficients showed a peak at 4 Hz in the $[eqn]$ direction, and a high peak at 15 Hz in the $[eqn]$ direction. We also witnessed an energy leakage at 4 Hz because of the close distance between two driving frequencies.

Although in this example, MF-TFCCA generated similar GC results consistent with xcorr and recovered the driving oscillatory frequencies. This observation could change depending on the relationship of driving frequency and sampling frequency. To illustrate this point, we also considered a special case when GC is unidirectional (either $[eqn]$ with $[eqn]$ and 4 Hz or $[eqn]$ with $[eqn]$ ; Methods). We assumed that $[eqn]$ had intrinsic oscillations at 4 Hz and 80 Hz, whereas $[eqn]$ had intrinsic oscillations at 2 Hz and 15 Hz (Supplementary Fig. 2a). First, we tested the direction of $[eqn]$ and down-sampled $[eqn]$ by a factor of 5 (Supplementary Fig. 2b). We computed the STFT of $[eqn]$ using a window size of 150 ms (Supplementary Fig. 2c). In this case, the xcorr profile showed an asymmetric shape between the original $[eqn]$ and subsampled $[eqn]$ (Supplementary Fig. 2d). We also recovered an asymmetric lag-CC profile, suggesting that the direction of information flow is from $[eqn]$ to $[eqn]$ (Supplementary Fig. 2e). The CCA coefficients showed two peaks at 4 Hz and 80 Hz (Supplementary Fig. 2f), which matched two dominant driving frequencies. Note that although the SGC $[eqn]$ was higher in 80 Hz than in 4 Hz, the causal relationship in the high-frequency band was suppressed after down-sampling of $[eqn]$ , resulting in a lower peak of CCA coefficients in 80 Hz than in 4 Hz. A follow-up “filter-one-frequency-out” analysis revealed the relative CC contributions of these two driving frequencies: a large negative CC gain suggests an important contribution of the putative driving frequency (i.e., 4 Hz) within the filtered-out frequency band, whereas a nearly zero CC gain suggests negligible contribution of the other putative driving frequency (Supplementary Fig. 2g). The reason why we didn’t observe a clear negative CC gain at 80 Hz is probably the spectral causality was main preserved in the low frequency range after down-sampling of $[eqn]$ .

Further, we tested the direction of $[eqn]$ and down-sampled $[eqn]$ by a factor of 5 (Supplementary Fig. 2h). In this case, xcorr showed a nearly symmetric shape (Supplementary Fig. 2i) and MF-TFCCA uncovered the true GC relationship and information flow from $[eqn]$ to $[eqn]$ . Importantly, comparing the $[eqn]$ and $[eqn]$ directions, the estimated lag-CC values (0.8 vs. 0.4; Supplementary Fig. 2j vs. Supplementary Fig. 2e) also quantitatively preserved the relationship of relative SGC values (5 vs. 2.5; Supplementary Fig. 2h vs. Supplementary Fig. 2a). The CCA coefficients recovered one dominant peak at the 15 Hz driving frequency and second peak at 80 Hz (false positive) (Supplementary Fig. 2k). Similarly, a subsequent “filter-one-frequency-out” analysis revealed the relative CC contribution between these two putative driving frequencies (Supplementary Fig. 2l).

Next, we investigated the impact of down-sampling factor of $[eqn]$ on the lag-CC estimate. Using the above unidirectional system as an example ( $[eqn]$ and 4 Hz), we systematically varied the sampling frequency of $[eqn]$ while keeping the sampling rate of $[eqn]$ intact and used MF-TFCCA based on different window sizes. To accommodate the change in ratio $[eqn]$ , we adapted the overlap of the moving window to assure the same data length between $[eqn]$ and $[eqn]$ . We found that the overall results were quite robust with respect to sampling frequency ratio and window size: the estimated lag-CC profiles remained consistent, and the shape looked qualitatively similar (Supplementary Fig. 3). The window size affected the shape of the lag-CC profile, as an increasing window size decreased the temporal resolution of the lag-CC profile. Meanwhile, all CC statistics computed from surrogate data were significantly smaller.

Bivariate GC systems with nonlinear coupling.

We further asked the question: how does the nonlinearity affect our estimation outcome in terms of identifiability? Specifically, we considered two variables that are Granger causally modulated by specific nonlinearity. We constructed two different types of nonlinear GC: one through phase-amplitude coupling (PAC), and the other through nonlinear amplitude coupling.

For PAC nonlinearity, we first generated raw signals $[eqn]$ and $[eqn]$ using a linear bidirectional causal system as Fig. 2a. We then produced a modulated signal $[eqn]$ that was coupled with the phase of $[eqn]$ , yielding a nonlinear coupling between $[eqn]$ and $[eqn]$ (Methods and Fig. 3a). Upon down-sampling of $[eqn]$ , we found that xcorr (Fig. 3b) or MF-TFCCA correlating complex-valued STFT of $[eqn]$ with $[eqn]$ (Fig. 3c) failed to discover the nonlinear causal relationship. In contrast, MF-TFCCA correlating the STFT amplitude of $[eqn]$ with $[eqn]$ could uncover the nonlinear GC relationship (Fig. 3d). The reason for such result discrepancy is that the causal information was contained in the spectral power of $[eqn]$ in the PAC system but not in its fluctuations. The periodic nature of $[eqn]$ and the random phases in the complex-valued STFT lead to a reduced lag-CC profile. In contrast, the real-valued STFT considers only the amplitude component of the signal spectrum and is capable of recovering the true causality in the PAC system. Additionally, the CCA coefficients (in absolute value) at both positive and negative lags showed peaks in the HF range where $[eqn]$ stored the causal information of $[eqn]$ (Fig. 3e).

For nonlinear amplitude coupling, we first generated raw signals $[eqn]$ and $[eqn]$ using a linear bidirectional causal VAR system (Supplementary Fig. 4a) and then produced an amplitude-modulated signal $[eqn]$ by applying a nonlinear memoryless function of $[eqn]$ (such as a sigmoid function or a periodic cosine function; Supplementary Fig. 4b,c; Methods). To detect such nonlinear causal relationship, we similarly applied MF-TFCCA correlating the complex-valued STFT ( $[eqn]$ ) with $[eqn]$ . Consequently, the lag-CC profile successfully recovered nonlinear GC relationship, which was significantly greater than the CC statistics of the surrogate data. As the amplitude modulation only enhances or suppresses the signals without altering the causal structure, MF-TFCCA with complex-valued STFT successfully revealed the causal relationship. On the opposite direction, when we applied nonlinear amplitude modulation to $[eqn]$ , and the results were similar (Supplementary Fig. 4d,e). Together, these results suggest that MF-TFCCA was capable of detecting nonlinear GC coupling between MF time series.

Trivariate systems with chain and parallel GC systems.

To generalize the bivariate system to a trivariate system, we first considered the chain system and assumed that $[eqn]$ with the known ground truth SGC profiles: $[eqn]$ and $[eqn]$ (Fig. 4a). More specifically, $[eqn]$ Granger causes $[eqn]$ at two driving frequencies (5 Hz and 80 Hz), and $[eqn]$ drives $[eqn]$ at one frequency (50 Hz), with the MF setting as $[eqn]$ and $[eqn]$ (Fig. 4b). Direct xcorr analysis suggested an asymmetric information flow between $[eqn]$ and $[eqn]$ (Fig. 4c). Provided that we only had access to $[eqn]$ and $[eqn]$ from a partially observed system, the inferred SGC would suggest a direct causality between $[eqn]$ and $[eqn]$ because of the latent variable $[eqn]$ ; however, the SGC inferred from the fully observed system would suggest negligible $[eqn]$ information flow. In comparison, MF-TFCCA was also able to distinguish these two cases, where the partial CC (i.e., $[eqn]$ ) was statistically insignificant (Fig. 4d,e).

Next, we considered the other forms of trivariate chain system $[eqn]$ , with the MF setting as $[eqn]$ and $[eqn]$ (Supplementary Fig. 5a and Methods). We computed the STFT for both $[eqn]$ and $[eqn]$ (Supplementary Fig. 5b,c). In this example, without considering $[eqn]$ , the lag-CC profile identified the directed information flow $[eqn]$ (Supplementary Fig. 5d); even partializing out the intermediate LF variable $[eqn]$ , the lag-CC profile still produced a similar result, yielding a false positive (Supplementary Fig. 5e). Alternatively, we could consider partializing out the STFT spectrum of $[eqn]$ and repeated the analysis procedure. As a result, the lag-CC magnitude reduced but the estimate remained to be a false positive due to the information loss of down-sampling of $[eqn]$ .

Furthermore, we considered a trivariate parallel system $[eqn]$ and $[eqn]$ with a similar setting (Supplementary Fig. 6a–c), with the MF setting as $[eqn]$ and $[eqn]$ . The xcorr profile suggested a directed information $[eqn]$ (Supplementary Fig. 6c). Additionally, the partial CC had the same estimate of standard CC (Supplementary Fig. 6d,e).

Two-species Logistic model and coupled Rössler-Lorenz system.

Next, we validated our MF-TFCCA method with a two-dimensional deterministic nonlinear system: two-species Logistic model^30^ (Ma et al., 2014), where the two generated time series can be either unidirectional ( $[eqn]$ , Supplementary Fig. 7a) or bidirectional ( $[eqn]$ , Supplementary Fig. 7e) depending on the model parameters (Methods). As comparison, we also computed the standard GC between $[eqn]$ and $[eqn]$ using the evenly-sampled time series (Supplementary Fig. 7b,c,f,g). We further down-sampled $[eqn]$ and selected the sampling frequency ratio $[eqn]$ as 5:1. In both uni- and bidirectional cases, MF-TFCCA succeeded in recovering the correct directionality of information flow in the presence of nonlinearity (Supplementary Fig. 7d,h).

We further considered a nonlinear coupled system^30^ (Ma et al., 2014), where 3-dimensional Rössler system $[eqn]$ drives the 3-dimensional Lorenz system $[eqn]$ , where both systems are deterministic chaotic and coupled through a unidirectional coupling term from $[eqn]$ to $[eqn]$ (Methods). The ground truth causality of the system is shown in Supplementary Fig. 8a. The goal was to infer the causality relationship between $[eqn]$ and down-sampled $[eqn]$ by a factor of 5 (Supplementary Fig. 8b). As a comparison, we first computed the pairwise SGC between using the evenly-sampled time series (Supplementary Fig. 8c,d), and found significant information flow in many pairwise directions, including $[eqn]$ and $[eqn]$ directions. Noticeably, the SGC strength was greater in $[eqn]$ than $[eqn]$ direction. Additionally, some false-positive casual relationships were reported. In the MF time series setting, we applied MF-TFCCA to infer the relationship between $[eqn]$ and $[eqn]$ using real-valued version of $[eqn]$ (Supplementary Fig. 8e). In general, we found strong information flow in $[eqn]$ and $[eqn]$ directions. Notably, false-positive causal relationships were also detected from $[eqn]$ to $[eqn]$ . This result might be due to (i) the high periodicity of the chaotic systems, which could mislead the linear method to infer a false causality direction; (ii) the high correlation between the variables of the chaotic systems, such as between $[eqn]$ and $[eqn]$ and between $[eqn]$ and $[eqn]$ . Together, these results suggest that the causality coupled through a highly nonlinear system may not be captured by the linear method.

Benchmark comparison

We further made systematic benchmark comparisons on GC inference (qualitatively) and computational efficiency (quantitatively) between our proposed MF-TFCCA and the state-of-the-art MF-VAR method. These results are summarized in Supplementary Table 2 and Supplementary Fig. 9.

First, more frequently, MF-VAR produced false-positive results in several systems. For example, MF-VAR mistakenly detected bidirectional causality in a truly unidirectional $[eqn]$ system, and mistakenly detected $[eqn]$ in a fully observed $[eqn]$ chain system. It even mistakenly detected an unreal bidirectional causality between $[eqn]$ and $[eqn]$ in a parallel system. In contrast, MF-TFCCA detected less false-positives. Moreover, MF-VAR failed to detect GC relationship in the presence of nonlinear coupling.

Second, MF-VAR can only detect GC in the time domain, whereas MF-TFCCA is capable of uncovering the driving frequency (i.e., SGC). When there are multiple driving frequencies present in the system, we can also assess the relative contribution of each driving frequency using the “filter-one-frequency-out” strategy.

More importantly, MF-TFCCA showed a greater advantage in computational efficiency compared to MF-VAR. In dealing with multiple-trial data structure, MF-TFCCA maintained an approximately linear time complexity, while the computational time of the MF-VAR inference method increased substantially (slower than quadratically) even in the simplest $[eqn]$ setting. In the general $[eqn]$ setting with a larger order $[eqn]$ , the computational bottleneck is much worse for MF-VAR; whereas MF-TFCCA is nonparametric and independent of the VAR model order. Moreover, MF-TFCCA is efficient when computing across a broad range of time lags, whereas MF-VAR can only handle one lag each time. For instance, it took only a few minutes for MF-TFCCA to analyze a 100-trial trivariate system across a broad range of lags. In contrast, MF-VAR required around an hour to analyze a simple 2-trial trivariate system with two HF variables and one LF variable at only one lag. Therefore, it is nearly infeasible for MF-VAR to detect GC relationship with an uncertain lag or detect bidirectional causality with different lags (the null result on the bidirectional system in Fig. 2 serves as an example).

Finally, we conducted a benchmark comparison on a real-world neuroscience dataset using concurrent local field potential (LFP) recordings of the rat primary somatosensory cortex (S1) and anterior cingulate cortex (ACC) during thermal pain experiments (Fig. 5a; Methods). Consistent with our prior investigation^31^ (Guo et al., 2020), we identified bidirectional SGC between the S1 and ACC in chronic pain-treated rats, where GC was prominent at theta (4–8 Hz) and gamma (50–70 Hz) bands at both $[eqn]$ and $[eqn]$ directions derived from a bivariate VAR(15) model (Fig. 5b). Treating that as a “ground truth”, we down-sampled the original sampling rate of LFP signals from one brain region by 4 folds (from 200 Hz to 50 Hz) and compared MF-TFCCA and MF-VAR results with the ground truth. Since our method is nonparametric, we also compared the parametric models with different orders: MF-VAR(1) and MF-VAR(5).

We first down-sampled S1 LFP signals and kept ACC LFP signals intact. We found that MF-TFCCA was able to identify bidirectional GC, whereas MF-VAR only identified unidirectional GC (Fig. 5c). MF-TFCCA was more computationally efficient than MF-VAR (by about 1–2 order of magnitude) at various settings (Fig. 5d). From the lag-CC curve (Fig. 5e), we inferred a bidirectional information flow. The CCA coefficients (absolute values) at both directions revealed two peaks, suggesting two driving frequencies (one at theta and another at gamma band). To further recover the dominant driving frequency, we systematically conducted the “filter-one-frequency-out” analysis and found a greater CC contribution at the theta-band driving frequency from its relative GC gains (Fig. 5f). Notably, both down-sampling (up to 50 Hz) and short-duration trials (1 second) and produced limited frequency resolution, preventing us to recover the gamma-band driving frequency. In contrast, the MF-VAR models were incapable of conducting similar analyses. At last, we reversed the direction by down-sampling ACC LFP signals and repeated the comparison, and the conclusion remained similar (Fig. 5g–j).

Together, these extensive benchmark comparisons pointed to the advantages of our proposed method compared to the state-of-the-art MF-VAR method. MF-TFCCA not only can correctly identify bidirectional spectral information flow, but also can uncover the driving frequency or quantify the relative contribution among multiple driving frequencies.

Real-world finance and weather data

To validate our method to real-world finance data, we studied public dataset that includes the year-to-year (YOY) growth rates of monthly Consumer Price Index (CPI), monthly OK WTI spot oil price (OIL) and quarterly real Gross Domestic Product (GDP) in the US (Fig. 6a). The CPI and OIL shares the same sampling rate, and standard VAR method detects GC from OIL to CPI (Fig. 6b). We then used TF-CCA to estimate the relationship between OIL and GDP and the relationship between CPI and GDP. The two lag-CC profiles revealed the causality from OIL to GDP (Fig. 6c), and a weaker causality from CPI to GDP (Fig. 6d). Next, we used partial CCA to assess the effect of CPI on the relationship between OIL and GDP, and found that causality relationship between OIL and GDP diminished (Fig. 6e). In conclusion, MF-TFCCA and GC methods reveal a similar causal chain from OIL to CPI to GDP. It is worth noting that MF-TFCCA can be used on a large range of lags without heavy computational burden, making it possible to detect potential long-term driving forces in the data.

Next, we applied MF-TFCCA to real-world climate data, which consisted of monthly and annual precipitation and temperature time series collected in Central Park, New York City from January 1869 to May 2024 (Fig. 7a). The goal here is to infer the causal relationship between monthly precipitation and monthly temperature, monthly precipitation and yearly temperature, and monthly temperature and yearly precipitation. Interestingly, there was only weak or non-significant GC between rainfall and temperature after removing the linear trend in the time series (Fig. 7b). By examining the MF time series, however, MF-TFCCA discovered a directed information flow between monthly rainfall and yearly temperature time series (Fig. 7c), but not between monthly temperature and yearly rainfall time series (Fig. 7d), suggesting that the detection result of true causal relationship of two random variables may depend on their relative temporal resolution.

Concurrent hippocampal-neocortical recordings

Advances in multi-region multimodal neural recordings (e.g., concurrent EEG-fMRI, or concurrent electrophysiology and calcium imaging, or concurrent electrophysiology and multifiber photometry) have enabled us to assess functional interactions between brain areas. Identifying directed information flow between upstream and downstream structures of the brain may help reveal state-dependent neural mechanisms.

We considered concurrent hippocampal-neocortical recordings that were publicly available^32^ (Abadchi et al., 2020), which consisted of mouse hippocampal CA1 electrophysiological activity such as LFPs (<300 Hz) and multiunit activity (MUA, >300 Hz), as well as widefield calcium imaging (WFCI) activity of the restrosplenial cortex (RSC), primary visual cortex (V1), and forelimb primary somatosensory cortex (FLS1) (Fig. 8a).

It has been well known that the hippocampal MUA often increases during ripple events, especially in the UP or UP-to-DOWN state^33^ (Buzsaki, 2015). First, we examined the relationship between hippocampal MUA and RSC activity, as the RSC is the closest downstream structure to the hippocampus. The hippocampal MUA was sampled at 2000 Hz, whereas the RSC activity extracted from GCaMP was sampled at 100 Hz, resulting a sampling frequency ratio 20:1 (Fig. 8b). Based on the previously detected hippocampal ripple events during non-rapid-eye-movement (NREM) sleep (n=753), we conducted STFT on the filtered hippocampal ripple band (100–250 Hz) activity using a 40-ms moving window with 30-ms overlapping window. Comparing the lag-CC profile and peak timing before, during and after the detected hippocampal ripple events (Fig. 8c, i–iii), we detected a bidirectional information flow between the RSC and hippocampus (HPC), with a slightly stronger strength in the RSC→HPC direction before and during hippocampal ripples. Additionally, the CC was significantly greater than the surrogate data (Supplementary Fig. 10a) and multiple driving frequencies were discovered (Supplementary Fig. 10b). The CC strength was the highest after the ripple events (Fig. 8c–iv), suggesting a possible cortico-hippocampal-cortical information loop covering before, during, and after ripple events. Among all tested trials, we sorted their CC scores (Fig. 8d) and separated the top 25% and bottom 25% trials (Fig. 8e). We further examined the trial averages of RSC activity and found the top 25% CC trials had higher amplitude (Fig. 8f). Next, we repeated the analysis for V1 and FLS1 and compared their lag-CC profiles during ripples. Similar to the bidirectional (but stronger RSC→HPC) information flow, we also found a bidirectional (but stronger V1→HPC) V1-HPC information flow but with a weaker strength than RSC-HPC (Fig. 8g). In contrast, the strength between FLS1 and HPC was the weakest (nearly indistinguishable from the surrogate data, Supplementary Fig. 10c). Therefore, the hippocampal-cortical information flow during ripple events varies according to the downstream cortical structures, which is consistent with the finding based on large-scale electrophysiological findings^34^ (Nitzan et al., 2022). We also repeated the analysis between the hippocampal LFP and RSC activity and found similar trends (Supplementary Fig. 10d).

DISCUSSION

Identifying GC among MF time series has been studied in the field of statistics and econometrics^11,13,15,17^ (Anderson et al., 2016; Ghysels et al., 2016; Tang et al., 2019; Shojaie and Fox, 2022). However, the current literature is limited to the discussion on the time domain or linearly coupled systems. In many physical and biological systems, directed information flow is associated with oscillatory frequency and it is more informative to interpret GC in the frequency domain and identify the driving frequency for the unidirectional or bidirectional GC system. Here we propose a new analysis pipeline to infer directed spectral information flow between MF time series. Our proposed MF-TFCCA is a nonparametric method for analyzing MF time series, which distinguishes itself from other parametric structural VAR methods in several aspects. First, parametric models require a data stationarity assumption, and often encounter problems in estimation biases (due to a limited sample size), model mismatch or identifiability; in contrast, MF-TFCCA requires no strong statistical assumption and is computationally efficient. Second, MF-TFCCA can not only recover the direction of information flow, but also reveal uni- or bidirectional driving frequencies and their relative strengths. Third, MF-TFCCA can potentially detect a wide range of nonlinear couplings between MF time series. Additionally, it can detect GC among more than two time series based on partialization of CCA. Therefore, our proposed method can serve as a hypothesis testing strategy for SGC during exploratory data analysis.

In MF-TFCCA estimation procedure, we convert the HF signal into a time-frequency representation before running CCA with the LF signals. Our choice of STFT as the time-frequency representation can be modified based on preference. Other types of time-frequency representations can be also considered, such as Wigner-Ville distribution (WVD), continuous wavelet spectrum (CWT), Stockwell-transformation (ST)^35^ (Stockwell et al., 1996), short-time linear canonical transform (STLCT)^36^ (Wei and Hu, 2021), and model-based parametric spectrum^37^ (Ba et al., 2014). Specifically, WVD has several advantages over the STFT in that it can produce a high degree of resolution in both time and frequency, especially for non-stationary signals, and can distinguish between closely spaced frequency components. However, computation of STFT is more efficient due to its $[eqn]$ complexity. To characterize potentially nonlinear interactions between random variables, it is also possible to generalize CCA to nonlinear CCA or combine CCA with different preprocessing^25,38^ (Zhuang et al., 2020; Wang et al., 2020). Importantly, MF-TFCCA is nonparametric and is highly computationally efficient compared to other parametric methods such as MF-VAR.

We have validated our method using extensive computer simulations (Supplementary Table 1) and compared it with the standard MF-VAR method on benchmark experiments (Supplementary Table 2). In nearly all tested examples, MF-TFCCA was capable of identifying true directed information flow even in the presence of nonlinear coupling. However, various factors may contribute to false positive detection. Specifically, the effectiveness of our method is influenced by multiple factors, such as signal duration, driving frequency (in relation to the sampling frequency) and driving strength. For instance, constrained by the low sampling rate of LF signals, it requires that spectral causality is preserved with the LF range. Therefore, our method will perform better when the driving frequency is also low, or when the interference of other signals is at HF range. In the case of chain system, our method can somewhat mediate the intermediate variable using partial CCA, but is still prone to detect the false positive depending on the role of intermediate variable. In the case of complex nonlinear systems with highly periodic and correlated variables (e.g., coupled Rössler-Lorenz system), MF-TFCCA also identifies some false-positives of causal relationship. A systematic investigation of other nonlinear CCA methods for reducing the false-positive detection is beyond the scope of current paper and will be the subject of future investigation. Common challenges remain for detecting the causal relationship, such as strong autocorrelation, time delays or periodicity, and nonlinearity^39^ (Runge et al., 2019). Combing our method with parametric approaches to detrend and reduce the autocorrelation in the time series could be a potential solution. Our method can be easily integrated with the existing causal inference framework for causal hypothesis testing.

We have also validated our method in a wide range of real-world data. The use of CCA for identifying inter-area brain communications from neural signals has a long history^25^ (Zhuang et al., 2020). Our proposed MF-TFCCA has replicated the GC results between the S1 and ACC LFP signals during pain experiments (Guo et al., 2020). Furthermore, MF-TFCCA can be viewed as a multimodal CCA method, which can analyze mixed-modal neural signals, such as concurrent electrophysiology and calcium imaging^40,41^ (Wei et al., 2020; Patel et al., 2020), or concurrent EEG and fMRI^42^ (Correa et al., 2010). Bidirectional hippocampal-neocortical communications are well known for sleep-dependent memory consolidation especially around hippocampal SWRs^33,43,44^ (Buzsaki, 2015; Sirota et al., 2003; Pedrosa et al., 2022). When assessing hippocampal MUA and neocortical calcium imaging activity, we found a bidirectional yet asymmetric information flow with a stronger strength in the cortico-hippocampal direction. The degree of information flow from cortical areas to the hippocampus varied during ripple events, with the highest strength in the intermediate hippocampal downstream structure (RSC→HPC), and the weaker strength in sensory cortical areas (V1→HPC and FLS1→HPC). Additionally, we found differences in the RSC→HPC information flow before, during, and after hippocampal ripple events, where the cortico-hippocampal information was stronger before and during ripples. Concurrent mouse RSC and hippocampal LFP recordings during sleep have shown that the RSC displays a pre-ripple activation associated with slow and fast oscillations, and putative RSC inhibitory and excitatory neurons increase and decrease firing activities, respectively immediately after ripples^45^ (Opalka et al., 2020). Recent rodent data have also shown the complexity of bidirectional hippocampal-neocortical communications during UP-DOWN states in NREM sleep^46^ (Feliciano-Ramos et al., 2023), and their couplings are brain-state dependent^47^ (Nitzan et al., 2020). It should be emphasized that the concept of Granger causality needs to be cautioned in self-organized systems such as the brain, where causes in such systems are often circular or multidirectional. Finally, the successes of MF-TFCCA in detecting directed causality in finance and weather data further demonstrated its broader applications and ability to discover weak directed GC relationship.

In general, inferring statistical dependency between MF variables presents a challenge for many well-studied estimation problems for evenly sampled data, including GC estimation^13^ (Ghysels et al., 2016), regression and forecasting^48,49^ (Ghysels et al., 2006; Bai et al., 2013), and stochastic control^50^ (Sinopoli et al., 2004). The information of directed Granger causality may provide an intuitive strategy to forecast low-frequency variables (e.g., weekly) at higher frequencies (e.g., daily). In the context of inferring generalized SGC, compared to the standard MF-VAR method, our proposed analysis framework provides an exploratory and computationally efficient nonparametric approach to quantify directed information flow of MF time series in the presence of complex and nonlinear interactions. In the future, we expect to see more investigations to new explorations in physical and biological sciences.

METHODS

Canonical correlation analysis (CCA).

Without the loss of generality, Given two multivariate random variables $[eqn]$ and $[eqn]$ (where the dimensionality of $[eqn]$ and $[eqn]$ may differ), CCA finds pairs of dimensions, such that the correlation between the projected activity onto these dimensions is maximally correlated: $[eqn]$ , where the vectors $[eqn]$ and $[eqn]$ have respective dimensions. CCA can be estimated by solving a generalized eigenvalue decomposition problem:

$$[eqn]$$

where $[eqn]$ denotes the sample covariance matrix for specific random variables, and $[eqn]$ corresponds to the canonical correlation (CC). The largest squared CC, $[eqn]$ , is used to measure the “set overlap.” The CC also represents the association between the set of $[eqn]$ and the set of $[eqn]$ after the within-set correlations have been removed, and CCA coefficients can be interpreted similar to the coefficients of principle components. CCA has also a natural link to multivariate regression. If the goal is to predict a target variable $[eqn]$ given an input vector $[eqn]$ with a minimum squared error (MSE); let $[eqn]$ , the optimal regression coefficient $[eqn]$ is given by cross-correlation of two whitened variables $[eqn]$ and $[eqn]$ ^51^ (Jendoubi and Streimmer, 2019). A high value of coefficients in $[eqn]$ indicates the relative importance of individual components in $[eqn]$ . In the special case when $[eqn]$ is univariate, the MSE reduces to variance-scaled coefficient of determination: $[eqn]$ . To impose a sparsity constraint onto the CCA coefficients (or equivalently the regression coefficient $[eqn]$ ), regularized CCA algorithms have been designed^25^ (Zhuang et al., 2020).

In time series analysis, traditional CCA assumes that two signals have zero lag. To generalize CCA to lagged CCA where two signals are temporally lagged, let $[eqn]$ , we can reformulate the above equation as follows

$$[eqn]$$

Partialization of CCA (also known as partial CCA) generalizes CCA and estimates a pair of linear projections onto a low dimensional space, where the correlation between two multivariate variables is maximized after eliminating the influence of a third variable^52,53^ (Rao, 1969; Mukuta and Harada, 2014). Let $[eqn]$ and $[eqn]$ be the target variables and $[eqn]$ be the third variable, partial CCA is computed by solving a generalized eigenvalue decomposition problem as follows

$$[eqn]$$

where $[eqn]$ . In CCA or partial CCA, it is common to standardize all the variables.

MF-TFCCA for three or more time series

To extend the analysis framework to the setting with multivariate $[eqn]$ (HF process) and univariate $[eqn]$ (LF process). For the sake of simple discussion, we assume that $[eqn]$ is two-dimensional, where $[eqn]$ and $[eqn]$ denote two signals have the same temporal resolution. In this case, we apply STFT separately to $[eqn]$ and $[eqn]$ , and then apply partial CCA to estimate the CC between $[eqn]$ and $[eqn]$ conditional on $[eqn]$ , or the CC between $[eqn]$ and $[eqn]$ conditional on $[eqn]$ .

In another setting with a univariate variable $[eqn]$ (HF process) and a multivariate variable $[eqn]$ (LF process). For the sake of simple discussion, we assume that $[eqn]$ is two-dimensional, where $[eqn]$ and $[eqn]$ denote two signals have the same temporal resolution. In this case, we apply STFT to $[eqn]$ and then apply partial CCA to assess the CC between $[eqn]$ and $[eqn]$ conditional on $[eqn]$ , or the CC between $[eqn]$ and $[eqn]$ conditional on $[eqn]$ .

Lagged cross-correlation

We use the MATLAB function “xcorr” to compute time-lagged normalized cross-correlation function between two time series $[eqn]$ and $[eqn]$ in the same sampling rate. In the case of MF time series, we compute the correlation between $[eqn]$ and $[eqn]$ for different lags, with the temporal resolution of the lag determined by the larger sampling frequency.

VAR and SGC estimate

For an arbitrary multivariate time series with the same sampling rate, we can model it using a $[eqn]$ -order linear $[eqn]$ system for the demeaned multivariate random variable $[eqn]$

$$[eqn]$$

where $[eqn]$ denotes $[eqn]$ -dimensional white Gaussian noise, $[eqn]$ denotes the $[eqn]$ -by- $[eqn]$ coefficient matrix for the $[eqn]$ -th lag. If the $[eqn]$ process is stable, then all the roots of the reverse characteristic polynomial are bigger than 1 in terms of the Euclidean norm (i.e., outside the unit circle). When $[eqn]$ is an even number, the $[eqn]$ can maximally produce $[eqn]$ oscillatory frequencies. Once $[eqn]$ are known or fully identified, the SGC can be analytically computed^10^ (Barnett and Seth, 2014). Note that when more than two variables are involved, the conditional SGC will be computed. In the special case of two variables, conditional GC and marginal GC are equivalent. The $[eqn]$ parameters can be identified from the least-squared estimation, following a model order selection for $[eqn]$ . Between time series with the same sampling frequency, we can compute GC statistics in either time or frequency domain using an established MATLAB toolbox (www.sussex.ac.uk/sackler/mvgc/).

MF-VAR and GC estimate

For multivariate time series with MF, we adopted the established MF-VAR model^13–15^ (Ghysels et al., 2016; Gotz and Hecq, 2019; Anderson et al., 2016). Specifically, MF-VAR estimates GC by fitting a $[eqn]$ model of a stacked vector $[eqn]$ :

$$[eqn]$$

where $[eqn]$ denotes the HF variable and $[eqn]$ denotes the LF variable. $[eqn]$ is the ratio between the high sampling frequency and the low sampling frequency. For implementation, we modified an existing MATLAB toolbox (http://www2.kobe-u.ac.jp/~motegi/Matlab_Codes.html) for computing single-trial MF-VAR into a multi-trial version.

Surrogate data

To generate surrogate data for time series, we conducted the Fourier transform of time series, and randomly shuffled the phase while keeping the magnitude unchanged; we then conducted the inverse Fourier transform using the magnitude and randomized phase to reconstruct the surrogate time series. The procedure was repeated 100–500 times.

Computer Simulations of VAR Processes

Bivariate unidirectional systems.

The two random variables were first generated from a bivariate VAR(4) model

$$[eqn]$$

The diagonal elements of $[eqn]$ determine the oscillatory frequencies of $[eqn]$ and $[eqn]$ , and the off-diagonal elements determine causality between $[eqn]$ and $[eqn]$ . In generating the raw signals with a sampling rate of 200 Hz and a GC relationship from $[eqn]$ to $[eqn]$ , we assumed $[eqn]$ had two oscillatory frequencies ( $[eqn]$ and $[eqn]$ ), and $[eqn]$ had two oscillatory frequencies ( $[eqn]$ and $[eqn]$ ), and there was SGC at 80 Hz and 4 Hz at $[eqn]$ direction (Supplementary Fig. 2a). The VAR(4) coefficient matrices were set as follows:

$$[eqn]$$

where $[eqn]$ , and $[eqn]$ and $[eqn]$ denotes the sampling interval. Note that the coefficients $[eqn]$ in $[eqn]$ in $[eqn]$ and $[eqn]$ in $[eqn]$ control the degree of causality.

In the bivariate system $[eqn]$ (Supplementary Fig. 2h), we assumed that $[eqn]$ and $[eqn]$ had the same oscillations as in the $[eqn]$ system. The $[eqn]$ coefficient matrices were set as follows:

$$[eqn]$$

Note that although the off-diagonal coefficients were smaller than the $[eqn]$ coefficients, the driving force in the $[eqn]$ system could be higher than that in the $[eqn]$ system, since the causality is also related with the spectrum of the transmitter^54^ (Stokes and Purdon, 2017).

In computer simulations, we down-sampled $[eqn]$ by a factor $[eqn]$ (e.g., $[eqn]$ ) while keeping the sampling rate of $[eqn]$ intact ( $[eqn]$ ). To avoid aliasing in down-sampling, we conducted a low-pass filtering before decimation (see Supplementary Note). To compute the STFT of $[eqn]$ , we chose a window size of 200 ms and the overlap for the moving window as 200–25=175 ms (i.e., window size minus sampling interval of $[eqn]$ ); consequently, the matrix $[eqn]$ has the same length in the time-axis as time series vector $[eqn]$ .

Bivariate bidirectional system.

We generated two random variables with a bidirectional causal relationship $[eqn]$ using a bivariate $[eqn]$ model. To discriminate the two driving directions, the VAR order was large as 41 to create the lags for each direction (Fig. 2a). We assumed that $[eqn]$ had an oscillatory frequency at $[eqn]$ and $[eqn]$ had an oscillatory frequency at $[eqn]$ . Both oscillatory frequencies were in the low frequency band in order to maintain the bidirectional causality after down-sampling of $[eqn]$ was sampled at 200 Hz and $[eqn]$ was sampled at 40 Hz. The first two VAR(41) coefficient matrices control the oscillations of $[eqn]$ and $[eqn]$ :

$$[eqn]$$

where $[eqn]$ and $[eqn]$ . Additionally, coefficient matrices $[eqn]$ to $[eqn]$ control the causality $[eqn]$ with a lag of $[eqn]$ . Coefficient matrices $[eqn]$ to $[eqn]$ control the causality $[eqn]$ with a lag of $[eqn]$ . The elements $[eqn]$ were set to be the $[eqn]$ element of a 10th-order band-pass filter: [0.000, 0.001, −0.014, −0.039, 0.026, 0.098, 0.026, −0.039, −0.014, 0.001, 0.000]. Other elements were all set as 0.

Trivariate chain system.

We generated three random variables using a trivariate VAR(4) model

$$[eqn]$$

In the chain system $[eqn]$ , we assumed that $[eqn]$ had two resonant frequencies at $[eqn]$ and $[eqn]$ had an oscillatory frequency at $[eqn]$ , and $[eqn]$ had two oscillatory frequencies at $[eqn]$ and $[eqn]$ (Fig. 4a). $[eqn]$ and $[eqn]$ were sampled at 200 Hz, whereas $[eqn]$ was sampled at 40 Hz. Further, we assumed that $[eqn]$ strongly and $[eqn]$ moderately using the following VAR(4) setup:

$$[eqn]$$

where $[eqn]$ . When all three variables were known, the system are completely observed; if only $[eqn]$ and $[eqn]$ can be accessed, the system is a partially observed chain system.

For the chain system $[eqn]$ , time-series data were generated using the same VAR model (Supplementary Fig. 5a), except for the intermediate variable sampled at 40 Hz instead of 200 Hz.

Trivariate parallel system.

In the parallel system, we assumed that $[eqn]$ weakly and $[eqn]$ strongly using the following new VAR(4) coefficients (Supplementary Fig. 6a):

$$[eqn]$$

where $[eqn]$ and $[eqn]$ had the same oscillatory frequencies as in the chain system. In the case of parallel system, the results of fully observed and partially observed systems are similar, since the missing variable will not affect the GC between the other two variables. However, the missing variable ( $[eqn]$ or $[eqn]$ ) may interfere with specific frequencies of $[eqn]$ , making it more challenging to identify SGC.

Trivariate VAR(3) model.

We also studied a generic trivariate VAR system (Supplementary Fig. 1d) in the previous study^54^ (Stokes and Purdon, 2017)

$$[eqn]$$

where $[eqn]$ .

The parameters are set as follows: $[eqn]$ , and $[eqn]$ (i.e., sampling frequency 120 Hz).

Phase amplitude coupling (PAC) system.

We used a two-step procedure to generate a pair of signals with PAC. In step 1, we generated raw signals $[eqn]$ and $[eqn]$ using a linear bidirectional causal system as described in Fig. 2a. In step 2, we viewed $[eqn]$ as the phase signal with a resonating frequency at $[eqn]$ Hz and generated a modulated signal with its amplitude coupling with the phase of $[eqn]$ as follows^55^ (Munia and Kviyente, 2019):

$$[eqn]$$

With $[eqn]$ and $[eqn]$ ; which thereby yielded a nonlinear PAC coupling between $[eqn]$ and $[eqn]$ (Supplementary Fig. 3a).

Two-species Logistic Model

The nonlinear two-species Logistic model was described in the following difference equations^30^ (Ma et al., 2014),

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are self-regulation parameters; $[eqn]$ and $[eqn]$ are two coupling constants. In the unidirectional case (Supplementary Fig. 7a), we set the coupling constants as $[eqn]$ and $[eqn]$ ; in the bidirectional case (Supplementary Fig. 7e), we set $[eqn]$ and $[eqn]$ . We used uniformly distributed random numbers in $[eqn]$ as the system’s initial conditions and discarded the initial 100 time points after the transient dynamics. To create MF time series, we subsequently down-sampled $[eqn]$ by a factor of 5.

Coupled Rössler-Lorenz System

The Rössler attractor and Lorenz attractor are two deterministic chaotic systems. To create a complex nonlinear coupled system, we used the same computer simulation setting as REF^30^ (Ma et al., 2014), in which a 3-dimensional nonlinear Rössler system $[eqn]$ drives a 3-dimensional nonlinear Lorenz system $[eqn]$ (Supplementary Fig. 8a), which is described by the following differential equations

$$[eqn]$$

where $[eqn]$ is a timescale constant and $[eqn]$ denotes the strength of unidirectional coupling. The initial values of the system were randomly chosen and the original sampling frequency for both systems was 100 Hz. To create MF time series, we subsequently down-sampled $[eqn]$ by a factor of 5.

Finance Data

The financial time series consists of the year-over-year (YOY) growth rates of monthly Consumer Price Index (CPI), monthly OK WTI spot oil price (OIL) and quarterly real Gross Domestic Product (GDP) in the United States. The CPI, OIL and GDP data are published by the U.S. Department of Labor, the Energy Information Administration, and the Bureau of Economic Analysis, respectively. These data were collected over the course of 30 years (ranging from January 1987 to January 2016), consisting of 360 monthly observations and 120 quarterly observations. The monthly CPI and OIL time series share the same sampling rate, whereas quarterly GDP time series has the lower sampling rate. The YOY growth rates for each value were computed as follow:

$$[eqn]$$

Climate Data

The climate data consist of average monthly and annual precipitation and temperatures (starting from January 1869 to May 2024) at the Central Park of New York City, which are publicly available at the National Weather Service (https://www.weather.gov/media/okx/Climate/CentralPark/monthlyannualprecip.pdf and https://www.weather.gov/media/okx/Climate/CentralPark/monthlyannualtemp.pdf). The time series consist of 1865 monthly observations and 155 yearly observations. Given this dataset, the estimation goal is to infer the causal relationships between monthly precipitation and monthly temperature, monthly precipitation and yearly temperature, and monthly temperature and yearly precipitation time series.

Concurrent Rat Cortical LFP Recordings

In a prior study^31^ (Guo et al., 2020), noxious pain stimuli were used for freely exploring male adult Sprague-Dale rats (250–300 g) in a plastic chamber of size 38×20×25 cm^3^ on top of a mesh table. In thermal pain experiments, a blue (473 nm diode-pumped solid-state) laser with 250 mW intensity was delivered to the rat’s hindpaw. The laser stimulation was delivered in repeated trials (25–40) during 30–45 min. The stimulation was terminated by animal’s paw withdrawal. To produce chronic inflammatory pain, 0.075–0.1 ml of Complete Freund’s adjuvant (CFA) was suspended in an oil-saline (1:1) emulsion, and injected subcutaneously into the plantar aspect of the hindpaw opposite to the paw that was stimulated by a blue laser. Namely, only a unilateral inflammation was induced, and nociceptive stimuli were delivered to the opposite paw of the injured foot. All experimental studies were performed in accordance with the National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals to ensure minimal animal use and discomfort, and were approved by the New York University Grossman School of Medicine Institutional Animal Care and Use Committee (IACUC).

We used silicon probes (NeuroNexus) with 3D printed drive to record multi-channel (up to 64 channels) neural activities from the rat ACC and S1 areas simultaneously. The probe implant was on the contralateral side of the paw that received noxious stimulation. The Plexon (Dallas, TX) data acquisition system was used to record in vivo extracellular neural signals at a sampling rate of 40 kHz. Local field potential (LFP) signals were further preprocessed and down-sampled to 200 Hz. Details have been presented elsewhere (Guo et al., 2020)^31^. For the purpose of benchmark comparison, we selected some preprocessed LFP data from a single CFA rat during 250 mW laser stimulation trials. Each laser stimulation trial lasted 1 second. Data used in our analysis are also publicly available (https://github.com/QiqiXian/MF-TFCCA).

Multimodal Mouse Hippocampal-Neocortical Recordings

In a prior study^32^ (Abadchi et al., 2020), concurrent hippocampal-neocortical recordings were collected from head-restrained mice during natural sleep and urethane anesthesia. Specifically, wide-field mesoscale neocortical activity was imaged based on optical imaging and voltage-sensitive dye (VSD) in the animal’s right hemisphere, and in vivo electrophysiological recordings, including LFPs and multi-unit activity (MUA), were simultaneously measured in the animal’s ipsilateral hippocampus. The VSD imaging employed a CCD camera and recorded 12-bit images at 100 Hz sampling rate. Sensory stimulation was used to determine the coordinates for the primary sensory areas, whereas the relative locations of associated areas were estimated using stereotaxic coordinates. The animal protocols were approved by the University of Lethbridge Animal Care Committee and were in accordance with guidelines set forth by the Canadian Council for Animal Care.

Muscle electromyogram (EMG) recordings and hippocampal delta-to-theta band power ratio to score the sleep state. VSD imaging signals were denoised, filtered and preprocessed to obtain the ΔF/F activity relative to the baseline. Hippocampal LFP was first down-sampled to 2 kHz, sharp-wave ripples (SWRs) were identified using ripple-band filter combined with an established threshold criterion. Condition on the detected hippocampal SWRs, that patterns of activity in neocortical regions were differentially modulated around hippocampal SWRs. Three regions of interest (ROIs) in neocortical areas, including the retrosplenial cortex (RSC), primary visual cortex (V1), and forelimb primary somatosensory cortex (FLS1), were identified and selected in our current analysis. Neural recordings were preprocessed in advance and are publicly available (https://doi.org/10.5061/dryad.qnk98sfbb). However, the original raw recordings were not available, limiting the scope of our analyses.

Supplementary Material

1