Cross-population amplitude coupling in high-dimensional oscillatory neural time series

Neural oscillations have long been considered important markers of interaction across brain regions, yet identifying coordinated oscillatory activity from high-dimensional multiple-electrode recordings remains challenging. We sought to quantify time-varying covariation of oscillatory amplitudes across two brain regions, during a memory task, based on local field potentials recorded from 96 electrodes in each region. We extended Canonical Correlation Analysis (CCA) to multiple time series through the cross-correlation of latent time series. This, however, introduces a large number of possible lead-lag cross-correlations across the two regions. To manage that high dimensionality, we developed rigorous statistical procedures aimed at finding a small number of dominant lead-lag effects. The method correctly identified ground truth structure in realistic simulation-based settings. When we used it to analyze local field potentials recorded from the prefrontal cortex and visual area V4, we obtained highly plausible results. The new statistical methodology could also be applied to other slowly varying high-dimensional time series.

1 Introduction

Contemporary technologies for recording neural activity can produce multiple time series in each of two or more brain regions simultaneously (e.g., Jun et al., 2017; Steinmetz et al., 2018), enabling identification of cross-regional interactions relevant to cognitive processes that support behavior. One such process is working memory (Miller et al., 2018; Schmidt et al., 2019). During a variety of experimental working memory tasks, strong neural oscillations in the beta range (16–30 Hz) have been observed, leading to the proposal that beta oscillations serve to coordinate activity across regions (Miller et al., 2018). To identify coordinated oscillatory activity, many statistical methods based on coherence, or phase coupling, have been developed, studied, and applied (Belluscio et al., 2012; Klein et al., 2020; Ombao and Pinto, 2022; Urban et al., 2023). As shown in Figure 1, however, oscillatory amplitudes (amplitude envelopes) also vary substantially and may well exhibit statistical dependence across regions. Our aim is to detect whether these slowly-varying amplitude envelopes tend to fluctuate together, and to pinpoint when during a task such coupling occurs. This form of amplitude-based coordination reflects a different aspect of neural communication than phase-based measures and requires statistical tools that can accommodate smooth, non-oscillatory signal components. In this study, we develop methods for analyzing groups of slowly-varying multiple time series, and apply them to uncover coordinated amplitude activity between brain regions.

Figure 1

Figure 1. Bandpass-filtered LFP amplitude envelopes. The top panels display beta-range filtered LFPs (gray) together with their amplitude envelopes (black) from a pair of electrodes in V4 (left) and PFC (right) on a single trial, during the delay period of a working memory task. (The filtering is described in Section 3.2.1) The bottom panels display amplitude envelopes of active electrodes in the two brain regions for two trials, those for one trial in black and those for the other trial in gray. The beta amplitude envelopes show a consistent temporal pattern within each region, for both trials, and, at 400 ms (red vertical line), the two groups of amplitude envelopes from the two brain regions illustrate correlated behavior across the two trials in the sense that the black amplitude curves for the first trial are elevated, compared with the gray for the second trial, in both V4 and PFC. This correlation, however, varies across time.

A common strategy for quantifying inter-areal association is to reduce multichannel recordings within each brain region to a small number of region-level summaries and then compute pairwise association measures between regions. For oscillatory neural signals, this is often done by computing amplitude envelope correlations after extracting instantaneous amplitudes via the Hilbert transform, based on region-level summaries (Bruns et al., 2000; Zamm et al., 2018; Raghavan et al., 2024). Constructing these summaries typically requires additional preprocessing, such as selecting a single representative channel per region (Bruns et al., 2000), applying spatial filters (Zamm et al., 2018), or designating channels based on expert judgment (Raghavan et al., 2024).

A similar dimension reduction principle underlies methods developed for non-oscillatory recordings such as fMRI. Time Delay (TD) analysis (Mitra et al., 2015, 2018; Raut et al., 2019) averages voxel-level signals within each region to obtain a single regional time-series. It then compares pairs of regional time series by measuring how much one signal is delayed relative to the other using cross-covariance, resulting in a single lag value for each region pair that can be further analyzed using PCA. Cyclicity analysis (Shahsavarani et al., 2020; Abraham et al., 2024) instead aggregates directionality measures across voxel-level interactions to produce a signed scalar measure that represents a region-level summary that captures region-level lead-lag relationships.

These approaches collapse multichannel activity into signals that mix task-relevant and task-irrelevant processes, which reduces sensitivity to between-region coupling and obscures the subpopulations that carry the interaction. We instead developed a method to learn low-dimensional and time-localized components directly from multichannel data, without spatial priors. By maximizing between-region dependence (relative to within-region variance), our approach increases sensitivity to inter-areal communication.

The data that motivated this work are local field potentials (LFPs; Buzsáki et al., 2012; Pesaran et al., 2018), recorded from prefrontal cortex (PFC) and visual area V4 during repeated trials of a working memory task (Khanna et al., 2020). PFC is generally considered to exert control aimed at areas involved in perceptual processing (Miller and Cohen, 2001), such as V4, which is active during the retention of higher-order visual information (e.g., color and shape) and attention to visual objects (Orban, 2008; Fries et al., 2001). Both regions exhibit strong oscillatory activity during visual working memory tasks, and their interaction, particularly through oscillatory coupling, has been implicated in supporting working memory function (Sarnthein et al., 1998; Liebe et al., 2012; Miller et al., 2018; Gregoriou et al., 2009).

The data consist of 96 LFP time series for each brain region (192 time series in total), offering millisecond temporal precision, and are relatively precise spatially, as well. These signals reflect neural activity across large populations of neurons within each region and exhibit substantial spatial correlation. To summarize the cross-regional interactions more succinctly, dimensionality reduction becomes a natural and effective step (Gallego-Carracedo et al., 2022). If we were to analyze the data at a single time point, we would have repeated pairs of 96-dimensional observations, for which canonical correlation analysis (CCA) would be a standard tool to assess statistical dependence. Here we build on a latent variable formulation of CCA known as probabilistic CCA (pCCA), and extend it to handle temporally structured data, as illustrated in Figure 2b. In the extended model described in Section 2.2, each of the two high-dimensional time series is driven by a latent univariate time series, with the correlation between the two latent time series representing the statistical dependence we seek to identify (see Figure 3). Theorem 3.3 demonstrates that this model yields estimates equivalent to those from a generalization of classical CCA, called multiset CCA (Kettenring, 1971). This model-based framework allows for flexible modeling of covariance structure and facilitates principled high-dimensional inference.

Figure 2

Figure 2. Graphical representation of pCCA models. (a) Model of Bach and Jordan (2005)), where X₁ and X₂ are random vectors and Z is a random variable. (b) A variation on (a) that facilitates the extension to the case when X₁ and X₂ are multivariate time series and (Z₁, Z₂) is a bivariate time series.

Figure 3

Figure 3. Extended pCCA model for two multivariate time series $X_1^{(t)}$ and $X_2^{(s)}$, t, s = 1, …, T. (a) Dynamic associations between vectors $X_1^{(t)}$ and $X_2^{(s)}$ are summarized by the dynamic associations between their associated 1-dimensional latent variables $Z_1^{(t)}$ and $Z_2^{(s)}$ represented by their cross-precision matrix Ω₁₂. (b) When a significant cross-precision entry is identified, for example, the red star in the expanded view of Ω₁₂, its coordinates and distance from the diagonal indicate at what time in the experiment connectivity between two brain areas occurs, and at what lead or lag. Here the red star is in the upper diagonal of Ω₁₂, which means that, at this particular time, region 1 leads region 2, or Z₁→Z₂ in short (a non-zero entry in the lower diagonal would mean Z₂→Z₁). We represent this association by the red arrow on the right-most plot, with a lag of two units of time for illustration.

Related latent variable methods include Gaussian Process Factor Analysis (GPFA; Lakshmanan et al., 2015; Semedo et al., 2020; Yu et al., 2009) and state-space models (Goris et al., 2014). GPFA assumes that latent time series follow a low-dimensional Gaussian process, while state-space models typically impose autoregressive dynamics. Although these methods were originally devised to capture shared low-dimensional structure within a single brain region, they can be embedded in a two-step approach to discover cross-regional interaction: (1) apply those methods separately to each brain region, and (2) compute cross-correlations between the resulting latent trajectories. While common in neuroscience, this strategy may lack statistical power when within-region activity does not align with cross-regional communication. For example, in recordings from macaque visual areas V1 and V2, Semedo et al. (2019)) found that communication occurred in a dedicated “communication subspace” distinct from the dominant within-region activity. Ignoring this distinction can lead to underestimating cross-regional coupling and reduced sensitivity (Semedo et al., 2020).

To address this, alternative methods have been proposed. Rodu et al. (2018)) extended a kernel version of CCA (KCCA; Hardoon et al., 2004) to develop Dynamic Kernel CCA (DKCCA), designed to detect time-varying connectivity between multivariate signals.

In this study, we introduce a more flexible alternative by leaving the latent covariance unspecified, allowing for a fully data-driven estimate of cross-regional correlations. This yields a partial correlation graph representation of neural interactions. Given the high dimensionality of the possible partial correlations across time lags, we adopt a sparse estimation framework (Friedman et al., 2008), resulting in our method: Latent Dynamic analysis via Sparse banded graphs (LaDynS). While we applied LaDynS to LFP recordings, it could also be used to analyze other slowly-varying multidimensional neural measurements made during repeated-trial cognitive tasks.

In Section 2, after defining the model and proving that it produces a time series generalization of CCA, we give details of the fitting procedure and discuss inference using clustered contiguous significant association based on a de-biased estimate of the precision matrix (Jankova et al., 2015) together with control of false discovery rate (Benjamini and Hochberg, 1995). We also use a state-space formulation to impose a local stationarity assumption (Ombao and Pinto, 2022), pointing out how this can produce a time-varying latent Granger causality assessment. Simulations in Section 3.2 show that our implementation of LaDynS is able to correctly identify the timing of interactions when applied to artificial data designed to be similar to those we analyzed, even in scenarios where existing methods exhibit low statistical power. The simulations support data-analytic results presented in Section 3.4, which show significant cross-regional interaction at roughly 400 ms after presentation of the stimulus. This timing coincides with the delay in visual cortical response seen after the onset of a visual stimulus in mice (Chen et al., 2022), humans (Del Cul et al., 2007; Yang et al., 2019), and macaques (Supèr et al., 2001; Schmolesky et al., 1998). The results reinforce the idea that PFC and V4 are involved, together, in working memory. We add a discussion in Section 4.

2 Methods

We begin by reviewing and reformulating probabilistic CCA (pCCA) in Section 2.1, which we then generalize to time series in Section 2.2 to yield a dynamic version of pCCA. We define the LaDynS model based on the log-likelihood function in Theorem 2.1. We go over the choice of regularization parameters in Section 2.3.1 and the fitting algorithm in Section 2.3.2. We discuss statistical inference in Section 2.4.

2.1 Probabilistic CCA for two random vectors

Given two random vectors $X_1\in {ℝ}^{d_1}$ and $X_2\in {ℝ}^{d_2}$, canonical correlation analysis (CCA; Hotelling, 1992) finds the sets of weights $w_1\in {ℝ}^{d_1}$ and $w_2\in {ℝ}^{d_2}$ that maximize Pearson's correlation between linear combinations $w_1^{\top }{X}_1$ and $w_2^{\top }{X}_2$. This can be rewritten as

$\begin{array}{ll}{\sigma }_{cc}={max}_{w_k,k=1,2:w_k^{\top }{\Sigma }_{kk}{w}_k=1}\left|w_1^{\top }{\Sigma }_{12}{w}_2\right|,& (1)\end{array}$

where Σ_kk = Var(X_k) is the covariance matrix of X_k, k = 1, 2, and Σ₁₂ = Cov(X₁, X₂) the cross-covariance matrix between X₁ and X₂. The sample estimator ${\widehat{\sigma }}_{cc}$ is obtained by replacing Σ_kk and Σ₁₂ with their sample analogs ${\bar{\Sigma }}_{kk}$ and ${\bar{\Sigma }}_{12}$ respectively. The maximizing weights ŵ_k and linear combinations ${Ẑ}_k={ŵ}_k^{\top }{X}_k$ are referred to as the canonical weights and canonical variables, respectively.

Probabilistic CCA assumes that X₁ and X₂ are driven by a common one dimensional latent variable Z:

