{'fraction_non_alphanumeric': 0.04577710251139512, 'fraction_numerical': 0.02407721869693449, 'mean_word_length': 4.331744972839036, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Latent Cross-population Dynamic Time-series Analysis of High-dimensional Neural Recordings

An important problem in analysis of neural data is to characterize interactions across brain regions from high-dimensional multiple-electrode recordings during a behavioral experiment. Lead-lag effects indicate possible directional flows of neural information, but they are often transient, appearing during short intervals of time. Such non-stationary interactions can be difficult to identify, but they can be found by taking advantage of the replication structure inherent to many neurophysiological experiments. To describe non-stationary interactions between replicated pairs of high-dimensional time series, we developed a method of estimating latent, non-stationary cross-correlation. Our approach begins with an extension of probabilistic CCA to the time series setting, which provides a model-based interpretation of multiset CCA. Because the covariance matrix describing non-stationary dependence is high-dimensional, we assume sparsity of cross-correlations within a range of possible interesting lead-lag effects. We show that the method can perform well in realistic settings and we apply it to 192 simultaneous local field potential (LFP) recordings from prefrontal cortex (PFC) and visual cortex (area V4) during a visual memory task. We find lead-lag relationships that are highly plausible, being consistent with related results in the literature.

1. Introduction. Contemporary technologies for recording neural activity can produce multiple time series in each of two or more brain regions (e.g., Jun et al., 2017;Steinmetz et al., 2018), enabling identification of interactions across regions that are relevant to behavior. In some situations, it may be possible to find lead-lag relationships, which could indicate directional flow of information. For stationary time series, this is the problem solved by Granger causality (Geweke, 1982). When recordings are made from an alert animal during a behavioral task, however, the neural activity is constantly evolving, and the statistical challenge is to find a method that allows for non-stationarity.
Various non-stationary methods have been applied to discovering simultaneous associations within brain regions (Buesing et al., 2014;Zhao and Park, 2017) or between regions (Gallagher et al., 2017), but their ability to find lead-lag relationships has not yet been established. The DKCCA method of  uses a dynammic kernel canonical correlation analysis to establish cross-region interaction, and it can estimate non-stationary lead-lag relationships, but it is unable to assess statistical significance of those relationships. In this paper we define and develop a framework that is able to detect lead-lag relationships across two brain regions, prefrontal cortex (PFC) and visual area V4, during a visual memory task, based on multiple time series of local field potential recordings (LFPs).
Visual cortical area V4 has been reported to retain higher order information (e.g. color and shape) and attention to visual objects (Orban, 2008;Fries et al., 2001), while prefrontal cortex (PFC) is considered to exert cognitive control in working memory (Miller and Cohen, 2001). Despite their spatial distance and functional differentiation, these regions have been presumed to cooperate during visual working memory tasks. For example, Sarnthein et al. (1998) discovered enhanced oscillatory activity in both regions during visual memory retention and Liebe et al. (2012) found delay-specific oscillatory coupling between V4 and lateral PFC. The nature and exact timing of this interaction, however, have remained unknown. One leading possibility is that PFC exerts control using oscillations in the beta range (16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) that are coupled across brain regions (Miller, Lundqvist and Bastos, 2018). The previous work of Klein et al. (2020) developed methods for analyzing multivariate phase coupling of such oscillations. Here, we turn our attention to coupling of the oscillatory power, while the subject was remembering the location of a visual target, with the goal being to identify lead-lag relationships in beta-band power between PFC and V4, observed in LFPs recorded simultaneously from the two regions.
LFPs result from neural currents generated near the electrode (roughly within 150-200 microns), involving large numbers of neurons (Buzsáki, Anastassiou and Koch, 2012;Einevoll et al., 2013;Pesaran et al., 2018), and they have been shown to correlate substantially with the BOLD fMRI signal (Logothetis et al., 2001;Magri et al., 2012). The recordings in this study come from a pair of 96-electrode arrays implanted in PFC and V4. The data were sampled every millisecond across a 500 millisecond time interval, with 3000 replications (3000 experimental trials). Each of the 3000 observations forms a 192-dimensional time series, in two groups of 96-dimensional time series. The method we developed, Latent Dynamic analysis via Sparse banded graphs (LaDynS), was able to identify non-stationary cross-area interactions from the PFC and V4 data. Our data-analytic findings are summarized in Figs. 6-9, and presented in full in Section 4.
According to the LaDynS model, each of the two multivariate time series is driven by a latent univariate time series, with the resulting latent bivariate time series potentially representing the time-varying cross-area dependence we seek to identify. Instead of assuming a parametric structure for the cross-dependence of the latent bivariate time series, however, we have left the correlation matrix unspecified. As we show in Section 2.2, assuming each multivariate observation, at each time point, depends linearly on a latent normal random variable, having mean 0 and variance 1, we then get a time series generalization of probabilistic CCA (pCCA; Bach and Jordan (2005)), as depicted graphically in Figure 1. Furthermore, maximum likelihood applied to the model gives the same data summary as multiset CCA ). This provides a nice interpretation, but it isn't yet practical because the latent time series correlation matrix, which must be estimated, has 2T (2T − 1)/2 parameters, where T is the number of time points. We therefore resorted to sparse estimation methods.
As we describe in Section 2, LaDynS uses L1 regularization on the cross-precision parameters within a band defined by a maximal lead or lag of the partial cross-correlations. There is also a maximal lead or lag of the partial auto-correlations, but the remaining partial auto-correlations (at less than the maximum lead or lag) are not regularized. Because the data we analyzed were band-pass filtered, the time series were unusually smooth. This led us to consider and apply additional regularization along the diagonals of the auto-correlation matrices. Statistical inference in this context can be based on false discovery rate (FDR) based on p-values computed from asymptotics for a desparsified version of the estimated precision matrix (Jankova et al., 2015). Because we are interested in epochs of correlated activity, we used a modification of the STAR method of finding clustered contiguous significant statistics (a) (b) Fig 1: (a) Graphical representation of the pCCA model of Bach and Jordan (2005), where X 1 and X 2 are random vectors and Z is a random variable. (b) A variation on (a) that facilitates extension to the case in which X 1 and X 2 are multivariate time series and (Z 1 , Z 2 ) is a bivariate time series. . Our modification allowed the method to find multiple, nonconvex regions. Simulations in Section 3 show that LaDynS is able to correctly identify the timing, relative to behavior, of interactions between regions, as illustrated in Fig. 2b, when applied to artificial data designed to be similar to those we analyzed. The simulations make credible our data-analytic results in Section 4. We add some discussion in Section 5.
2. Methods. We begin by reviewing and reformulating probabilistic CCA (pCCA) in Section 2.1, then generalize to time series, in Section 2.2, obtaining to form a dynamic version of pCCA. Theorem 2.3 establishes an equivalence between the GENVAR version of multi-set CCA ) and maximum likelihood applied to our dynamic pCCA. We define LaDynS based on the loglikelihood function in Eq. (16) of Section 2.3. We go over choice of regularization parameters in Section 2.3.1, application to smooth time series in Section 2.3.2, and our algorithm for solving the penalized maximum likelihood problem in Section 2.3.3. We discuss statistical inference in Section 2.4.
2.1. Probabilistic CCA for two random vectors. Given two random vectors X 1 ∈ R d1 and X 2 ∈ R d2 , canonical correlation analysis (CCA) (Hotelling, 1992) finds the sets of weights w 1 ∈ R d1 and w 2 ∈ R d2 that maximize Pearson's correlation between linear combinations w 1 X 1 and w 2 X 2 . This can be rewritten as (1) σ cc = max wk,k=1,2:w k Σkkwk=1 w 1 Σ 12 w 2 where Σ kk = Var(X k ) is the covariance matrix of X k , k = 1, 2, and Σ 12 = Cov(X 1 , X 2 ) the cross-covariance matrix between X 1 and X 2 . The sample estimator σ cc is obtained by replacing Σ kk and Σ 12 with their sample analogs Σ kk and Σ 12 respectively. The maximizing weights w k and linear combinations w k X k are referred to as the canonical weights and canonical variables, respectively. Probabilistic CCA assumes that X 1 and X 2 are driven by a common one dimensional latent variable Z: where µ k ∈ R dk and β k ∈ R dk are mean vectors and factor loadings, respectively, and k indep ∼ MVN(0, Φ k ) (Bach and Jordan, 2005). Fig. 1a depicts the dependence of X 1 and X 2 on Z. The parameters in Eqs. (1) and (2) relate as follows: THEOREM 2.1 (Bach and Jordan, 2005, Theorem 2). The maximum likelihood estimators (MLEs) β 1 and β 2 in Eq. (2) based on N observed vector pairs X 1, [n] , X 2,[n] n=1,2,...,N are equivalent to the CCA solution w 1 , w 2 , and σ cc in Eq. (1) according to: (3) β k = Σ kk w k m k , where m 1 m 2 = σ cc and |m k | ≤ 1, k = 1, 2.
Theorem 2.1 proves that the original CCA setting and the generative pCCA model both yield the same estimate of σ cc . The advantage of the model-based approach is that maximum likelihood estimation allows inferences to be made about σ cc or other parameters that may be of interest. It also provides a framework for generalization: we introduce a pCCA extension that allows distinct latent variables for X 1 and X 2 , as depicted in Fig. 1b. Specifically, we assume that where µ k ∈ R dk , β k ∈ R dk and Φ ∈ R dk×dk are defined as in Eq.
(2), and (Z 1 , Z 2 ) are bivariate normally distributed: In this model, the canonical correlation σ cc is equal to the MLE of σ 12 , as stated in Theorem 2.2 below, so inference about σ cc is easier than inference based on model Eq.
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 (t) 1 ∈ R d1 and X (t) 2 ∈ R d2 , t = 1, 2, . . . , T . We use Eq. (4) to model the dependence of X (t) k on its associated latent variable Z (t) k at each time t: k are defined as in Eq. (4), and , and we impose the same constraints as in Eq. (7) to ensure parameters identifiability: Then for each t, we could define a parameter σ (t) 12 as in Eq. (5), to capture population-level association between X (t) 1 and X (t) 2 at t. But because we are also interested in lagged associations between X where Σ captures jointly all simultaneous and lagged associations within and between the two time series. Fig. 2a  is the partial-correlation between the two time-series at times t and s, respectively. Thus, when an element of Ω 12 is non-null, depicted as the red star in the expanded display within Fig. 2b, 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 the connection within the behavior and the direction of information flow between two brain regions, respectively, which are useful in identifying the functional nature of the connection.
Theorem 2.2 provided an equivalence between a non-distributional method (CCA) and its probabilistic representation (pCCA). We now derive a similar connection between the multiset generalization of CCA introduced by  and the dynamic pCCA model Eqs. (9), (10), (11) and (12). Multi-set CCA applied to 2T random vectors {X 2 : t = 1, . . . , T }. In particular, the GENVAR extension minimizes the generalized variance of these linear combinations, defined as the determinant of their correlation matrix (Wilks, 1932): where Var denotes the sample variance-covariance matrix and the weights w (t) k are scaled so that every diagonal entry of the matrix is 1. 2 , and estimated by their cross-precision matrix Ω 12 . (b) When a significant cross-precision entry is identified, e.g., the red star in the expanded view of Ω 12 , 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 Ω 12 , which means that, at this particular time, region 1 leads region 2, or Z 1 → Z 2 in short (a non-zero entry in the lower diagonal would mean Z 2 → Z 1 ). We represent this association by the red arrow on the right-most plot, with a lag of two units of time for illustration.