$\begin{array}{ll}\begin{array}{l}{X}_k|Z={\mu }_k+Z·{\beta }_k+{ϵ}_k,k=1,2,\\ Z~N(0,1),\end{array}& (2)\end{array}$

where ${\mu }_k\in {ℝ}^{d_k}$ and ${\beta }_k\in {ℝ}^{d_k}$ are mean vectors and factor loadings, respectively, and ${ϵ}_k\stackrel{\mathrm{\text{indep}}}{~}MVN(0,{\Phi }_k)$ (Bach and Jordan, 2005). Figure 2a depicts the dependence of X₁ and X₂ on Z. The parameters in Equations 1, 2 both yield the same estimate of σ_cc; see Theorem 3.1 for more details.

Next, we introduce an alternative pCCA extension that assigns distinct latent variables for X₁ and X₂, as depicted in Figure 2b. Specifically, we assume that

$\begin{array}{ll}{X}_k|Z_k={\mu }_k+Z_k·{\beta }_k+{ϵ}_k& (3)\end{array}$

where ${\mu }_k\in {ℝ}^{d_k}$, ${\beta }_k\in {ℝ}^{d_k}$ and ${ϵ}_k\stackrel{\mathrm{\text{indep}}}{~}MVN(0,{\Phi }_k)$ are defined as in Equation 2, and (Z₁, Z₂) are bivariate normally distributed:

$\begin{array}{ll}\left(\begin{array}{c}{Z}_1\\ Z_2\end{array}\right)~MVN\left(\left(\begin{array}{c}0\\ 0\end{array}\right),\left(\begin{array}{cc}1& {\sigma }_{12}\\ {\sigma }_{12}& 1\end{array}\right)\right).& (4)\end{array}$

Like the original method, the alternative pCCA yields the same estimate of σ_cc. We prove this equivalence in Theorem 3.2.

2.2 Probabilistic CCA for two time series of random vectors

Suppose now that we are interested in the correlation dynamics between two times series of random vectors $X_1^{(t)}\in {ℝ}^{d_1}$ and $X_2^{(t)}\in {ℝ}^{d_2}$, t = 1, 2, …, T. For each time t, we use Equation 3 to model the dependence of $X_k^{(t)}$ on its associated latent variable $Z_k^{(t)}$:

$\begin{array}{ll}{X}_k^{(t)}|Z_k^{(t)}={\mu }_k^{(t)}+{\beta }_k^{(t)}·Z_k^{(t)}+{ϵ}_k^{(t)},k=1,2,& (5)\end{array}$

where ${\mu }_k^{(t)}$, ${\beta }_k^{(t)}$ and ${ϵ}_k^{(t)}\stackrel{\mathrm{\text{indep}}}{~}MVN(0,{\Phi }_k^{(t)})$ are defined as in Equation 3.

Then for each t we could define a parameter ${\sigma }_{12}^{(t)}$ as in Equation 4 to capture population-level association between $X_1^{(t)}$ and $X_2^{(t)}$ at t. But because we are also interested in lagged associations between $X_1^{(t)}$ and $X_2^{(s)}$ for s≠t, we replace the bivariate Equation 4 for $Z_1^{(t)}$ and $Z_2^{(t)}$ for a given t by a global model for all t = 1, …, T:

$\begin{array}{ll}{\left({\left(Z_1^{(t)}\right)}_{t=1,\dots ,T},{\left(Z_2^{\left(t\right)}\right)}_{t=1,\dots ,T}\right)}^{\top }~MVN\left(0,\Sigma \right),diag\left(\Sigma \right)=\text{1},& (6)\end{array}$

where Σ captures all simultaneous and lagged associations within and between the two time series jointly. Equation 3a illustrates the dependence structure of this model. to highlight the auto-correlations Σ_kk, k = 1, 2, within and cross-correlations Σ₁₂ between the times series, and denote each element of Σ_ij by Σ_ij(t, s), (t, s)∈[T]². Then Σ₁₂(t, t) has the same interpretation as σ₁₂ in EqXXXX and when t≠s, Σ₁₂(t, s) measures the lead-lag association between the time series at latency |t−s|. We decompose Σ and its inverse Ω as

$\begin{array}{ll}\Sigma =\left(\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{12}^{\top }& {\Sigma }_{22}\end{array}\right)\Omega =\left(\begin{array}{cc}{\Omega }_{11}& {\Omega }_{12}\\ {\Omega }_{12}^{\top }& {\Omega }_{22}\end{array}\right)& (7)\end{array}$

to highlight the auto-correlations Σ₁₁ and Σ₂₂ within and cross-correlations Σ₁₂ between the time series, and denote by ${\Sigma }_{kl}^{(t,s)}$, (t, s)∈[T]², the elements of Σ_kl. Then ${\Sigma }_{12}^{(t,t)}$ for some fixed t has the same interpretation as σ₁₂ in Equation 2. Further, $-{\Omega }_{12}^{(t,s)}/\sqrt{{\Omega }_{11}^{(t,t)}{\Omega }_{22}^{(s,s)}}$ is the partial correlation between the two latent time series at times t and s. Thus, when an element of Ω₁₂ is non-null, depicted by the red star in the expanded display in Figure 3b, its coordinates (t, s) and distance (t−s) from the diagonal indicate at what time in the trial a connectivity happens between two time series, and at what lead or lag, respectively. In our neuroscience application, they represent the timing of connections and the direction of information flow between two brain regions. As is common with latent variable models, the signs of the latent covariance Σ and therefore the latent precision Ω are not identifiable. This non-identifiability is benign for inference on the directionality in lead-lag connectivity, which is encoded in the support of Ω₁₂, rather than by their signs (see Figure 3b). The directionality remains interpretable even though the signs of the partial correlations are unidentifiable. Therefore, in the following sections, we focus on the estimation and inference for the magnitude of associations, rather than their sign.

As the equivalence between a non-distributional method (CCA) and its probabilistic representation (pCCA), there exists a similar connection between the multiset generalization of CCA introduced by Kettenring (1971)) and the dynamic pCCA model in Equations 5, 6. Multiset CCA applied to 2T random vectors $\left\{X_1^{(t)},X_2^{(t)}:t=1,\dots ,T\right\}$ finds weights $\left\{w_1^{(t)},w_2^{(t)}:t=1,\dots ,T\right\}$ that maximize a notion of correlation among linear combinations $\left\{w_1^{(t)\top }{X}_1^{(t)},w_2^{(t)\top }{X}_2^{(t)}:t=1,\dots ,T\right\}$. In particular, the GENVAR extension minimizes the generalized variance of these linear combinations, defined as the determinant of their correlation matrix (Wilks, 1932), which we refer to as the canonical correlation matrix:

$\begin{array}{ll}\begin{array}{l}{\widehat{w}}_1^{(1)},\dots ,{\widehat{w}}_2^{(T)}=\arg {\min }_{w_1^{(1)},\dots ,w_2^{(T)}}\det (\overline{\text{Var}}[{\left(w_1^{(t)\top }{X}_1^{(t)}\right)}_{t=1,\dots ,T},\hfill \\ {\left(w_2^{(t)\top }{X}_2^{(t)}\right)}_{t=1,\dots ,T}])\hfill \end{array}& (8)\end{array}$

where $\overline{Var}$ denotes the sample variance-covariance matrix and the weights $w_k^{(t)}$ are scaled so that every diagonal entry of the matrix is 1. In Theorem 3.3 we show that the maximum likelihood estimators from the dynamic pCCA model recover the same canonical weights $\left\{w_1^{(1)},\dots ,w_2^{(T)}\right\}$ and the canonical correlation matrix $\widehat{\Sigma }$ as the GENVAR formulation of multiset CCA. Under this equivalence, the MLEs minimize

$\begin{array}{ll}logdet(\Sigma )+tr({\Sigma }^{-1}\bar{\Sigma }),& (9)\end{array}$

An objective that motivates the Latent Dynamic Analysis via Sparse Banded Graphs (LaDynS) developed in the next section. The objective above does not involve any component of the observation model Equation 5; thus, the estimators do not depend on a Gaussian assumption for $X_k^{(t)}\mid Z_k^{(t)}$. They only require approximate normality of the latent factors $Z_k^{(t)}$, which is often plausible since each $Z_k^{(t)}$ is a weighted sum of many coordinates of $X_k^{(t)}$.

2.3 Latent dynamic analysis via sparse banded graphs (LaDynS)

Our goal is to estimate the association dynamics between two multivariate time series of length T using the covariance matrix Σ of their associated latent time series in Equation 6. However, the prohibitive number of parameters in Σ means its estimation is prone to errors, especially when T is large. We reduce their number by regularizing Ω = Σ⁻¹, rewriting logdet(Σ) = logdet(Ω⁻¹) = −logdet(Ω), and assuming that Ω has the banded structure depicted in Figure 4.

Figure 4

Figure 4. The elements of Ω_kk, k = 1, 2, and Ω₁₂ are set to zero outside of the gray bands of widths (1+2d_auto) and (1+2d_cross), respectively.

Definition 0.0.1 (LaDynS). Given N simultaneously recorded pairs of multivariate time series {_{X1[n], X2[n]}n = 1, …, N}, and a 2T×2T sparsity matrix Λ with element ${\Lambda }_{kl}^{(t,s)}$ regularizing $\left|{\Omega }_{kl}^{(t,s)}\right|$, k, l = 1, 2, LaDynS finds weights $\left\{{ŵ}_k^{(t)},t=1,2,\dots ,T,k=1,2\right\}$ and precision matrix $\widehat{\Omega }$ that minimize the penalized negative log-likelihood:

$\begin{array}{ll}-logdet(\Omega )+tr(\Omega \bar{\Sigma })+\left|\right|\Lambda \odot \Omega |{|}_1,& (10)\end{array}$

where $\bar{\Sigma }=\overline{Var}\left[w_1^{(1)\top }{X}_1^{(1)},\dots ,w_2^{(T)\top }{X}_2^{(T)}\right]$ satisfies $diag(\bar{\Sigma })=\text{1}$, ⊙ denotes the Hadamard product operator such that (A⊙B)_ij = A_ij×B_ij, $\left|\right|A|{|}_1=\sum _{i,j}|A_{ij}|$, and

${\Lambda }_{kl}^{(t,s)}=\{\begin{array}{ll}{\lambda }_{\text{cross}},\hfill & k\ne l\text{and}0<|t-s|\le d_{\text{cross}},\hfill \\ {\lambda }_{\text{auto}},\hfill & k=l\text{and}0<|t-s|\le d_{\text{auto}},\hfill \\ {\lambda }_{\text{diag}},\hfill & t=s,\hfill \\ \infty ,\hfill & \text{otherwise},\hfill \end{array}$

which constrains auto-precision and cross-precision elements within a specified range.

In our neuroscience application, in particular, it is reasonable to assume that lead-lag relationships occur with delay less than temporal bandwidth d_cross, which can be determined by the maximal transmission time in synaptic connections between two brain regions under study. We thus set ${\Lambda }_{12}^{(t,s)}=\infty$ when |t−s|>d_cross to force the corresponding cross-precision elements to zero and thus impose a banded structure on Ω₁₂. We apply sparsity constraint ${\Lambda }_{12}^{(t,s)}={\lambda }_{\mathrm{\text{cross}}}>0$ on the remaining off-diagonals of Ω₁₂ to focus our discovery of sparse dominant associations and reduce the effective parameter size. We proceed similarly with the auto-precision matrices Ω₁₁ and Ω₂₂, using penalty λ_auto and temporal bandwidth d_auto. Unless domain knowledge is available, we recommend that d_auto be set to the largest significant auto-correlation across all observed time series $X_{k,i}^{(t)}$, k = 1, 2, i = 1, …, N, and impose no further sparsity (λ_auto = 0) unless there is reason to expect it.

Notice that, for simplicity, we grouped the elements of Λ into diagonal and off-diagonal elements and assigned the same penalties, λ_cross, λ_auto = 0, and λ_diag, within each group.

Again, LaDynS estimates the magnitude of the associations, but their signs remain non-identifiable. Accordingly, our results in Section 3 report and display the absolute values of the estimated variance and precision elements.

2.3.1 Choosing regularization parameters

In graphical LASSO (gLASSO) problems, where the aim is to recover correct partial correlation graphs, penalties are often chosen to minimize the predictive risk (Shao, 1993; Zou et al., 2007; Tibshirani et al., 2012). Our aim is different: only the partial cross-precision matrix Ω₁₂ is of substantive interest, and because minimizing the predictive risk does not select models consistently (Shao, 1993; Zhu and Cribben, 2018) and may thus fail to retrieve non-zero elements of Ω₁₂, we choose instead a value of λ_cross that controls the number of false cross-precision discoveries. We proceed by permuting the observed time series in one brain region to create a synthetic dataset that contains no cross-region correlation, then applying LaDynS to that data for a range of values of λ_cross and recording the resulting number of significant partial correlation estimates, which are necessarily spurious. We use the smallest λ_cross that yields fewer false discoveries than a chosen threshold. We expect this regularization to make similarly few false discoveries on experimental data.

Finally, if $\widehat{\Sigma }$ cannot be inverted, as is the case for the band-pass filtered experimental data we analyze in Section 3, we penalize its diagonal by λ_diag>0. We explain the specific calibration We used the analyzed datasets and studied the properties in Section 3.2.2.

2.3.2 Fitting LaDynS using coordinate descent

Equation 10 is not a convex function of the weights and precision elements (although it is not impossible that it may be for some particular Σ), and its convex relaxation is unknown, so it is difficult to find its global minimum. The coordinate descent algorithm described below finds a minimum, possibly local, since it may be sensitive to the choice of initial parameter values. However, the simulation in Supplementary Section S2.2 suggests that the algorithm's sensitivity to initial values has little impact on its ability to recover the correct connectivity.

Assuming that all canonical weights $w_k^{(t)}$ are fixed, Equation 10 reduces to the gLASSO problem:

$\begin{array}{ll}{\displaystyle \underset{\Omega }{argmin}}-logdet(\Omega )+tr(\Omega \bar{\Sigma })+\left|\right|\Lambda \odot \Omega |{|}_1,& (11)\end{array}$

which we can solve efficiently using a number of existing algorithms; here we use the P-gLASSO algorithm of Mazumder and Hastie (2012)). Then, assuming that all parameters are fixed but a single weight $w_k^{(t)}$, Equation 10 can be rearranged as the linear problem:

$\begin{array}{ll}{\displaystyle \underset{w_k^{(t)}}{argmin}}{\sum }_{(l,s)\ne (k,t)}{w}_k^{(t)\top }\overline{Cov}\left[X_k^{(t)},X_l^{(s)}\right]w_l^{(s)}{\Omega }_{kl}^{(t,s)}& (12)\end{array}$

$\begin{array}{l}\text{ s.t. }{w}_k^{(t)\top }\overline{Var}\left[X_k^{(t)}\right]w_k^{(t)}=1,\end{array}$

for which an analytical solution is available. That is, our algorithm alternates between updating Ω and the weights $w_k^{(t)}$ until the objective function in Equation 10 converges. Its computational cost is inexpensive: a single iteration on our cluster server (with 11 Intel(R) Xeon(R) CPU 2.90 GHz processors) took, on average, less than 0.8 s, applied to the experimental data in Section 3.4. A single fit on the same data took 47 iterations for around 33.57 s until the objective function converged at threshold 0.001.

Algorithm 1 presents a high-level pseudocode sketch of the coordinate descent algorithm. See Supplementary Section S2.1 for details and Python package ladyns on github.com/HeejongBong/ladyns.

Algorithm 1

Algorithm 1. Coordinate descent algorithm to fit LaDynS.

2.4 Inference for associations between two vector time series

Let $\widehat{\Omega }$ and ${ŵ}_k^{(t)}$, t = 1, …, T, k = 1, 2, be the LaDynS estimates of canonical precision matrix and canonical weights, and $\bar{\Sigma }=\overline{Var}\left[{ŵ}_1^{(1)\top }{X}_1^{(1)},\dots ,{ŵ}_2^{(T)\top }{X}_2^{(T)}\right]$ be the empirical covariance of the estimated latent variables, defined in Equation 10. Note that $\widehat{\Omega }\ne {\bar{\Sigma }}^{-1}$ since $\widehat{\Omega }$ is constrained to be sparse. Based on these estimates, we want to identify the non-zero partial cross-correlations in Ω₁₂, which correspond to the epochs of association between the two time series.

Formal inference methods for Ω based on its LaDynS estimate (Equation 10) are not available, but because LaDynS reduces to graphical LASSO (gLASSO) when the weights $w_k^{(t)}$ in Equation 11 are fixed, we co-opt gLASSO inference methods. Specifically, Jankova et al. (2015)) suggested de-sparsifying the gLASSO estimate $\widehat{\Omega }$ according to

$\begin{array}{ll}\tilde{\Omega }=2\widehat{\Omega }-\widehat{\Omega }(\bar{\Sigma }+{\lambda }_{\text{diag}}{I}_T)\widehat{\Omega },& (13)\end{array}$

and proved that, under mild assumptions and as n → ∞, each entry of $\tilde{\Omega }$ satisfies the Central Limit Theorem with center the true precision Ω:

$\begin{array}{ll}\forall (t,s),\frac{({\tilde{\Omega }}_{12}^{(t,s)}-{\Omega }_{12}^{(t,s)})}{\sqrt{Var\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]}}\stackrel{d}{\to }N(0,1).& (14)\end{array}$

The assumptions include independence across experimental trials, irrepresentability condition, and bounded eigenvalues of Ω [see Assumptions (A1) and (A2) in Jankova et al. (2015))]. While some of these conditions may be strong and difficult to verify in practice, they are commonly imposed in the literature on sparse graph estimation with ℓ₁ regularization. In Supplementary Section S3.2, we examine the asymptotic normality of the de-sparsified LaDynS estimate using simulation.

We applied this result to the de-sparsified LaDynS estimate of Ω, even though we do not quite have a gLASSO setup, and we verified by simulation that its elements are indeed approximately normal in Section 3.2.2. Jankova et al. (2015)) also proposed an estimator of $Var\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$ that is likely to be downward-biased in our framework since estimating the canonical weights $w_k^{(t)}$ induces extra randomness. Instead, we use the bootstrap estimate $\widehat{Var}\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$ described at the end of this section. We then rely on Equation 14 to obtain p-values

$\begin{array}{ll}{p}^{(t,s)}=2-2\Phi \left(\left|{\tilde{\Omega }}_{12}^{(t,s)}\right|/\sqrt{\widehat{Var}\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]}\right)& (15)\end{array}$

to test $H_0^{(t,s)}:{\Omega }_{12}^{(t,s)}=0$, for each (t, s)∈[T]² within d_cross of the diagonal of Ω₁₂.

Permutation bootstrap estimate of $Var\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$: a permutation bootstrap sample ${\left\{X_{1\left[n\right]}^{*},X_{2\left[n\right]}^{*}\right\}}_{n=1,\dots ,N}$ is generated by permuting separately {_{X1[n]}n = 1, …, N} and {_{X2[n]}n = 1, …, N}. Hence, applying LaDynS to ${\left\{X_{1\left[n\right]}^{*},X_{2\left[n\right]}^{*}\right\}}_{n=1,\dots ,N}$ yields estimates of canonical precision matrix ${\widehat{\Omega }}^{*}$, canonical weights ${ŵ}_k^{*}(t)$s, empirical covariance of the estimated latent variables ${\bar{\Sigma }}^{*}=\overline{Var}({ŵ}_1^{*(1)\top }{X}_1^{*(1)},\dots ,{ŵ}_2^{*(T)\top }{X}_2^{*(T)})$, and de-sparsified precision matrix estimate ${\tilde{\Omega }}^{*}$ (Equation 13) under the global null hypothesis of no correlated activity. Repeating the bootstrap simulation B times produces B bootstrap values ${\widehat{\Omega }}^b$, ${ŵ}_k^{b(t)}$, ${\bar{\Sigma }}^b$, and ${\tilde{\Omega }}^b$, b = 1, …, B. We then estimate $Var\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$ with $\widehat{Var}\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$, the sample standard error of ${\left\{{\tilde{\Omega }}_{12}^{b,(t,s)}\right\}}_{b=1,2,\dots ,B}$. Notice that $\widehat{Var}\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$ is obtained under the global null $H_0^{(t,s)}:{\Omega }_{12}^{(t,s)}=0$ simultaneously for all (t, s), because it is not trivial to simulate bootstrap data that satisfy a specific $H_0^{(t,s)}$ without assuming that all other elements of Ω₁₂ are also null. We garnered from simulations that $\widehat{Var}\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$ is thus likely to slightly underestimate $Var\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$, which makes for slightly sensitive p-values.

Control of false discoveries: because we perform tests for many entries of Ω₁₂, we cap the false discovery rate

$\begin{array}{ll}\text{FDR}=E\left[\text{FDP}\right],\text{where}\text{FDP}=\frac{\#\left\{\text{falsely discovered entries}\right\}}{\#\left\{\text{discovered entries}\right\}\vee 1}& (16)\end{array}$

below a pre-specified level α_BH using the procedure of (Benjamini and Hochberg, 1995; BH). To proceed, let $p^{\left[1\right]}\le \cdots \le p^{\left[n_{\mathrm{\text{roi}}}\right]}$ denote the ordered permutation bootstrap p-values p^{(t, s)} that correspond to n_roi cross-precision elements in the region of interest. Then, we find the maximum k_BH satisfying $p^{\left[k_{\mathrm{\text{BH}}}\right]}\le \frac{k_{\mathrm{\text{BH}}}}{n_{\mathrm{\text{roi}}}}{\alpha }_{\mathrm{\text{BH}}}$ and reject $H_0^{(t,s)}$ with p^{(t, s)} smaller than $\frac{k_{\mathrm{\text{BH}}}}{n_{\mathrm{\text{roi}}}}{\alpha }_{\mathrm{\text{BH}}}$. The FDR guarantee is established by Benjamini and Hochberg (1995)) as long as the p^{(t, s)}'s are independent and valid p-values.

Cluster-wise inference by excursion test: as a further safeguard against falsely detecting correlated activity between brain areas, we obtain p-values for each identified connectivity epoch using the excursion test of Ventura et al. (2005)), as follows.

After computing p-values p^{(t, s)} for each precision entry (Equation 15), we define a cluster as any contiguous set of entries whose p-values fall below $\frac{k_{\text{BH}}}{n_{\text{roi}}}{\alpha }_{\text{BH}}$. No additional cluster-forming threshold or minimum cluster size is imposed. For each cluster k identified by the BH procedure, we calculate the test statistic:

$\begin{array}{ll}{T}_k:=-2{\displaystyle \sum _{(t,s)\in \text{cluster}k}}log{p}^{(t,s)},& (17)\end{array}$

which is reminiscent of Fisher's method Fisher (1925)) for testing the global significance of multiple hypotheses. Large values of T_k provide evidence against cross-area connectivity in cluster k, so we calculate the corresponding p-value as $\underset{T_k}{\overset{\infty }{\int }}{f}_{T_{max}}(t)dt$, where f_Tmax is the null distribution of $T_{max}:=\underset{j}{max}{T}_j$ under the global null hypothesis of no connectivity anywhere. We use f_Tmax rather than the respective null distributions of each T_k to control the family-wise type I error rather than the type I error for each cluster. We approximate f_Tmax by the previous permutation bootstrap: for each permuted dataset b = 1, …, B, we estimate the cross-precision matrix and corresponding p-values, identify all clusters of p-values below $\frac{k_{\mathrm{\text{BH}}}}{n_{\mathrm{\text{roi}}}}{\alpha }_{\mathrm{\text{BH}}}$, calculate the corresponding test statistics in Equation 17, and let $T_{max}^b$ be their maximum. The B values $T_{max}^b$ are samples from f_Tmax, which we use to approximate the p-value for cluster k by the sampling proportion:

$\begin{array}{l}\frac{1}{B}{\displaystyle \sum _{b=1}^B}I(T_{max}^b\ge T_k).\end{array}$

2.5 Locally stationary state-space model and local Granger causality

Our model in Equations 5, 6 can be formulated as a state-space model by rewriting the joint multivariate Gaussian model for the latent vectors in Equation 6 as the set of all conditional distributions