Latent Dynamic Analysis via Sparse Banded Graphs (LaDynS).
Our goal is to estimate the association dynamics between two multivariate time series using the covariance matrix Σ of their associated latent time series in Eq. (14). 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 Ω = Σ −1 in Eq. (14), rewriting log det(Σ) = log det(Ω −1 ) = − log det(Ω), and assuming that Ω has the banded structure depicted in Fig. 3.  k , t = 1, 2, . . . , T, k = 1, 2 and precision matrix Ω that minimize the penalized negative log-likelihood: 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 Λ (t,s) 12 = ∞ when |t − s| > d cross to force the corresponding cross-precision elements to zero and thus impose a banded structure on Ω 12 . We apply sparsity constraint Λ (t,s) 12 = λ cross > 0 on the remaining off-diagonals of Ω 12 to focus our discovery on sparse dominant associations and reduce the effective parameter size. We proceed similarly with the auto-precision matrices Ω 11 and Ω 22 , 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 (t) k,i , k = 1, 2, i = 1, . . . , N , and impose no further sparsity (λ auto = 0) unless there is reason to expect it.
Notice that to facilitate the choice of Λ, we grouped its elements into diagonal and offdiagonal elements and assigned the same penalties, λ cross , λ auto = 0 and λ diag , within each group.
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 Ω 12 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 Ω 12 , we choose instead Fig 3: The elements of Ω kk , k = 1, 2, and Ω 12 are set to zero outside of the gray bands of widths (1 + 2d auto ) and (1 + 2d cross ), respectively. 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. We study the properties of this method in Section 3.

2.3.2.
LaDynS on smooth time-series data. When Σ is not invertible, as is the case for the band-pass filtered experimental data we analyze in Section 4, we penalize its diagonal by λ diag > 0. We pick the smallest λ diag such that when we add λ diag I T to each block of Σ it is invertible and the resulting precision matrix starts to describe a proper auto-correlation structure in the following sense. Let X . . , N } be the observed time series with the largest significant auto-correlation, and S k * ,i * ∈ R T ×T be its auto-correlation matrix. Because X (t) k * ,i * is auto-correlated, we should see a band of negative values around the main diagonal of S −1 k * ,i * , but that band is not visible if S −1 k * ,i * is degenerate. We take λ diag to be the smallest value such that (S k * ,i * + λ diag I T ) −1 appears banded.
2.3.3. Fitting LaDynS. Eq. (16) 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 following coordinate descent algorithm finds a minimum, possibly local.
Assuming that all canonical weights w (t) k are fixed, Eq. (16) reduces to the gLASSO problem: which we can solve efficiently using a number of existing algorithms; here we use the P-gLASSO algorithm of . Then assuming that all parameters are fixed but a single weight w (t) k , Eq. (16) can be re-arranged as the linear problem: for which an analytical solution is available. That is, our algorithm alternates between updating Ω and the weights w