$\begin{array}{ll}\begin{array}{c}{Z}_1^{(t)}={\displaystyle \sum _{s=1}^{d_{\text{auto}}}}{\alpha }_{11,s}^{(t)}{Z}_1^{(t-s)}+{\displaystyle \sum _{s=1}^{d_{\text{cross}}}}{\alpha }_{12,s}^{(t)}{Z}_2^{(t-s)}+{\eta }_1^{(t)},\\ Z_2^{(t)}={\displaystyle \sum _{s=1}^{d_{\text{auto}}}}{\alpha }_{22,s}^{(t)}{Z}_2^{(t-s)}+{\displaystyle \sum _{s=1}^{d_{\text{cross}}}}{\alpha }_{21,s}^{(t)}{Z}_1^{(t-s)}+{\eta }_2^{(t)},\\ \end{array}& (18)\end{array}$

where d_auto and d_cross are the maximal time delays in within-region and cross-region connections, ${\eta }_k^{(t)}$, k = 1, 2, are independent $N(0,{\varphi }_k^{(t)})$ random variables, and the ${\alpha }_{kl,s}^{(t)}$'s are vector auto-regressive coefficient parameters for the auto-correlation within region if k = l, k = 1, 2, and cross-correlation between regions if k≠l with time lag s.

This state-space formulation is convenient to impose local stationarity on the latent time series, which we do by fitting stationary state-space models in moving windows of time. Local stationarity is justified because the functional connectivity within and between brain regions changes relatively slowly over time. The state-space formulation is also convenient for calculating the Granger causality between regions: Z₂ is said to Granger-cause Z₁ at time t if some ${\alpha }_{12,s}^{(t)}$ are non-zero (Ombao and Pinto, 2022; and conversely if some ${\alpha }_{21,s}^{(t)}$ are non-zero). The coefficient of partial determination (partial R²) between $(Z_2^{(t-d_{\text{cross}})},\dots ,Z_2^{(t-1)})$ and $Z_1^{(t)}$, conditional on $Z_1^{(t-d_{\text{auto}})},\dots ,Z_1^{(t-1)}$, calculated as

$\begin{array}{ll}{R}_{2\to 1}^2(t)=1-\frac{Var\left[\text{residual of Regression 1}\right]}{Var\left[\text{residual of Regression 2}\right]},& (19)\end{array}$

where

$\begin{array}{ll}\begin{array}{l}\text{Regression 1}:Z_1^{(t)}~Z_1^{(t-d_{\text{auto}})}+\cdots +Z_1^{(t-1)}+Z_2^{(t-d_{\text{cross}})}\\ +\cdots +Z_2^{(t-1)},\\ \text{Regression 2}:Z_1^{(t)}~Z_1^{(t-d_{\text{auto}})}+\cdots +Z_1^{(t-1)},\end{array}& (20)\end{array}$

may therefore be considered a test statistic for local Granger causality at time t. To allow a physiologically meaningful minimum connection time τ₁ from brain regions 2 to 1, we can also replace the second regression with

$\begin{array}{l}\begin{array}{l}\text{Regression 2}:Z_1^{(t)}~Z_1^{(t-d_{\text{auto}})}+\cdots +Z_1^{(t-1)}+Z_2^{(t-{\tau }_1+1)}\\ +\cdots +Z_2^{(t-1)}\\ +\text{1}\left\{{\tau }_2<d_{\text{cross}}\right\}(Z_2^{(t-d_{\text{cross}})}+\cdots +Z_2^{(t-{\tau }_2-1)}),\\ \end{array}\end{array}$

where τ₂ is the maximum connection time from brain regions 2 to 1, τ₁ ≤ τ₂ ≤ d_cross. We set τ₂ = d_cross by default, unless there is reason to consider shorter connection times. A plug-in estimator of $R_{2\to 1}^2(t)$ is easily obtained from the estimated covariance matrix of $(Z_1^{(t)},\dots ,Z_1^{(t-d_{\text{auto}})},Z_2^{(t-1)},\dots ,Z_2^{(t-d_{\text{cross}})}),$ without actually running the regressions (Anderson-Sprecher, 1994).

Autocorrelations in the latent time series can inflate R² values. We therefore test the statistical significance of $R_{2\to 1}^2(t)$ (or $R_{1\to 2}^2(t)$) by comparing its observed value to its null distribution, obtained by repeatedly permuting the trials in one region and calculating $R_{2\to 1}^2(t)$ in the permuted data. We used 2,000 permutations in Section 3.4. The permuted data satisfy the null hypothesis of no cross-region connection and exhibit the same autocorrelation structure as the original latent time series.

3 Results

We have introduced LaDynS to estimate the dynamic connectivity between two or more multivariate time series, and we proposed inference procedures to identify when connectivity is statistically significant. We apply LaDynS to experimental data in Section 3.4, but first we present the theoretical results on the equivalence between generative pCCA models and model-free CCA methods. In particular, Theorem 3.3 establishes an equivalence between the GENVAR version of multi-set CCA (Kettenring, 1971) and maximum likelihood applied to our dynamic pCCA. Next we examine its performance on simulated data that have properties similar to the experimental data. The data are simulated from the shared oscillatory driver model described in Section 3.2, which is unrelated to the LaDynS model. Supplementary Section S3.2 contains results based on simulated data that are consistent with the LaDynS model. In Section 3.3, we compare the performance of LaDynS to other existing methods. The reproducible code scripts for the simulations and experimental data analyses are provided at github.com/HeejongBong/ladyns.

3.1 Equivalence between dynamic pCCA model and multiset CCA

We first revisit the connection between classical CCA and the probabilistic CCA formulation of Bach and Jordan (2005)). Let ${\left\{(X_{1\left[n\right]},X_{2\left[n\right]})\right\}}_{n=1}^N$ be independent observations drawn from the joint distribution in Equation 2, and denote the maximum likelihood estimates by $({\widehat{\beta }}_1,{\widehat{\beta }}_2)$. The result below restates Theorem 2 of Bach and Jordan (2005)).

Theorem 0.0.2 (Bach and Jordan, 2005, Theorem 2). The maximum likelihood estimators (MLEs) $({\widehat{\beta }}_1,{\widehat{\beta }}_2)$ in Equation 2 based on N observed vector pairs {_{X1[n], X2[n]}n = 1, 2, …, N} are equivalent to the CCA solution $({ŵ}_1,{ŵ}_2,{\widehat{\sigma }}_{cc})$ in Equation 1 according to:

$\begin{array}{ll}{\widehat{\beta }}_k={\bar{\Sigma }}_{kk}{ŵ}_k{m}_k,\text{where }{m}_1{m}_2={\widehat{\sigma }}_{cc}\text{and }\left|m_k\right|\le 1,k=1,2.& (21)\end{array}$

Theorem 3.1 proves that the original CCA setting and the generative pCCA model both yield the same estimate of σ_cc.

We now state an equivalence similar to Theorem 3.1 between the original CCA and the alternative pCCA model.

Theorem 0.0.3. The MLEs $({\widehat{\beta }}_1,{\widehat{\beta }}_2,{\widehat{\sigma }}_{12})$ in Equations 3, 4 based on N observed vector pairs {_{X1, [n], X2, [n]}n = 1, 2, …, N} are equivalent to the CCA solution $({ŵ}_1,{ŵ}_2,{\widehat{\sigma }}_{cc})$ according to:

$\begin{array}{ll}{\widehat{\beta }}_k={\bar{\Sigma }}_{kk}{ŵ}_k{m}_k,\text{where }{m}_1{m}_2{\widehat{\sigma }}_{12}={\widehat{\sigma }}_{cc}\text{and }\left|m_k\right|\le 1,k=1,2.& (22)\end{array}$

In practice we take m₁ = m₂ = 1 out of all possible solutions, because then $Z_k|X_k={ŵ}_k^{\top }{X}_k$ is the canonical variable almost surely, and σ₁₂ = Cov[Z₁, Z₂] equals the canonical correlation σ_cc. This means that σ₁₂ is an interpretable parameter, and one for which inference is simpler than for the canonical correlation in the other model, because, in Equation 2, ${\widehat{\sigma }}_{cc}$ is an indirect function of the maximum-likelihood parameter estimates (see Theorem 3.1). The interpretability property also persists when we extend Equation 4 to capture the lagged association between two vector time series (see Equation 6). Finally, the choice m₁ = m₂ = 1 implies that the MLEs $({\widehat{\beta }}_1,{\widehat{\beta }}_2,{\widehat{\sigma }}_{12})$ do not depend on the Gaussian assumption in Equation 3, an assumption that is questionable if, for example, the X's are positive variables like LFP power envelopes or discrete variables like spike counts. Theorem 3.2 is recovered as a special case of Theorem 3.3 below, which extends the result to multiset CCA for vector time series.

We now derive a similar connection between the multiset generalization of CCA introduced by Kettenring (1971)) and the dynamic pCCA model in Equations 5 6. The proof is provided in Supplementary Section S1.

Theorem 0.0.4. Suppose that ${\widehat{\beta }}_k^{(t)},k=1,2,t=1,\dots ,T$, and $\widehat{\Sigma }$ are the MLE in Equations 5, 6 based on N observed pairs of vector time series $\left\{X_{1\left[n\right]}^{(t)},X_{2\left[n\right]}^{(t)}:t=1,\dots ,T\right\}$, n = 1, …, N. Then, they are equivalent to the GENVAR multiset CCA solution according to:

$\begin{array}{ll}{\widehat{\beta }}_k^{(t)}=\overline{Var}\left[X_k^{(t)}\right]{ŵ}_k^{(t)}{m}_k^{(t)}\text{and }{\widehat{\Sigma }}_{kl}^{(t,s)}=& (23)\end{array}$

$\{\begin{array}{ll}1,\hfill & k=landt=s,\hfill \\ \overline{\text{Var}}\left[{\widehat{Z}}_k^{(t)},{\widehat{Z}}_l^{(s)}\right],\hfill & \text{elsewhere,}\hfill \end{array}$

where the canonical variable is

$\begin{array}{ll}{Ẑ}_k^{(t)}={ŵ}_k^{(t)\top }{X}_k^{(t)}=\frac{1}{m_k^{(t)}}{\widehat{\beta }}_k^{(t)\top }{\overline{Var}}^{-1}\left[X_k^{(t)}\right]X_k^{(t)},& (24)\end{array}$

and $\left|m_k^{(t)}\right|\le 1$ for k = 1, 2 and t∈[T]. Furthermore, if all $m_k^{(t)}=1$, then the MLE minimizes

$\begin{array}{ll}logdet(\Sigma )+tr({\Sigma }^{-1}\bar{\Sigma }),& (25)\end{array}$

where $\overline{\Sigma }=\overline{\text{Var}}[{\left({\beta }_1^{(t)\top }{\overline{\text{Var}}}^{-1}[X_1^{(t)}]X_1^{(t)}\right)}_{t=1,\dots ,T},$ ${\left({\beta }_2^{(t)\top }{\overline{\text{Var}}}^{-1}[X_2^{(t)}]X_2^{(t)}\right)}_{t=1,\dots ,T}]$, with β₁ and β₂ scaled such that $diag\bar{\Sigma }=\text{1}$.

3.2 LaDynS performance on simulated data from a shared oscillatory driver model with time delay

We used a probabilistic model of the inter-areal coherence in local field potentials to generate LFP time series with a dynamic lead-lag relationship between two brain areas. The shared oscillatory driver model with time delay in Ombao and Pinto (2022)) assumes that the coherence between two univariate LFP time series L₁ and L₂ is driven by a univariate latent oscillation L₀ at frequency f₀:

$\begin{array}{ll}\begin{array}{c}{L}_k^{(t)}={\beta }_k·L_0^{(t-{\tau }_k)}+{\eta }_k^{(t)},k=1,2,\end{array}& (26)\end{array}$

where ${\eta }_k^{(t)}$ is regional baseline noise, and τ_k is the lead-lag from L₀ to L_k. This model generalizes the Synaptic-Source-Mixing (SSM) model of Schneider et al. (2021)), which was shown to be capable of modeling the dynamic cross-regional coherence observed in experimental LFP data between frontal and parietal cortices, as well as between LGN and the visual cortex.

For our simulation, we extended Equation 26 to generate two sets of multi-dimensional LFP time series L_k of dimension d_k, k = 1, 2, and we used three latent oscillations L_{0, j}, j = 1, 2, 3, to allow non-stationary lead-lag relationships between the two brain areas. The resulting model was a special case of the Latent Dynamic Factor Analysis model in Bong et al. (2020)):