Inference for associations between two vector time series.
Let Ω and w (t) k , t = 1, . . . , T , k = 1, 2, be the LaDynS estimates of canonical precision matrix and canonical weights, and Σ = Var w be the empirical covariance of the estimated latent variables, defined in Eq. (16). Note that Ω = Σ −1 since Ω is constrained to be sparse. Based on these estimates, we want to identify the non-zero partial cross-correlations in Ω 12 , that is to identify the epochs of association between the two time series.
Formal inference methods for Ω based on its LaDynS estimate (Eq. (16)) are not available, but because LaDynS reduces to graphical LASSO (gLASSO) when the weights w (t) k in Eq. (17) are fixed, we co-opt gLASSO inference methods. Specifically, Jankova et al. (2015) suggested de-sparsifying the gLASSO estimate Ω according to and proved that, under mild assumptions and as n → ∞, each entry of Ω satisfies the Central Limit Theorem with center the true precision Ω: 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 (Bong et al., 2021a, Section S5). Jankova et al. (2015) also proposed an estimator of σ (t,s) 12 , but it is likely to be downward biased in our framework since estimating the canonical weights w to test H (t,s) 0 : Ω (t,s) 12 = 0, for each (t, s) ∈ [T ] 2 within d cross of the diagonal of Ω 12 . If we are specifically interested in excitatory or inhibitory connections, we use one-sided tests.
Control of false discoveries Because we perform tests for many entries of Ω 12 , we control the number of spurious detections by applying a procedure to guarantee the false discovery rate , where FDP = #{falsely discovered entries} #{discovered entries} ∨ 1 is below a nominal value. We focus on the STAR-SeqStep method of Lei, Ramdas and Fithian (2017) (STAR for short), designed to identify contiguous epochs of correlated activity rather than isolated effects, which we expect for the data we analyze here.  provided an algorithm to detect a single convex area of significant effects that we extended to detect multiple such areas; see Section S3 of the Supplementary Material (Bong et al., 2021a).
As a last step and further safeguard against falsely detecting correlated activity between brain areas, we obtain p-values for each connectivity epoch identified by STAR using the excursion test of Xu et al. (2011) and . See Section S4 of the Supplementary Material for details (Bong et al., 2021a). ..,N . The resulting sample contains no correlated activity. Hence, applying LaDynS to the bootstrap sample yields estimates of canonical precision matrix Ω * , canonical weights w * k (t)s, empirical covariance of the estimated latent variables , and de-sparsified precision matrix estimate Ω * (Eq. (19)) under the global null hypothesis of no correlated activity. Repeating the bootstrap simulation B times produces B bootstrap values is obtained under the global null hypothesis -i.e. under H (t,s) 0 : Ω (t,s) 12 = 0 simultaneously for all (t, s) -because it is not trivial to simulate bootstrap data that satisfy a specific H (t,s) 0 without assuming that all other elements of Ω 12 are also null. We garnered from simulations that σ (t,s) 12 is thus likely to slightly underestimate σ (t,s) 12 , which makes for slightly sensitive p-values.
Directional lagged influence: For each time point t, and lags ranging from τ 1 to τ 2 prior to t, (and vice versa) can be obtained from the estimated covariance matrix Σ and precision matrix Ω using standard formulas. Let us assume Z (t) 1 and Z (t−τ1:t−τ2) 2 are centered and write the covariance matrix in block form with the precision matrix being in corresponding block form. Then, writing It is possible to get spurious large R 2 values due to autocorrelation. To test significance we compare the estimated R 2 2→1 (t) and R 2 1→2 (t) with their null distributions under independence, based on permutation tests obtained by permuting trials in one region. The trialpermuted values of R 2 are based on null time series that exhibit the autocorrelation structure of the original latent time series. These R 2 values are, therefore, in the spirit of Granger causality, except that we do not need to assume stationarity: we can avoid it because we have repeated trials.
3. LaDynS performance on simulated data. We have introduced LaDynS to estimate the dynamic connectivity between two or more multivariate time series, and proposed inferential procedures to identify when connectivity is statistically significant. We apply LaDynS to experimental data in Section 4, but first, we examine its performance using simulated datasets that have properties similar to the experimental data. Code for simulations and experimental data analyses is provided at Bong et al. (2021b).
3.1. Simulated datasets. One simulated dataset consists of N = 3000 i.i.d. vector timeseries X 1 and X 2 of dimensions d 1 = d 2 = 96 and durations T = 100, simulated from Eq. (9). The latent time series Z 1 and Z 2 in Eqs. (9) and (11) have zero mean vectors and covariance matrix Σ = Ω −1 , with where Ω 12 is the cross-precision matrix of interest. The elements of the auto-precision matrices of Z k , for k = 1, 2, were simulated from the squared exponential function: (27) Σ (t,s) 0,k = exp −c 0,k (t − s) 2 , with c 0,1 = 0.148 and c 0,2 = 0.163 chosen to match the observed LFPs autocorrelations in the experimental dataset. The diagonal regularizer λI T was added to ensure that Σ 0,1 and Σ 0,2 are invertible, and we set λ = 0.5 to match the regularization λ diag we used in the experimental data analysis in Section 4. For Ω 12 , we considered the connectivity scenario depicted in Fig. 4d, where the two latent times series connect in three epochs, the first with no latency, the second with series 2 preceding series 1, and the third with series 1 preceding series 2. We accordingly set the cross-precision matrix elements to where r measures the intensity of the connection. To force the simulated data to be as smooth as the experimental data analyzed in Section 4, we then subtracted cI T from Σ = Ω −1 , with c chosen so that the smallest eigenvalue of (Σ − cI T ) equals the smallest eigenvalue of the correlation matrix of the experimental data. Finally, we rescaled Σ to have diagonal elements equal to one. Once the latent time series Z 1 and Z 2 were generated, we simulated a pair of observed time series according to for k = 1, 2 and t = 1, . . . , T , where Y  2 to be the two multivariate time-series neural recordings analyzed in Section 4, which we permute to remove all cross-correlations, and β k ∈ R dk×T to be the factor loadings estimated in Section 4. The resulting latent time series β (t) k Z (t) k and noise baseline vector Y (t) k in Eq. (29) have comparable scales and autocorrelations by construction, for k = 1, 2, to ensure that the simulated observed time series k have scales and auto-correlations comparable to the experimental data in Section 4. Moreover, the simulation parameters satisfy the identifiability conditions in Eq. (10), due to the second term on the right hand side of Eq. (29).
3.2. LaDynS estimation details. The simulated time series X (t) k were very smooth, so we needed to add the the regularizer λ diag to the diagonal entries of the estimated correlation matrix Σ so it could be inverted. This introduced inevitable bias to both sparsified and desparsified LaDynS' precision estimates, Ω and Ω. The other hyperparameters were set to d auto = 10, d cross = 10, and λ auto = 0. The penalty on the cross-correlation elements, λ cross , was automatically tuned at every run to control false discoveries (see Section 2.3).
3.3. Results. Figure 4a  cluster discovery rates, (c) false non-discovery rate and (d) false cluster non-discovery rate of STAR followed by excursion test at significance level α = 0.05 to identify non-zero partial correlations, for target FDR ∈ [0, 10]%, under the connectivity scenario in Fig. 4d, for connectivity intensities r = 0.002, 0.004, 0.006, 0.008 and 0.010. The dotted line is a (0,1) line. STAR does not control FDR, due to the bias introduced by λ diag -regularization. However, FCDR is controlled for all r and nominal FDR. Additionally, FNR and FCNR decrease to zero as the nominal FDR increases, suggesting that STAR has high power to detect all connectivity epochs, unless r is very small.
analyzed in Section 4. Figure 4b shows the bootstrap p-values for the entries of the desparsified cross-precision estimate Ω 12 (Eq. (21) with bootstrap simulation size B = 200; see Section 2.4). Small p-values concentrate near the locations of true non-zero cross-precision entries and are otherwise scattered randomly. We applied the STAR procedure with target FDR 5% (Section 2.4) and plotted the discovered epochs of connectivity in Fig. 4c: they match the true epochs well. Notice that, by design, STAR discovered only contiguous epochs and no scattered small effects. The three connectivity epochs identified by STAR were further found to be statistically significant using the excursion test (Bong et al., 2021a, Section S4), with p-values less than 0.005. We calculated false cluster discovery and non-discovery rates Perone Pacifico et al. (2004), estimating them by the corresponding proportions over 60 simulated datasets, with true connectivity in Fig. 4d, for a range of connectivity intensities r (Eq. (28)) that span from half to two times the largest estimated partial correlation value in the LFP experimental data of Section 4. We define a cluster to be falsely discovered if it contains no true effect, and falsely non-discovered if no estimated cluster overlaps with it. We further applied the excursion test to each discovered clusters and eliminated clusters with p-values larger than α = 0.05.
In Fig. 5, the discovery proportions are plotted against nominal FDR values between 0 and 10%. Although STAR does not control FDR in this simulation setting (panel (a)), it does control the proportion of false positive clusters, as the estimated FCDR is capped by the nominal FDR (panel (b)). FNR and FCNR (panels (c) and (d)) decrease quickly towards zero as the nominal FDR increases, suggesting that, unless the true connectivity is very weak (small r), LaDynS has high power to detect all connectivity epochs.
The precursor to LaDynS is DKCCA , with both methods aiming to identify epochs of lead-lag relationships between brain regions in non-stationary situations. The two methods are compared for simulated data with ground truth in Fig. 4d. In the Supplementary Material (Bong et al., 2021a), we present in Figure S6.3 the simulated data true cross-correlation matrix, obtained by inverting the true precision matrix, and the estimated cross-correlation matrices using by LaDynS and DKCCA. Both methods recover the essential structure of the true cross-correlations, although LaDynS does a better quantitative job. Importantly, unlike DKCCA, LaDynS also provides detailed assessment of statistical significance.
One possible source of problems with FDR control for the data simulated as in Section 3.1 could be sensitivity of the bootstrap p-values in Eq. (21) (see the comment at the end of Section 2.4), but the results are more likely due the use of the penalty on diagonal elements of Σ to guarantee invertibility, which induced bias in the estimate of Ω (see Section 3.2). That is, because the data were smooth, additional regularization was used in this simulation, and some zero entries in Ω were likely estimated to be significantly different from zero. Regularization also induces correlations between p-values, which STAR is known to be sensitive to . Section S6.2 of the Supplementary Material (Bong et al., 2021a) contains results based on simulated data that are not as smooth and therefore required no λ diag -regularization to run LaDynS. In this case, we found that STAR almost perfectly controls the FDR. We presented the more difficult case of very smooth data, here, because the experimental data we analyse in Section 2.4 is similarly smooth.

Experimental Data Analysis.
We applied LaDynS to LFPs recorded from PFC and V4 during a memory-guided saccade task, using two Utah arrays implanted in a macaque monkey, as described in Khanna, Scott and Smith (2020). Each trial of the task started with a monkey fixating its eyes toward the center of the screen. A visual cue was given for 50ms to indicate a target location, which was randomly chosen from eight locations. The monkey had to remember the target location while maintaining eye fixation for a delay period of 500ms. After the delay period, the monkey reoriented its eyes toward the remembered position, and reward was given on successful trials. As in Khanna, Scott and Smith (2020), we analyzed the time series during the delay period, based on 3000 successful trials. Because beta oscillations are often associated with communication across brain areas (Klein et al., 2020;Miller, Lundqvist and Bastos, 2018), we first filtered LFP recordings using complex Morlet wavelets at beta oscillation frequency 18Hz (which was the most active frequency in the magnitude scalograms), and obtained the beta oscillation power envelopes as the absolute values of the filtered signals. We note that the target frequencies were different from those in Bong et al. (2020), yielding different data and somewhat different results. After downsampling them to 200Hz, we applied LaDynS. Because the filtered data were very smooth, we used a regularization λ diag > 0 on the diagonal of Σ, as in Section 2.3. Figure 6 displays the LaDynS cross precision estimate Ω 12 , the bootstrap p-values for the entries of the desparsified cross-precision estimate Ω 12 , and the epochs of connectivity identified by STAR at target FDR 5%. Time t = 0 corresponds to when the visual stimulus (c) Discovered regions of cross-precision using STAR at nominal FDR 5%, suggesting that activity in PFC preceded that in V4 immediately post-stimulus and around 0.1 seconds thereafter (blue blobs), and that activity in V4 preceded that in PFC around time points 0.2 and 0.4 seconds post-stimulus (red blobs). These four epochs had excursion test p-values of p = .024, p = .0015, p = .001 and p < .0005, respectively. The two smaller blobs were not significant (p > 0.5) and are therefore likely spurious.
(the memory cue) was turned off. Three highly significant contiguous regions of the precision matrix (p ≈ .001 or smaller), discovered by STAR, are shown in Fig. 6. These correspond to beta power in V4 being correlated with earlier beta power in PFC, after conditioning on beta power at all other times and lags, about 100 millisecond after the cue was extinguished (blue blobs), and beta power in PFC being correlated with earlier beta power in V4 (after conditioning on beta power at all other times and lags), about 200 and 400 milliseconds post-stimulus (red blobs). It is important to keep in mind that the sparse estimation method finds a relatively small number of large effects; smaller effects, some of which may remain interesting, do not appear.
To better understand the lead-lag relationships, we used the estimated covariance matrix Σ, which is the empirical covariance matrix of the latent time series (see Section 2.4), to compute a series of R 2 values from interpretable regressions, as described in Section 2.4. We used τ 1 = .03 and τ 2 = .06 seconds. The null distribution was based on 3000 trial permutations. Fig. 7 displays two sets of R 2 values, the first set, in panel (a), corresponds to regressions of Y (t) on X (t−.03:t−.06) , where Y was the latent time series for V4 and X the latent time series for PFC. Thus, PFC leads V4 by about 45 milliseconds in these regressions. Panel (b) displays R 2 values with Y being the latent time series for PFC and X the latent time series for V4.
There are two large excursions of R 2 above the null values in the PFC → V4 plot and one in the V4 → PFC plot. The first excursion in the PFC → V4 plot, beginning around 100 ms, and the single excursion in the V4 → PFC plot centered around 400 ms, are consistent with the results in Fig. 6. The second excursion in the PFC → V4 plot, combined with the excursion in the V4 → PFC plot, shows mutually predictive activity.
These results are consistent with the idea that PFC exerts top-down control in a taskdependent manner, and V4 responds strongly to visual stimuli. PFC modulates V4, and PFC receives input from the visual processing stream, which includes V4, but we would expect the directional influences to differ across time, during the task, with PFC influencing V4 early in the task. Fig. 8 displays the estimated total beta power due to the latent time series, across time, for both PFC and V4. The latent beta power in PFC is largest at the beginning of the task, then drops to a roughly constant amplitude, while the latent beta power in V4 initially drops slightly but then starts to rise just before 100 ms, peaks around 250 ms, and then declines but remains large. Taking Fig. 7 and Fig. 8 together, these results show that taskdependent beta power in PFC leads beta power in V4, as identified using linear regression with time lags of 30 to 60 ms, and this becomes visible in Fig. 7 once beta power in V4 starts to increase, at around 100 ms. Later in the task, starting a little after 300 ms, each of the two latent beta power envelopes is able to predict the other, 30 to 60 ms in the future, suggesting the two areas are participating in bidirectionally linked beta oscillations. It is also possible to get information from the normalized factor loadings, across the electrode arrays, which are displayed in Fig. 9 for times t = 150 and t = 400 ms. (The loadings have been divided by their maximal value across the array.) Most noticeably, the spatial pattern of activity in PFC at 400 ms, compared with that at 150 ms, becomes more concentrated in the upper left corner of the electrode array. Bootstrap confidence intervals suggest this change in spatial pattern within PFC is not due solely to noise in estimation, see Figure S7.1 in the Supplementary Material (Bong et al., 2021a). Thus, apparently, the early and late periods of the task may be distinguished not only by very different prediction effects but, in addition, by spatial shifts in lag-correlated activity. (An animation over the complete timeline is available at Bong et al., 2021b.) 5. Discussion. Motivated by the problem of describing interactions across a pair of brain regions, based on high-dimensional neural recordings, we have provided a time-series extension of probabilistic CCA together with a novel sparse estimation methodology. According to our Equation (9), 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 (11), the latent bivariate time series is a discrete Gaussian process but its correlation matrix is unrestricted, allowing for non-stationarity. The repeated trial structure enabled us to estimate the resulting high-dimensional covariance matrix by applying sparse estimation and inference methods. Substantively, as displayed in Fig. 7, we found interesting directional relationships between beta-band LFPs in PFC and V4. These appeared during two distinct epochs, one early and one late in the task. In the earlier epoch PFC was leading V4, while in the later one the relationship was bidirectional. These were based on R 2 values, computed from the estimated covariance matrices, corresponding to lagged regressions of one latent time series on the other. The analysis in Fig. 7 is in the spirit of  9)) is plotted across time t for V4 (blue) and PFC (orange).  Fig. 7a (0.15 sec) and (right) the second bump of Fig. 7a or the first bump of Fig. 7b (0.40 sec). Contours at fractions .25, .5, and .75 of maximal power have been added.
Granger causality, but differs from it by allowing for non-stationarity, so that we could obtain the time-varying results.
In addition to making the analysis possible, the repeated trial structure suggests substantive interpretation based on trial-to-trial variability. Neurophysiological experiments like the one generating the data analyzed here take pains to make the experimental setting nearly the same on each trial. The inevitable, relatively small fluctuations in the way the subject interacts with the environment, together with changes in the subject's underlying state (involving fluctuations 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 Fig. 7. To interpret it, we acknowledge there could be some task-irrelevant, trial-varying pulses of inputs that drive beta power in V4 and PFC, having just the right differential time lags to produce the correlated activity picked up by the R 2 values we computed in the first peak of Fig. 7a, which would render those results merely coincidental. Could such task-irrelevant pulses of activity change across time, within repetitions of the task, in such a way as to produce, in addition, the second peak in Fig. 7a and also the peak in Fig. 7b? It is possible, but it would be surprising, especially when we consider contemporary ideas about beta oscillations during working memory tasks (Miller, Lundqvist and Bastos, 2018) along with the well-identified distinction between early and late visual processing, e.g., Yang et al. (2019). The alternative we mentioned, that PFC and V4 are involved, together, in goaldirected visual processing and memory, with PFC influencing V4 early and the two areas acting bidirectionally later, seems a credible possibility.
Figs. 5c and 5d shows LaDynS has good power to detect epochs of connectivity under realistic conditions. The identifiability constraint in Eq. (10), however, can lower statistical power when applied to data that do not meet that constraint. This is a feature shared by many methods that require constraints on parameters, but it should be kept in mind. As an illustration, in Section S6.1 of the Supplementary Material (Bong et al., 2021a), we simulated a dataset as described in Section 3.1, except that we omitted the second term β in the right hand side of Eq. (29), so that the resulting data did not satisfy the identifiability constraint. Then the effect of imposing the constraint to run LaDynS was to bias to various degrees the elements of Ω 12 , as seen in Figure S6.1 of Bong et al. (2021a). We cannot quantify this bias, but extensive simulations (not shown) suggest that effects are "smeared" (dampened and spread), and weak, typically statistically insignificant ripples can appear around the effects. In this particular simulation, the procedure was unable to detect the second cluster entirely, and the third cluster partially. Of course, as with all statistical procedures, larger samples and stronger effects yield more power. Also, in practice, data may or may not meet the constraint, so statistical power is not necessarily lost.
There are many ways to extend the ideas developed here. For band-pass filtered data, such as those analyzed in Section 4, phase analysis (Klein et al., 2020) could be combined with amplitude analysis. Multiple frequencies could be considered (along the lines of Gallagher et al. (2017)), as well. In Bong et al. (2020) we have described an extension of LaDynS, where the within-region noise vectors k were allowed to have general time series structure and the latent time series driving each brain region were multidimensional. That brief report, however, does not provide the details given here, nor does it include inferential procedures. Thus, an important future step will be to show how inference can be carried out for the general models in Bong et al. (2020). A different direction for additional research would be to simplify the version of LaDynS we have used here by imposing suitable spatiotemporal structure on the latent time series. While we hope such approaches will be fruitful, we believe the general framework of LaDynS can be of use whenever interest focuses on non-stationary interactions among groups of repeatedly-observed multivariate time series.