$\begin{array}{ll}{L}_k^{(t)}={\displaystyle \sum _{j=1}^3}{\beta }_{kj}·L_{0,j}^{(t-{\tau }_{kj})}+{\eta }_k^{(t)},k=1,2,& (27)\end{array}$

where the factor loadings ${\beta }_{kj}\in {ℝ}^{d_k}$ were vectors of dimension d_k = 25. We virtually arranged the 25 components of $L_1^{(t)}$ and $L_2^{(t)}$ on 5 × 5 regular two-dimensional arrays so we could create simulated data exhibiting spatial correlations similar to those in the experimental data in Section 3.4. To do that we let ${\eta }_k^{(t)}$ be spatially correlated Gaussian noise whose temporal power spectrum scales as 1/f^α for α>1. More specifically, for each temporal frequency f, the respective Fourier coefficients ${\widehat{\eta }}_{ki}(f)$ of region k and electrode i has spatial correlation: $Cov\left[{\widehat{\eta }}_{ki}(f),{\widehat{\eta }}_{kj}(f)\right]=f^{-\alpha }exp(-\frac{\text{dis}{\text{t}}^2(i,j)}{2{\sigma }_{\text{spatial}}^2})$, where dist(i, j) was the distance between electrodes i and j on one array, α = 1.4, and σ_spatial = 0.8.

Next, we let the factor loadings β_kj have components ${\beta }_{kji}=\gamma ·exp(-\frac{\text{dis}{\text{t}}^2(i,p_{kj})}{2{\sigma }_{\text{spatial}}^2})$, i = 1, …d_k, where p_kj was a randomly chosen location on an array and γ took values in {0.055, 0.063, 0.071, 0.077, 0.084, 0.089}. This set spanned {0.5, 0.66, 0.84, 1.0, 1.16, 1.34} times the observed signal-to-noise ratio (SNR) in the experimental data in Section 3.4. More specifically, the observed data exhibited a signal-to-noise ratio (defined using signal power) of approximately 0.75. Finally, the three latent oscillations L_{0, j} were set at 18 Hz and designed to induce three epochs of lead-lag relationships between L₁ and L₂ around experimental times 80, 200, and 400 ms; L₁ leads L₂ by 30 ms during the first epoch, and L₁ lagged L₂ by 30 ms during the other two, as depicted in Figure5a.

Figure 5

Figure 5. Output and inference of LaDynS applied to simulated datasets from the shared oscillatory driver model. (a) True cross-precision matrix Ω₁₂ for the connectivity scenario described in Section 3.2 with γ = 0.077. (b) Average over 60 simulation datasets of LaDynS de-sparsified precision estimates ${\tilde{\Omega }}_{12}$. There is a good match to the true Ω₁₂ in (a). (c) Cross-precision estimate ${\widehat{\Omega }}_{12}$ for one simulated dataset. It matches (a) up to random error. (d) Permutation bootstrap p-values for the de-sparsified estimate ${\tilde{\Omega }}_{12}$. (e) Discovered non-zero cross-precision estimates by the BH procedure at nominal FDR 5%. The cluster-wise p-values of the three discovered clusters by the excursion test were all smaller than 0.5%. Panels (a), (b), and (c) share the same color bar in (c).

One simulated dataset consisted of N = 1, 000 trials of 500 ms long d₁ = 25 and d₂ = 25-dimensional time series $\left\{L_1^{(t)}\in {ℝ}^{d_1},L_2^{(t)}\in {ℝ}^{d_2}\right\}$ generated from Equation 27 at sampling frequency 1, 000 Hz. We filtered the simulated LFP recordings using the complex Morlet wavelet at frequency f₀ = 18 Hz and bandwidth 50 ms, the same we applied to the data in Section 3.4, and collected the beta oscillation amplitude envelopes as the absolute values of the filtered signals. The complex Morlet wavelet is a complex sinusoidal with a Gaussian envelope (Gabor, 1946), where the bandwidth refers to the Gaussian standard deviation.

(The filtered amplitude at a given time t is essentially the amplitude that would be obtained using a harmonic regression with cosine and sine damped by the Gaussian kernel centered at t; the wavelet formulation is computationally different and more efficient.)

After downsampling the power envelopes to 100 Hz, we applied LaDynS to the resulting data X₁ and X₂ with T = 50 time points.

3.2.1 LaDynS estimation details

The simulated amplitude time series was very smooth, similar to the experimental data. We thus added the regularizer λ_diag to the diagonal entries of the estimated correlation matrix $\widehat{\Sigma }$ so it could be inverted. Our calibration strategy for λ_diag was as follows. Let $X_{k,i}^{(t)}$ be an observed time series and $S_{k,i}\in {ℝ}^{T\times T}$ be its auto-correlation matrix, k∈{1, 2}, i∈{1, …, d_k}. Band-pass filtering the LFP data induced auto-correlations in $X_{k,i}^{(t)}$, which we should observe in $S_{k,i}^{-1}$, unless $S_{k,i}^{-1}$ is degenerate. We thus took λ_diag to be such that ${(S_{k,i}+{\lambda }_{\mathrm{\text{diag}}}{I}_T)}^{-1}$ displayed the expected auto-correlation. Practically, we chose λ_diag automatically to minimize the ℓ₂ distance between the off-diagonal entries of ${(S_{k,i}+{\lambda }_{\mathrm{\text{diag}}}{I}_T)}^{-1}$ and the band-pass filter kernel-induced auto-correlation, after a scale adjustment, summed over k and i. Penalizing the diagonal introduced inevitable bias to the sparsified and desparsified LaDynS' precision estimates, $\widehat{\Omega }$ and $\tilde{\Omega }$. The other hyperparameters were set to d_auto = d_cross = 10, λ_auto = 0, and λ_cross, the penalty on the cross-correlation elements, was determined as per Section 2.3.1.

3.2.2 Results

The simulated data true cross-precision matrix Ω₁₂ is unknown so we estimated it by simulation. Given a simulated dataset $(X_{1,\left[n\right]}^{(t)},X_{2,\left[n\right]}^{(t)})$, for each trial n = 1, …, N and time t = 1, …, T, we used Equation 24 to recover the true latent factors $(Z_{1,\left[n\right]}^{(t)},Z_{2,\left[n\right]}^{(t)})$ from $(X_{1,\left[n\right]}^{(t)},X_{2,\left[n\right]}^{(t)})$ and known factor loadings $({\beta }_1^{(t)},{\beta }_2^{(t)})$. We then obtained the empirical covariance matrix of the true latent factors: ${\Sigma }^o:=\frac{1}{N}\sum Z_{\left[n\right]}{Z}_{\left[n\right]}^{\top }$, where $Z_{\left[n\right]}=(Z_1^{(1)},\dots ,Z_2^{(T)})$, and calculated the regularized precision matrix ${\Omega }^o:={({\Sigma }^o+{\lambda }_{\text{diag}}I)}^{-1}$. We estimated the true precision matrix with the average of 200 repeats of Ω^o. Figure 5a shows the cross-regional component ${\Omega }_{12}^o$ of this estimate. We see lead-lag relationships between simulated X₁ and X₂ around 80, 200, and 400 ms, as specified in Section 3.2.

Figure 5c displays the LaDynS cross precision estimate ${\widehat{\Omega }}_{12}$ fitted to one dataset simulated under the connectivity scenario depicted in Figure 5a, with connection strength γ = 0.077. Figure 5d shows the permutation bootstrap p-values in Equation 15 (with permutation bootstrap simulation size B = 200) for the entries of the desparsified cross-precision estimate ${\tilde{\Omega }}_{12}$. Small p-values concentrate near the locations of true non-zero cross-precision entries and are otherwise scattered randomly. We then identified the significant connections by first applying the BH procedure with target FDR 5% and next by applying the excursion test at significance level 5% to all discovered clusters (Section 2.4). The significant clusters are plotted in Figure 5e. They match approximately the true clusters in Figure 5a, although they exhibit random variability, as we should expect. To average this random variability out, we estimated ${\tilde{\Omega }}_{12}$ for each of 60 simulated datasets, and plotted their average in Figure 5b. The average LaDynS estimate is a close match to the true cross-precision matrix in Figure 5a.

Figure 6 displays the estimated partial R² obtained from the locally stationary state-space model described in Section 2.5. The pink shaded areas are pointwise 95th percentiles of the null partial R² distribution under the assumption of independence between the two time series. The estimated effects closely reflect the true lead-lag relationships specified in Section 3.2. Specifically, the left panel shows that time series 1 exerts a strong Granger causal influence on time series 2 during the first connectivity epoch around 80 ms, and the right panel shows that time series 2 Granger causes time series 1 during the later connectivity epochs around 200 and 400 ms.

Figure 6

Figure 6. Estimated partial R² from locally stationary state-space model

Figures 5, 6 used connection strength γ = 0.077. We repeated the simulation for a range of γ such that the SNR in the simulated data beta band ranged from zero to twice the estimated SNR of the experimental data in Section 3.4. Instead of showing graphs, we summarized the performance of LaDynS using false cluster discovery and non-discovery rates (Perone Pacifico et al., 2004; FCDR, FCNR) in the estimated precision matrix. We defined a cluster of estimated non-zero precision elements to be falsely discovered if it contained no true effect, and a true cluster was deemed falsely non-discovered if no estimated cluster overlapped with it. We estimated these rates with the corresponding proportions in 60 simulated datasets. They are plotted in Figure 7 against nominal significance levels of the excursion test between 0 and 10%, where the target FDR for the BH procedure was fixed at 5%. The FCDR remained below the nominal level across all tested significance values and connectivity intensities (all or almost all discovered clusters were true clusters) and for large enough γ, the FCNR was zero (all clusters were discovered).

Figure 7

Figure 7. Performance of cluster-wise inference after excursion test. (a) False cluster discovery rates and (b) false cluster non-discovery rate of BH at target FDR 5% followed by excursion test at significance level α∈[0, 0.10] to identify non-zero partial correlations, under the connectivity scenario in Figure 5a, for the simulated range of connectivity intensities.

Perhaps reporting FDR and FNR would have been more informative, but we reported FCDR and FCNR because they were less sensitive to Ω₁₂ being approximated rather than known in this particular example. See Supplementary Section S3.2 for an example where Ω₁₂ is known, so FDR and FNR could be calculated.

Remark. A defining feature of LaDynS is that it is capable of estimating amplitude lead-lag relationships between two high-dimensional oscillatory time-series even if there is no other source of correlations, such as coherence or phase locking, between them. We demonstrated this in Supplementary Section S3.1 using simulated data from the model in this section, modified so it would induce amplitude lead-lag relationships without inducing phase correlations between the simulated pair of time-series.

Checking the assumptions in Equations 14, 15. We assumed that, under the null hypothesis $H_0^{(t,s)}:{\Omega }_{12}^{(t,s)}=0$, the desparsified precision entries ${\tilde{\Omega }}_{12}^{(t,s)}$ were normally distributed and $Var\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$ was well approximated by the permutation bootstrap variance $\widehat{Var}\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$ . To check the former, we compared the empirical distribution of R = 60 repeat estimates ${\tilde{\Omega }}_{12}^{(t,s)}/\sqrt{\widehat{Var}\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]}$ to the standard normal distribution via QQ-plots. Figure 8 shows QQ-plots for three randomly chosen time pairs (t, s) that satisfy ${\Omega }_{12}^{(t,s)}=0$. The good agreement between empirical and theoretical quantiles suggests that the normal assumption is appropriate.

Figure 8

Figure 8. Null distributions of three representative entries of ${\tilde{\Omega }}_{12}^{(t,s)}/\sqrt{\widehat{Var}\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]}$, obtained from R = 60 simulated datasets (Section 3.2) and compared to the standard Gaussian distribution via QQ-plots. The good agreement suggests that the Normal assumption in Equation 13 is appropriate.

Next, for each (t, s) such that $H_0^{(t,s)}:{\Omega }_{12}^{(t,s)}=0$, we checked the validity of the permutation bootstrap variance estimate $\widehat{Var}\left[{\tilde{\Omega }}_{12}^{(t,s)}\right]$ by comparing it to the empirical variance of the R = 60 estimates ${\tilde{\Omega }}_{12}^{(t,s)}$.

Figure 9 shows the Q-Q plot of their ratios against the quantiles of the F(B−1, R−1) distribution, which suggests that the two estimates are equal up to random error.

Figure 9

Figure 9. Standard deviations of desparsified precision elements. F-statistics of ratios between bootstrap and empirical variances for the null entries of Ω₁₂, showing good agreement.

3.3 Comparison to existing methods

In this section, we compare the performance of LaDynS on simulated data generated from the shared oscillatory driver model (see Section 3.2) with three existing methods: Dynamic Kernel Canonical Correlation Analysis (DKCCA; Rodu et al., 2018), Delayed Latents Across Groups (DLAG; Gokcen et al., 2022), and Gaussian Process Factor Analysis (GPFA; Yu et al., 2009). The simulated data are generated under two distinct connectivity scenarios:

• Three-factor scenario: three pairs of latent factors drove the connections, each active in a distinct time epoch, as depicted by the true cross-covariance matrix in Figure 10a, left panel. Each pair had a unidirectional influence with a fixed lead-lag.

• Single-factor scenario: the epochs and lead-lags of the connections were the same as in the three-factor scenario (Figure 11a, left panel), but this time a single pair of latent factors drove all three epochs of connections. As a result, the connection between the pair of latent factors was bidirectional, and the lead-lag changed over time.

Figure 10

Figure 10. Comparison of LaDynS and competing methods applied to one simulated dataset from the shared oscillatory driver model under the three-factor scenario. (a) (left) The true cross-covariance matrix. (center and right) Estimated cross-covariance matrices from LaDynS and DKCCA, respectively. (b) The cross-covariance estimates obtained by DLAG, which pairs the three factors in region 1 with the three factors in region 2, yielding three cross-covariance plots (one for each factor pair). (c–e) The cross-covariance estimates obtained by GPFA, which considers all possible between-region factor pairs among the three factors in region 1 and the three factors in region 2, yielding nine cross-covariance plots.

Figure 11

Figure 11. Comparison of LaDynS and competing methods applied to one simulated dataset from the shared oscillatory driver model under the single-factor scenario. (a) (left) The true cross-covariance matrix. (center and right) Estimated cross-covariance matrices from LaDynS and DKCCA, respectively. (b) The cross-covariance estimates obtained by DLAG, which pairs the three factors in region 1 with the three factors in region 2, yielding three cross-covariance plots (one for each factor pair). (c–e) The cross-covariance estimates obtained by GPFA, which considers all possible between-region factor pairs among the three factors in region 1 and the three factors in region 2, yielding nine cross-covariance plots.

LaDynS models each latent factor time series $Z_k^{(t)}$ to have time-varying factor loadings ${\beta }_k^{(t)}$, accommodating changes in the activated latent factors across time t. Similarly, DKCCA allows each canonical component time series to have time-varying canonical weights. As a result, both LaDynS and DKCCA can capture multiple factor-driven connections using only a single pair of latent factors or canonical components, provided these underlying factors are active in distinct time intervals (e.g., the three-factor scenario). By contrast, DLAG and GPFA use fixed factor loadings over time, and thus require three latent factors per set of observed time series to capture the same three-factor scenario. Accordingly, we estimated the interaction between two sets of simulated time series, $X_1^{(t)}$ and $X_2^{(t)}$, under the following configurations:

• LaDynS: we followed the estimation details in Section 3.2.1. Although LaDynS represents the between-sets interaction via the cross-precision matrix Ω₁₂, we use the cross-covariance matrix Σ₁₂ here to facilitate comparison with the other methods, which estimate cross-covariance matrices only.

• DKCCA: we estimated the cross-covariance between a single pair of canonical components defined by a sliding time window of width 9 (g = 4 in the original notation).

• DLAG: we used three between-group and two within-group factors for each set of time series. In DLAG, each between-group factor in one set is exclusively paired with a counterpart in the other set, yielding three cross-covariance matrix estimates.

• GPFA: because GPFA inherently models a single set of time series, it does not directly provide a cross-connection component. We therefore used a two-step approach: first, we fitted GPFA with three latent factors to each of $X_1^{(t)}$ and $X_2^{(t)}$; then, we computed cross-covariance matrices for all 3 × 3 pairs of latent factors, producing nine cross-covariance estimates.

For all other tuning hyperparameters of DKCCA, DLAG, and GPFA, we used the default settings provided in the example scripts of the respective code packages. For details, we refer the reader to the example code in our ladyns package.

Figure 10 presents the simulation results for the three-factor scenario. Both LaDynS and DLAG successfully identified all three interaction epochs (Figure 10a, center panel, and Figure 10b), whereas GPFA only captured the first two (Figure 10c–e), and DKCCA failed to detect the interaction entirely (Figure 10a, right panel). This outcome for GPFA is unsurprising, given that GPFA was designed primarily to capture dominant within-set dynamics rather than cross-set interactions. Consequently, if cross-regional interactions are relatively weak compared to the within-region variability (as in this three-factor scenario), GPFA may fail to detect all active connections.

Under the single-factor scenario, where the interaction between a single pair of latent factors is sufficiently strong, GPFA was able to recover all three interaction epochs (Figure 11c–e), matching the performance of LaDynS (Figure 11a, center panel). By contrast, DLAG struggled to capture the bidirectional nature of the interaction (Figure 11b), consistent with its assumption that each pair of between-group latent factors has a fixed lead-lag relationship. We additionally ran DLAG with a single between-group latent factor to match the true latent dimensionality in this scenario; however, its qualitative behavior remained unchanged, and DLAG continued to struggle to capture the bidirectional nature of the interaction. DKCCA again did not detect any meaningful interaction in this setting.

Overall, LaDynS was the only method that consistently identified all three interaction epochs in both connectivity scenarios. DKCCA showed weak performance in both scenarios, while DLAG and GPFA had diminished power in at least one of the scenarios. The performances of DLAG and GPFA would further degrade if we did not specify the correct number of latent variables.

3.4 Experimental data analysis

We applied LaDynS to local field potentials (LFPs), collected in the experiment described in Khanna et al. (2020)), from two Utah arrays implanted in a Macaque monkey's prefrontal cortex (PFC) and V4 during a memory-guided saccade task. Each trial of the task started with a monkey fixating its eyes on a dot at the center of the screen. A visual cue was given for 50 ms to indicate a target location, which was randomly chosen from 40 candidate locations (eight directions and five amplitudes) that tiled the display screen. The cue was turned off and the monkey had to remember the target location while maintaining eye fixation for a delay period of 500 ms. After the delay period, the monkey made a saccade to the remembered position, and a liquid reward was provided on successful trials. As in Khanna et al. (2020)), we analyzed the time series during the delay period, based on 3, 000 successful trials. Time t = 0 corresponds to the start of that period. The data are available in Snyder et al. (2022)).

Because beta oscillations are often associated with communication across brain areas (Klein et al., 2020; Miller et al., 2018), we filtered LFP recordings at a beta oscillation frequency 18 Hz and obtained the beta oscillation power envelopes as described in Section 3.2.1. We chose 18 Hz because it was the frequency having the largest power within the range 12–40 Hz (see Supplementary Figure S11). After downsampling the power envelopes to 200 Hz, we applied LaDynS with the regularizer λ_diag on the diagonal of Σ, as in Section 3.2.1, because the filtered data were very smooth.

Figure 12a shows the only epoch of significant contiguous region of the precision matrix identified by our method (FDR at 5%, p < 0.0005 by the excursion test in Section 2.4). This result provides strong evidence that the 18 Hz beta amplitudes in PFC and V4 were correlated (after conditioning on beta amplitudes at all other times and lags) around 400 ms after the start of the delay period.

Figure 12

Figure 12. Experimental data analysis. (a) LaDynS inference. Discovered region of cross-precision using BH at nominal FDR 5% and the excursion test (p < 0.0005). The light gray area shows the region of time considered (one area leading the other by at most 100 ms). The blue blob suggests that activity in V4 preceded that in PFC around 400 ms in the delay period. (b) Estimated partial R² from locally stationary state-space model for (left) PFC → V4 and (right) V4 → PFC. The pink shaded areas indicate the range of values that fall below the 95th percentile under the null hypothesis of independence between V4 and PFC.

To better understand this relationship, we used the estimated latent time series to compute partial R² values under an assumption of locally stationarity, applying the model in Equation 18, Section 2.5. The frequency 18 Hz corresponds to a period of 55 ms, and it is hard to detect non-stationarity at time scales finer than a few periods. We therefore used a moving window of 100 ms (i.e., larger than the period) to calculate partial R², allowing connection delays between V4 and PFC in the τ₁ = 15 to τ₂ = 30 ms range. The partial R² from PFC to V4 and from V4 to PFC are shown in Figure 12b. There are large excursions of R² above the null values in both plots, suggesting that, at around 400 ms post stimulus, the two areas are involved in a bidirectional network: the power of the activity in each area predicts the power of the oscillation in the other, following a short delay (after conditioning on the power at all other times and lags). To see that this is not due solely to the activity passing through the visual stream, in Figure 13, for each of the two latent time series, we plotted the estimated total beta power as a function of time. The power in V4 increases dramatically much earlier than 400 ms, starting around 100 ms and reaching a peak just after 200 ms; it then remains substantial during the remainder of the delay period. The results in Figure 12 can not be explained by those in Figure 13 alone.

Figure 13

Figure 13. Estimated beta powers of electrophysiological activity driven by latent factors in V4 and PFC as functions of time. The summed beta power in the data $X_k^{(t)}$ attributable to the latent time series at time t was estimated by the ℓ₂ norm of the factor loading vector ${\beta }_k^{(t)}$; see Equation 5.

It is also possible to get spatial information from the normalized factor loadings across the electrode arrays, which are displayed in Figure 14 for experimental time 400 ms. The loadings were rescaled by their maximal value across the array. (An animation over the complete timeline is available at github.com/HeejongBong/ladyns.) It is apparent that a relatively small proportion of the electrodes, residing in limited portions of the recording areas, contribute most of the beta oscillatory power identified by the bivariate time series.

Figure 14

Figure 14. Factor loadings of V4 and PFC, spatially smoothed, normalized, and color coded over the electrode coordinates (μm) on (left) V4 and (right) PFC at late delay period (400 ms). Contours at 0.25, 0.5 and 0.75 of the maximal power were added.

4 Discussion

To describe cross-population interactions of oscillatory amplitudes in high-dimensional neural field potentials, we developed a novel and intuitive procedure: in LaDynS, cross-population interactions are described using the cross-correlation of latent drivers. To arrive at a statistically rigorous approach, we provided a time-series extension of probabilistic CCA together with a novel sparse estimation methodology. According to our Equation 5, each of the two multivariate time series is driven by a single latent time series, with the cross-dependence of these two latent time series representing cross-region interaction. According to Equation 6, the latent bivariate time series is a discrete Gaussian process, but its correlation matrix is unrestricted, allowing for time-varying dynamics, i.e., statistical non-stationarity. The repeated trial structure enabled us to estimate the resulting high-dimensional covariance matrix by applying sparse estimation and inference methods. We found, and displayed in Figure 12b, an interesting interaction between PFC and V4, involving beta power, that appeared during the late delay period, with the interaction being bidirectional. The results were based on partial R² values, computed from the estimated covariance matrices, corresponding to lagged regressions of one latent time series on the other. The analysis in Figure 12b is in the spirit of Granger causality, but differs from it by allowing for non-stationarity, so that we could obtain the time-varying results.

The partial R² values in Figure 12b are highly statistically significant, but they are very small. We note that our diagonal regularization of the estimated latent precision matrix artificially reduced these values, so that the scale is no longer interpretable in familiar terms. The cross-correlation estimates are several orders of magnitude larger (see Supplementary Figure S12). However, they remain smaller than .1, as are the raw correlations at 400 and 425 ms, respectively, across individual electrodes in V4 and PFC (not shown). This reflects the highly inhomogeneous nature of the LFP voltage recordings together with the dominance of low frequencies (unfiltered LFP spectra typically follow 1/f^α trends, where f is frequency and α is roughly 1.5 to 2; see Kass et al. (2023))). The large number of repeated trials allowed us to find and document the results.