SUPPLEMENTARY MATERIAL
Supplement A: Supplement to "Latent Cross-population Dynamic Time-series Analysis of High-Dimensional Neural Recordings" (doi: COMPLETED BY THE TYPESETTER; .pdf). A PDF manuscript providing with supplemental descriptions of proofs, algorithms, and extensive details regarding arguments conveyed in the main text.
Supplement B: LaDynS Code Package with Reproducible Examples and Data to "Latent Cross-population Dynamic Time-series Analysis of High-Dimensional Neural Recordings" (doi: COMPLETED BY THE TYPESETTER; .zip). LaDynS Python package with tutorials based on the reproducible scripts generating the results given in the main text. The latest version will be updated at github.com/HeejongBong/ladyns. S1. Proof of Bong et al. (2021) 2 ) and δ (t,s) kl = 1, k = l and t = s, 0, otherwise, for t, s ∈ [T ] and k, l = 1, 2. Then Eq. (10) in Bong et al. (2021) can be rewritten into kl . Due to the orthogonality between u (t) k and Ψ (t) k , the calculation of det(R) and R −1 is straightforward: det(R) = det(Ω)/ k,t pdet(Ψ (t) k ) and Q = R −1 consists of sub-matrices where Ω = Σ −1 is the precision matrix and pdet(A) and A + are the pseudo-determinant and Moore-Penrose pseudo-inverse of a positive semi-definite matrix A. Notice that Ψ and hence pdet(Ψ (t) k ) = 1. In turn, the negative log-likelihood under the model (Bong et al., 2021, Eqs. (9), (10) and (11) k , X for all t ∈ [T ] and k = 1, 2. In this case, . Plugging it into Eq. (S1.1), the maximum likelihood estimation reduces to minimizing where Σ (t,s) l , diag(Σ ) = 1, and Ω = DΩD for a diagonal matrix D with D (t,t) kk . It is the same with finding w (t) k to minimize log det(Σ ) under the same constraints, which is the GENVAR procedure of .
S3. Modified STAR procedure. Selectively Traversed Accumulation Rules , or STAR in short, is an iterative procedure designed for user interactive multiple testing in structured settings. The method guarantees FDR control by limiting the analyst's knowledge about the p-values in Eq. (21) of Bong et al. (2021): indeed the analyst sees p-values only after the masking function g(p) = min{p, 1 − p} is applied to them. Based on the masked p-values together with our expectation of contiguous rather than isolated correlated activity (this is the structured settings aspect of the method), at iteration τ we propose a candidate rejection set R τ of non-zero entries (this is the user interactive aspect of the method; details on constructing R τ are provided below). We estimate the FDR for the current rejection set as: where h(p) = 2p and |R τ | is the number of entries in R τ . Then if FDR τ ≤ α for some prechosen α, we stop the algorithm and take R τ as the final rejection set. Otherwise, STAR reveals the unmasked p-values outside of R τ , and we propose the next rejection set R τ +1 as detailed below, such that R τ +1 ⊂ R τ . We iterate this procedure until FDR τ drops below α or R τ = ∅. Algorithm 3 provides full details.
Algorithm 3 STAR  Input: P-values p (t,s) , (t, s) ∈ K, and target FDR level α. In our problem, the constraint set K consists of the p-values that are within dcross of the diagonal of Ω 12 . Output: Rejection set R * with FDR ≤ α. if ( FDRτ ≤ α and Rτ ∈ K) or Rτ = ∅ then 5: Return Rτ and end algorithm. 6: end if Construct R τ +1 ⊆ Rτ based on g p (t,s) : (t, s) ∈ Rτ ∪ p (t,s) : (t, s) ∈ R C τ as per Algorithm 4. 7: end for Suppose we ended iteration τ . We now explain how to construct R τ +1 . We consider all p-values p (s,t) , both masked inside and unmasked outside of R τ , and, as in , we assume that each of them arose from the Beta distribution: with mean 1 + θ (s,t) −1 . We further assume that parameters θ (s,t) are spatially smooth (which we expect since changes in the brain over time are smooth), that is we assume that θ (s,t) = h(s, t) for a smooth function h. We take h to be a Gaussian kernel regression function, and we fit h by maximum likelihood assuming that p-values are independent given their respective θ (s,t) , so that the likelihood is the product of the Beta distributions for the unmasked p-values times the product of the distributions of transformed Beta distributed p-values, with transformation the masking function g(·). We use the quasi EM algorithm in Algorithm 4 to maximize the likelihood. Finally, we obtain R τ +1 by eliminating δ% of p-values in R τ , chosen to be the p-values that are connected to the boundary of R τ or to another eliminated p-value and that have the largest expectations 1 + θ (s,t) −1 , since the corresponding entries Ω (t,s) 12 are more likely than others inside R τ to be zero. The differences between this algorithm and that of  are: (i) we replace the spline basis with Gaussian kernels to model h and (ii) we allow multiple areas of significance instead of a single convex area. the p-value for cluster k by the sampling proportion: Then when a cluster identified by the excursion test overlaps with a cluster identified by STAR, we assign the p-value of the former to the latter. If no excursion test cluster overlaps a STAR cluster, we give it p-value one. In Bong et al. (2021), the p-values for the 3 clusters in Fig. 4c are all smaller than 0.005 and for the 6 clusters in Fig. 6c are 0.024, 0.0015, 0.001, < 0.0005, and the two smallest clusters have p-values > 0.5.
As a check on the validity of this test, Fig. S4.1b shows the quantiles of null p-values versus the theoretical quantiles of the uniform distribution on [0, 1]: the points lie on the diagonal, as expected of a valid test. Fig. S4.1a shows the receiver operating characteristic (ROC) curve of this excursion test; the probability of detecting all connectivity epochs (y-axis) -estimated in 300 datasets simulated as in Section 3.1 of Bong et al. (2021) -are plotted against the cluster family wise error rate (x-axis). The area under the curve is 0.92, meaning that the test has very high power to detect all connectivity epochs. S5. Empirical distributions of Ω 12 entries. In Bong et al. (2021), the p-values in Eq. (21) depend on the normal assumption in Eq. (19). To investigate its validity, we simulated R = 60 datasets, as described in Section 3.1 therein. Fig. S5.1a displays normal Q-Q plots of the R = 60 repeat estimates of Ω (t,s) 12 (Bong et al., 2021, Eq.(19)) for three randomly chosen representative time pairs (t, s), which confirms that Ω S6.1. Effect of the identifiability constraint. Here, we illustrate that the identifiability constraint in Eq. (10) of Bong et al. (2021) that is necessary to run LaDynS can result in lower statistical power to detect connectivity epochs, when applied to data that do not meet the constraint; see Section 3 and Fig. 4 therein. This is not a flaw of the method, only a feature that is shared by many methods that require constraints on parameters. Note also that experimental data may or may not meet the constraint, so statistical power is not necessarily lost. .
We simulated a dataset exactly as described in Section 3.1 of Bong et al. (2021), with r = 0.006, except that we omitted the second term β in the right hand side of Eq. (29) therein, so that the resulting data do not automatically satisfy the identifiability constraint. We ran LaDynS with parameters set as in Section 3.2 of Bong et al. (2021) and display the estimates and inferences in Fig. S6.1. Fig. 4 of Bong et al. (2021) and Fig. S6.1 are similar, except the first is based on data that meet Eq. (10) therein and the second on data that do not. The effect of imposing the constraint is to bias to various degrees the elements of Ω 12 and thus to degrade inference. We cannot quantify this bias, but extensive simulations (not shown) suggest that effects are "smeared" (dampened and spread), and weak, typically statistically insignificant ripples can appear around the effects. In this particular example, the procedure was unable to detect the second cluster entirely, and the third cluster partially. Of course, as with all statistical procedures, larger samples and stronger effects yield more power. (c) Discovered non-zero cross-precision estimates by STAR at nominal FDR 5%.
S6.2. Effect of regularizing the diagonal of Σ. Here, we illustrate that if observed time series are not so smooth as to require regularization of the diagonal of Σ, then STAR controls the FDR much better than when the data are very smooth; see Section 2.4 and Fig. 4 of Bong et al. (2021).
We simulated datasets as described in Section 3.1 of the same paper, except we used the LFP data of Brincat and Miller (2016) for the baseline time series Y 1 and Y 2 in Eq. (29) therein. This dataset consists of multi-electrode recordings from pre-frontal cortex and hippocampus during an associate memory task. We subsampled the 1.25 seconds long trials of this dataset into T = 50 times bins, so that the resulting baselines would be less autocorrelated than the baseline time series (Section 3.1) and experimental data (Section 4) in Bong et al. (2021). We could therefore set λ diag = 0 to run LaDynS. We also used different simulation parameters so the datasets simulated here would have similar properties as the simulated datasets in Bong et al. (2021). Specifically, we set c 0,1 = c 0,2 = 0.025 for the latent autocorrelations (Bong et al., 2021, Eq.(27)), and connectivity intensities r = 0.4, 0.8, 1.2, 1.6, 2.0 (Bong et al., 2021, Eq.(28)). We ran LaDynS with hyperparameters d auto = 5, d cross = 10 and λ auto = λ diag = 0. The penalty on cross-correlation elements, λ cross , was automatically tuned at every run to control false discoveries (see Bong et al., 2021, Section 2.3). Results are in Fig. S6.2, to be compared with Fig. 5 of Bong et al. (2021), which is based on simulated data that require λ diag = 0 to run LaDynS. When the observed time series are not exceedingly smooth, we set λ diag = 0, which does not induce bias in the LaDynS' estimate of Ω 12 or correlations amongst the p-values in Eq. (21) of Bong et al. (2021). In turn, STAR better controls the FDR (see Fig. S6.2a) and has little chance of discovering non-existent clusters (see Fig. S6.2b).  Bong et al. (2021), when LaDynS is applied to very smooth time series. FCDR remains controlled, and FNR and FCNR decrease to zero as the nominal FDR increases, suggesting that STAR has high power to detect all connectivity epochs, unless the data are very weakly connected.
S6.3. Comparison of LaDynS and DKCCA on simulated data. Here we compare La-DynS' output to the output of DKCCA  for data simulated as described in Section S6.2, with ground truth in Fig. 4d of Bong et al. (2021) and connection intensity r = 1.2. For DKCCA, we used the default hyperparameter setting and tuning provided in the package of . Figure S6.3 shows the simulated data true cross-correlation matrix, obtained by inverting the true precision matrix (the simulated time series are not too smooth so the correlation matrix can be inverted without diagonal regularization), and the estimated cross-correlation matrices using by LaDyns and DKCCA. Both methods recover the essential structure of the true cross-correlations; LaDyns does a better quantitative job.
S7. Inference on the change of the factor loadings. The left column of Fig. S7.1 shows the estimated change (the difference) in the factor loadings for the two brain regions from the first bump (150 ms) and the second bump (400 ms) in Fig. 7a of Bong et al. (2021). To check the significance of the distinction, we obtained 95% confidence intervals based on 200 bootstrap estimates. In the bootstrap, N = 3000 trials was sampled from {X 1,[n] , X 2,[n] } n=1,...,N with replacement. Unlike the permutation bootstrap in Section 2.4 of Bong et al. (2021), the resulting dataset maintains not only the autocorrelation within each of X 1 and X 2 but also the correlated activity between them. The 97.5 and 2.5 percentiles of the bootstrap estimated differences are shown at the middle and right column of Fig. S7.1, respectively. In PFC, at 400 ms compared with 150 ms, several electrodes in the upper left corner have substantially larger weights, while the electrodes in the lower middle of the array have substantially smaller weights, and the CIs are clearly separated from 0. There appear to be differences in V4, but they are smaller.
(a) (b) Fig S7.1: Inference on the the factor loading change from 150ms to 400ms. The estimated difference in the factor loadings for (a) V4 and (b) PFC from 150ms to 400ms alongside the 97.5 and 2.5 percentiles of the bootstrap distribution. The .05 and -.05 contours have been added to the plots to show regions, blue above .05, red below -.05, where the differences may be distinguishable from zero. The differences are more dramatic in PFC: there is a region in the upper left corner of the array where the 2.5 percentile is well above .05, and thus separated from zero, while the differences in the patch of (1000, 1500) × (500, 1500) have 97.5 percentile well below -.05.