In addition to making the analysis possible, the repeated trial structure suggests substantive interpretation based on trial-to-trial variability. Although investigators take pains to make the experimental setting nearly the same on each trial, the inevitable small fluctuations in the way the subject interacts with the environment, together with changes in the subject's underlying state (involving alterations in motivational drive, for example), lead to observable fluctuations in behavior and in the recorded neural activity. Although the network sources of trial-to-trial variability in the PFC and V4 data are unknown, they produce the kind of correlated activity revealed in Figure 12b. To interpret it, we acknowledge there could be some potentially confounding trial-to-trial variation — such as changing motivational state, stimulus timing, bandpass-filter distortions, and volume conduction — that drive beta power in V4 and PFC, having just the right differential time lags to produce the correlated activity picked up by the partial R² values, as cautioned by Lintas et al. (2021)). Could such task-irrelevant pulses of activity change across time, within repetitions of the task, in such a way as to produce the peaks in Figure 12b? It is possible, but it would be surprising, especially when we consider contemporary ideas about beta oscillations during working memory tasks (Miller et al., 2018) along with the well-identified distinction between early and late visual processing that, presumably, corresponds to a distinction between feedforward and recurrent (bidirectional) flow of information (Supèr et al., 2001; Buschman and Miller, 2007; Del Cul et al., 2007; Yang et al., 2019; Chen et al., 2022; Mejias et al., 2016). The alternative we mentioned, that PFC and V4 are involved, together, in goal-directed visual processing and memory, with the two areas acting bidirectionally around 400 ms post-stimulus, seems likely.

There are many ways to extend the ideas developed here. While we applied LaDynS to LFP recordings, the framework is applicable more broadly to other slowly varying multidimensional time series. In particular, the proposed methodology can be extended to other neural modalities such as EEG, MEG, and fMRI. Although data from these modalities often violate global stationarity assumptions underlying classical approaches such as cross-correlograms (Guan et al., 2020; Klonowski, 2009; Chang and Glover, 2010), local stationarity remains reasonable when functional connectivity depends on cognitive/behavioral states that evolve on time scales longer than the sampling interval. While we have not yet explored such extensions empirically, LaDynS is particularly appealing in this regime, as it explicitly accommodates slowly varying connectivity and cross-correlation over time. At the same time, signals in EEG and MEG often reflect mixtures of activity from multiple brain structures and switches between metastable neural states, which may introduce additional confounding and complicate the interpretation of inferred associations.

LaDynS may perform less favorably in settings where within-region noise autocorrelations dominate the signal. The simulation results in Figure 4 of Bong et al. (2020)) illustrate the resulting loss of specificity under such adversarial conditions. In these cases, methods that explicitly model within-region autocorrelated noise, such as DLAG (Gokcen et al., 2022), can outperform LaDynS. In Bong et al. (2020)), we proposed an extension of LaDynS that addresses this issue by allowing the within-region noise processes ϵ_k to follow general time series structures and by modeling the latent processes driving each brain region as multidimensional. That brief report, which focused on data filtered at a different frequency band, did not provide the methodological details or inferential procedures developed here. An important direction for future study is to generalize the present framework to the setting considered in Bong et al. (2020)).

Another important direction for future research is a theoretical analysis of how regularization hyperparameters affect recovery of the underlying lag structure. In this study, we focused on proposing a regularization-based approach to functional connectivity estimation together with a data-driven hyperparameter calibration scheme and demonstrated its performance through simulation studies, while establishing theoretical guarantees for hyperparameter selection. Optimality remains an interesting topic for future investigation. In addition, developing a theory that accounts for dependence across experimental trials is of particular interest, as such dependence may arise from long-term memory or adaptation effects in neural recordings and could impact the validity of the proposed inference procedure.

Furthermore, for band-pass filtered data, such as those analyzed in Section 3.4, phase analysis (Klein et al., 2020) could be combined with amplitude analysis based on the complex normal distribution, as in Urban et al. (2023)) (which developed a latent variable model at a single time point). Multiple frequencies and non-linear association, such as cross-bicoherence, could be considered, along the lines of Tort et al. (2010)); Gallagher et al. (2017)); Aliramezani et al. (2024)); Abe et al. (2024)), as well. A different direction for additional research would be to simplify the version of LaDynS we have used here by imposing a suitable spatiotemporal structure on the latent time series. Regardless of whether such approaches are fruitful, the general framework of LaDynS could be of use whenever interest focuses on time-varying interactions among groups of repeatedly-observed multivariate neural time series.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Ethics statement

The animal study was approved by Institutional Animal Care and Use Committee of the University of Pittsburgh. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

HB: Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. VV: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Supervision, Validation, Writing – original draft, Writing – review & editing. EY: Conceptualization, Data curation, Investigation, Validation, Writing – original draft, Writing – review & editing. MS: Data curation, Investigation, Validation, Writing – original draft, Writing – review & editing. RK: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. Bong, Ventura and Kass are supported in part by NIMH grant R01 MH064537. Yttri is supported by NIH grant (1R21EY029441-01) and the Whitehall Foundation. Smith is supported by NIH (R01EY022928, R01MH118929, R01EB026953, P30EY008098) and NSF (NCS 1734901) grants.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2026.1703722/full#supplementary-material

References

Abe, T., Asai, Y., Lintas, A., and Villa, A. E. P. (2024). Detection of quadratic phase coupling by cross-bicoherence and spectral granger causality in bifrequencies interactions. Sci. Rep. 14:8521. doi: 10.1038/s41598-024-59004-8

PubMed Abstract | Crossref Full Text | Google Scholar

Abraham, I., Shahsavarani, S., Zimmerman, B., Husain, F. T., and Baryshnikov, Y. (2024). Hemodynamic cortical ripples through cyclicity analysis. Netw. Neurosci. 8, 1105–1128. doi: 10.1162/netn_a_00392

PubMed Abstract | Crossref Full Text | Google Scholar

Aliramezani, M., Farrokhi, A., Constantinidis, C., and Daliri, M. R. (2024). Delta-alpha/beta coupling as a signature of visual working memory in the prefrontal cortex. iScience 27:110453. doi: 10.1016/j.isci.2024.110453

PubMed Abstract | Crossref Full Text | Google Scholar

Anderson-Sprecher, R. (1994). Model comparisons and R². Am. Stat. 48, 113–117. doi: 10.2307/2684259

Crossref Full Text | Google Scholar

Bach, F. R., and Jordan, M. I. (2005). A probabilistic interpretation of canonical correlation analysis. Technical Report 688, Department of Statistics, University of California, Berkeley, Berkeley, CA.

Google Scholar

Belluscio, M. A., Mizuseki, K., Schmidt, R., Kempter, R., and Buzsáki, G. (2012). Cross-frequency phase-phase coupling between theta and gamma oscillations in the hippocampus. J. Neurosci., 32, 423–435. doi: 10.1523/JNEUROSCI.4122-11.2012

Crossref Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B. 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

Crossref Full Text | Google Scholar

Bong, H., Liu, Z., Ren, Z., Smith, M. A., Ventura, V., and Kass, R. E. (2020). Latent dynamic factor analysis of high-dimensional neural recordings. Adv. Neural Inf. Process. Syst. 33, 16446–16456.

PubMed Abstract | Google Scholar

Bruns, A., Eckhorn, R., Jokeit, H., and Ebner, A. (2000). Amplitude envelope correlation detects coupling among incoherent brain signals. Neuroreport 11, 1509–1514. doi: 10.1097/00001756-200005150-00029

PubMed Abstract | Crossref Full Text | Google Scholar

Buschman, T. J., and Miller, E. K. (2007). Top-down versus bottom-up control of attention in the prefrontal and posterior parietal cortices. Science 315, 1860–1862. doi: 10.1126/science.1138071

PubMed Abstract | Crossref Full Text | Google Scholar

Buzsáki, G., Anastassiou, C. A., and Koch, C. (2012). The origin of extracellular fields and currents–EEG, ECoG, LFP and spikes. Nat. Rev. Neurosci. 13, 407–420. doi: 10.1038/nrn3241

PubMed Abstract | Crossref Full Text | Google Scholar

Chang, C., and Glover, G. H. (2010). Time-frequency dynamics of resting-state brain connectivity measured with fMRI. NeuroImage 50, 81–98. doi: 10.1016/j.neuroimage.2009.12.011

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, Y., Douglas, H., Medina, B. J., Olarinre, M., Siegle, J. H., and Kass, R. E. (2022). Population burst propagation across interacting areas of the brain. J. Neurophysiol. 128, 1578–1592. doi: 10.1152/jn.00066.2022

PubMed Abstract | Crossref Full Text | Google Scholar

[Dataset] Snyder, A., Johnston, R., and Smith, M. (2022). Utah array recordings from visual cortical area V4 and prefrontal cortex with simultaneous EEG.

Google Scholar

Del Cul, A., Baillet, S., and Dehaene, S. (2007). Brain dynamics underlying the nonlinear threshold for access to consciousness. PLoS Biol. 5:e260. doi: 10.1371/journal.pbio.0050260

PubMed Abstract | Crossref Full Text | Google Scholar

Fisher, R. A. (1925). Statistical Methods for Research Workers. Number 3. London: Oliver and Boyd.

Google Scholar

Friedman, J., Hastie, T., and Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics 9, 432–441. doi: 10.1093/biostatistics/kxm045

PubMed Abstract | Crossref Full Text | Google Scholar

Fries, P., Reynolds, J. H., Rorie, A. E., and Desimone, R. (2001). Modulation of oscillatory neuronal synchronization by selective visual attention. Science 291, 1560–1563. doi: 10.1126/science.1055465

PubMed Abstract | Crossref Full Text | Google Scholar

Gabor, D. (1946). Theory of communication. Part 1: the analysis of information. J. Inst. Electr. Eng. 93, 429–441. doi: 10.1049/ji-3-2.1946.0074

Crossref Full Text | Google Scholar

Gallagher, N., Ulrich, K. R., Talbot, A., Dzirasa, K., Carin, L., and Carlson, D. E. (2017). Cross-spectral factor analysis. Adv. Neural Inf. Process. Syst. 30, 6842–6852.

Google Scholar

Gallego-Carracedo, C., Perich, M. G., Chowdhury, R. H., Miller, L. E., and Gallego, J. A. (2022). Local field potentials reflect cortical population dynamics in a region-specific and frequency-dependent manner. eLife 11:e73155. doi: 10.7554/eLife.73155

PubMed Abstract | Crossref Full Text | Google Scholar

Gokcen, E., Jasper, A. I., Semedo, J. D., Zandvakili, A., Kohn, A., Machens, C. K., et al. (2022). Disentangling the flow of signals between populations of neurons. Nat. Comput. Sci. 2, 512–525. doi: 10.1038/s43588-022-00282-5

PubMed Abstract | Crossref Full Text | Google Scholar

Goris, R. L. T., Movshon, J. A., and Simoncelli, E. P. (2014). Partitioning neuronal variability. Nat. Neurosci. 17, 858–865. doi: 10.1038/nn.3711

PubMed Abstract | Crossref Full Text | Google Scholar

Gregoriou, G. G., Gotts, S. J., Zhou, H., and Desimone, R. (2009). High-frequency, long-range coupling between prefrontal and visual cortex during attention. Science 324, 1207–1210. doi: 10.1126/science.1171402

PubMed Abstract | Crossref Full Text | Google Scholar

Guan, S., Jiang, R., Bian, H., Yuan, J., Xu, P., Meng, C., et al. (2020). The profiles of non-stationarity and non-linearity in the time series of resting-state brain networks. Front. Neurosci. 14:493. doi: 10.3389/fnins.2020.00493

PubMed Abstract | Crossref Full Text | Google Scholar

Hardoon, D. R., Szedmak, S., and Shawe-Taylor, J. (2004). Canonical correlation analysis: an overview with application to learning methods. Neural Comput. 16, 2639–2664. doi: 10.1162/0899766042321814

PubMed Abstract | Crossref Full Text | Google Scholar

Hotelling, H. (1992). “Relations between two sets of variates,” in Breakthroughs in Statistics, eds. S. Kotz and N. L. Johnson (New York, NY: Springer), 162–190. doi: 10.1007/978-1-4612-4380-9_14

Crossref Full Text | Google Scholar

Jankova, J., and Van De Geer, S. (2015). Confidence intervals for high-dimensional inverse covariance estimation. Electron. J. Stat. 9, 1205–1229. doi: 10.1214/15-EJS1031

Crossref Full Text | Google Scholar

Jun, J. J., Steinmetz, N. A., Siegle, J. H., Denman, D. J., Bauza, M., Barbarits, B., et al. (2017). Fully integrated silicon probes for high-density recording of neural activity. Nature 551, 232–236. doi: 10.1038/nature24636

PubMed Abstract | Crossref Full Text | Google Scholar

Kass, R. E., Bong, H., Olarinre, M., Xin, Q., and Urban, K. N. (2023). Identification of interacting neural populations: methods and statistical considerations. J. Neurophysiol. 130, 475–496. doi: 10.1152/jn.00131.2023

PubMed Abstract | Crossref Full Text | Google Scholar

Kettenring, J. R. (1971). Canonical analysis of several sets of variables. Biometrika 58, 433–451. doi: 10.1093/biomet/58.3.433

Crossref Full Text | Google Scholar

Khanna, S. B., Scott, J. A., and Smith, M. A. (2020). Dynamic shifts of visual and saccadic signals in prefrontal cortical regions 8Ar and FEF. J. Neurophysiol. 124, 1774–1791. doi: 10.1152/jn.00669.2019

PubMed Abstract | Crossref Full Text | Google Scholar

Klein, N., Orellana, J., Brincat, S. L., Miller, E. K., Kass, R. E., et al. (2020). Torus graphs for multivariate phase coupling analysis. Ann. Appl. Stat. 14, 635–660. doi: 10.1214/19-AOAS1300

PubMed Abstract | Crossref Full Text | Google Scholar

Klonowski, W. (2009). Everything you wanted to ask about eeg but were afraid to get the right answer. Nonlinear Biomed. Phys. 3:2. doi: 10.1186/1753-4631-3-2

PubMed Abstract | Crossref Full Text | Google Scholar

Lakshmanan, K. C., Sadtler, P. T., Tyler-Kabara, E. C., Batista, A. P., and Yu, B. M. (2015). Extracting low-dimensional latent structure from time series in the presence of delays. Neural Comput. 27, 1825–1856. doi: 10.1162/NECO_a_00759

PubMed Abstract | Crossref Full Text | Google Scholar

Liebe, S., Hoerzer, G. M., Logothetis, N. K., and Rainer, G. (2012). Theta coupling between V4 and prefrontal cortex predicts visual short-term memory performance. Nat. Neurosci. 15:456. doi: 10.1038/nn.3038

PubMed Abstract | Crossref Full Text | Google Scholar

Lintas, A., Sánchez-Campusano, R., Villa, A. E. P., Gruart, A., and Delgado-García, J. M. (2021). Operant conditioning deficits and modified local field potential activities in parvalbumin-deficient mice. Sci. Rep. 11:2970. doi: 10.1038/s41598-021-82519-3

PubMed Abstract | Crossref Full Text | Google Scholar

Mazumder, R., and Hastie, T. (2012). The graphical lasso: new insights and alternatives. Electron. J. Stat. 6:2125. doi: 10.1214/12-EJS740

PubMed Abstract | Crossref Full Text | Google Scholar

Mejias, J. F., Murray, J. D., Kennedy, H., and Wang, X.-J. (2016). Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex. Sci. Adv. 2:e1601335. doi: 10.1126/sciadv.1601335

PubMed Abstract | Crossref Full Text | Google Scholar

Miller, E. K., and Cohen, J. D. (2001). An integrative theory of prefrontal cortex function. Annu. Rev. Neurosci. 24, 167–202. doi: 10.1146/annurev.neuro.24.1.167

PubMed Abstract | Crossref Full Text | Google Scholar

Miller, E. K., Lundqvist, M., and Bastos, A. M. (2018). Working memory 2.0. Neuron 100, 463–475. doi: 10.1016/j.neuron.2018.09.023

PubMed Abstract | Crossref Full Text | Google Scholar

Mitra, A., Kraft, A., Wright, P., Acland, B., Snyder, A. Z., Rosenthal, Z., et al. (2018). Spontaneous infra-slow brain activity has unique spatiotemporal dynamics and laminar structure. Neuron 98, 297–305.e6. doi: 10.1016/j.neuron.2018.03.015

PubMed Abstract | Crossref Full Text | Google Scholar

Mitra, A., Snyder, A. Z., Blazey, T., and Raichle, M. E. (2015). Lag threads organize the brain's intrinsic activity. Proc. Nat. Acad. Sci. 112, E2235–E2244. doi: 10.1073/pnas.1503960112

PubMed Abstract | Crossref Full Text | Google Scholar

Ombao, H., and Pinto, M. (2022). Spectral dependence. Econometrics and Statistics.

Google Scholar

Orban, G. A. (2008). Higher order visual processing in macaque extrastriate cortex. Physiol. Rev. 88, 59–89. doi: 10.1152/physrev.00008.2007

PubMed Abstract | Crossref Full Text | Google Scholar

Perone Pacifico, M., Genovese, C., Verdinelli, I., and Wasserman, L. (2004). False discovery control for random fields. J. Am. Stat. Assoc. 99, 1002–1014. doi: 10.1198/0162145000001655

Crossref Full Text | Google Scholar

Pesaran, B., Vinck, M., Einevoll, G. T., Sirota, A., Fries, P., Siegel, M., et al. (2018). Investigating large-scale brain dynamics using field potential recordings: analysis and interpretation. Nat. Neurosci. 21, 903–919. doi: 10.1038/s41593-018-0171-8

PubMed Abstract | Crossref Full Text | Google Scholar

Raghavan, M., Pilet, J., Carlson, C., Anderson, C. T., Mueller, W., Lew, S., et al. (2024). Gamma amplitude-envelope correlations are strongly elevated within hyperexcitable networks in focal epilepsy. Sci. Rep. 14:17736. doi: 10.1038/s41598-024-67120-8

PubMed Abstract | Crossref Full Text | Google Scholar

Raut, R. V., Mitra, A., Marek, S., Ortega, M., Snyder, A. Z., Tanenbaum, A., et al. (2019). Organization of propagated intrinsic brain activity in individual humans. Cerebral Cortex, 30, 1716–1734. doi: 10.1093/cercor/bhz198

PubMed Abstract | Crossref Full Text | Google Scholar

Rodu, J., Klein, N., Brincat, S. L., Miller, E. K., and Kass, R. E. (2018). Detecting multivariate cross-correlation between brain regions. J. Neurophysiol. 120, 1962–1972. doi: 10.1152/jn.00869.2017

PubMed Abstract | Crossref Full Text | Google Scholar

Sarnthein, J., Petsche, H., Rappelsberger, P., Shaw, G., and Von Stein, A. (1998). Synchronization between prefrontal and posterior association cortex during human working memory. Proc. Nat. Acad. Sci. 95, 7092–7096. doi: 10.1073/pnas.95.12.7092

PubMed Abstract | Crossref Full Text | Google Scholar

Schmidt, R., Ruiz, M. H., Kilavik, B. E., Lundqvist, M., Starr, P. A., and Aron, A. R. (2019). Beta oscillations in working memory, executive control of movement and thought, and sensorimotor function. J. Neurosci. 39, 8231–8238. doi: 10.1523/JNEUROSCI.1163-19.2019

PubMed Abstract | Crossref Full Text | Google Scholar

Schmolesky, M. T., Wang, Y., Hanes, D. P., Thompson, K. G., Leutgeb, S., Schall, J. D., et al. (1998). Signal timing across the macaque visual system. J. Neurophysiol. 79, 3272–3278. doi: 10.1152/jn.1998.79.6.3272

PubMed Abstract | Crossref Full Text | Google Scholar

Schneider, M., Broggini, A. C., Dann, B., Tzanou, A., Uran, C., Sheshadri, S., et al. (2021). A mechanism for inter-areal coherence through communication based on connectivity and oscillatory power. Neuron 109, 4050–4067. doi: 10.1016/j.neuron.2021.09.037

PubMed Abstract | Crossref Full Text | Google Scholar

Semedo, J. D., Gokcen, E., Machens, C. K., Kohn, A., and Yu, B. M. (2020). Statistical methods for dissecting interactions between brain areas. Curr. Opin. Neurobiol. 65, 59–69. doi: 10.1016/j.conb.2020.09.009

PubMed Abstract | Crossref Full Text | Google Scholar

Semedo, J. D., Zandvakili, A., Machens, C. K., Yu, B. M., and Kohn, A. (2019). Cortical areas interact through a communication subspace. Neuron 102, 249–259. doi: 10.1016/j.neuron.2019.01.026

PubMed Abstract | Crossref Full Text | Google Scholar

Shahsavarani, S., Abraham, I. T., Zimmerman, B. J., Baryshnikov, Y. M., and Husain, F. T. (2020). Comparing cyclicity analysis with pre-established functional connectivity methods to identify individuals and subject groups using resting state fmri. Front. Comput. Neurosci. 13:94. doi: 10.3389/fncom.2019.00094

PubMed Abstract | Crossref Full Text | Google Scholar

Shao, J. (1993). Linear model selection by cross-validation. J. Am. Stat. Assoc. 88, 486–494. doi: 10.1080/01621459.1993.10476299

Crossref Full Text | Google Scholar

Steinmetz, N. A., Koch, C., Harris, K. D., and Carandini, M. (2018). Challenges and opportunities for large-scale electrophysiology with neuropixels probes. Curr. Opin. Neurobiol. 50, 92–100. doi: 10.1016/j.conb.2018.01.009

PubMed Abstract | Crossref Full Text | Google Scholar

Supèr, H., Spekreijse, H., and Lamme, V. A. (2001). Two distinct modes of sensory processing observed in monkey primary visual cortex (v1). Nat. Neurosci. 4, 304–310. doi: 10.1038/85170

PubMed Abstract | Crossref Full Text | Google Scholar

Tibshirani, R. J., and Taylor, J. (2012). Degrees of freedom in lasso problems. Ann. Stat. 40, 1198–1232. doi: 10.1214/12-AOS1003

Crossref Full Text | Google Scholar

Tort, A. B. L., Komorowski, R., Eichenbaum, H., and Kopell, N. (2010). Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. J. Neurophysiol. 104, 1195–1210. doi: 10.1152/jn.00106.2010

PubMed Abstract | Crossref Full Text | Google Scholar

Urban, K. N., Bong, H., Orellana, J., and Kass, R. E. (2023). Oscillating neural circuits: phase, amplitude, and the complex normal distribution. Can. J. Stat. 51, 824–851. doi: 10.1002/cjs.11790

PubMed Abstract | Crossref Full Text | Google Scholar

Ventura, V., Cai, C., and Kass, R. E. (2005). Statistical assessment of time-varying dependency between two neurons. J. Neurophysiol. 94, 2940–2947. doi: 10.1152/jn.00645.2004

PubMed Abstract | Crossref Full Text | Google Scholar

Wilks, S. S. (1932). Certain generalizations in the analysis of variance. Biometrika 24, 471–494. doi: 10.1093/biomet/24.3-4.471

Crossref Full Text | Google Scholar

Yang, Y., Tarr, M. J., Kass, R. E., and Aminoff, E. M. (2019). Exploring spatiotemporal neural dynamics of the human visual cortex. Hum. Brain Mapp. 40, 4213–4238. doi: 10.1002/hbm.24697

PubMed Abstract | Crossref Full Text | Google Scholar

Yu, B. M., Cunningham, J. P., Santhanam, G., Ryu, S. I., Shenoy, K. V., and Sahani, M. (2009). Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity. J. Neurophysiol. 102, 614–635. doi: 10.1152/jn.90941.2008

PubMed Abstract | Crossref Full Text | Google Scholar

Zamm, A., Debener, S., Bauer, A.-K. R., Bleichner, M. G., Demos, A. P., and Palmer, C. (2018). Amplitude envelope correlations measure synchronous cortical oscillations in performing musicians. Ann. N. Y. Acad. Sci. 1423, 251–263. doi: 10.1111/nyas.13738

PubMed Abstract | Crossref Full Text | Google Scholar

Zhu, Y., and Cribben, I. (2018). Sparse graphical models for functional connectivity networks: best methods and the autocorrelation issue. Brain Connect. 8, 139–165. doi: 10.1089/brain.2017.0511

PubMed Abstract | Crossref Full Text | Google Scholar

Zou, H., Hastie, T., and Tibshirani, R. (2007). On the “degrees of freedom” of the lasso. Ann. Stat. 35, 2173–2192. doi: 10.1214/009053607000000127

Crossref Full Text | Google Scholar