Integrative learning for population of dynamic networks with covariates.

Although there is a rapidly growing literature on dynamic connectivity methods, the primary focus has been on separate network estimation for each individual, which fails to leverage common patterns of information. We propose novel graph-theoretic approaches for estimating a population of dynamic networks that are able to borrow information across multiple heterogeneous samples in an unsupervised manner and guided by covariate information. Specifically, we develop a Bayesian product mixture model that imposes independent mixture priors at each time scan and uses covariates to model the mixture weights, which results in time-varying clusters of samples designed to pool information. The computation is carried out using an efficient Expectation-Maximization algorithm. Extensive simulation studies illustrate sharp gains in recovering the true dynamic network over existing dynamic connectivity methods. An analysis of fMRI block task data with behavioral interventions reveal sub-groups of individuals having similar dynamic connectivity, and identifies intervention-related dynamic network changes that are concentrated in biologically interpretable brain regions. In contrast, existing dynamic connectivity approaches are able to detect minimal or no changes in connectivity over time, which seems biologically unrealistic and highlights the challenges resulting from the inability to systematically borrow information across samples.

Introduction

There has been a steady development of graph-theoretic approaches to compute dynamic functional connectivity (FC) that is fueled by an increasing agreement that the brain network does not remain constant across time and instead undergoes temporal changes resulting from endogenous and exogenous factors (Filippi et al., 2019). For example, task-related imaging studies have shown that the brain networks will re-organize when the subjects undergo different modulations of the experimental tasks during the scanning session (Chang and Glover, 2010; Lukemire et al., 2020). Similarly, dynamic FC has also been observed during resting-state experiments (Bullmore and Sporns, 2009). These, and other recent studies, have found increasing evidence of underlying neuronal bases for temporal variations in FC which is linked with changes in cognitive and disease states (Hutchison et al., 2013).

Dynamic connectivity approaches involve time-varying correlations derived via graph-theoretic methods, and may be broadly classified into the following categories: (i) change point methods (Cribben et al., 2013; Kundu et al., 2018) that assume stable phases of connectivity inter-spersed with connectivity jumps at unknown locations, which results in piecewise constant connectivity; (ii) Hidden Markov Models (HMMs) involving fast transient networks that are reinforced or revisited over time, which have been applied to electrophysiological data (Quinn et al., 2018) and more recently to fMRI data (Warnick et al., 2018); and (iii) sliding window approaches that enforce temporally smooth correlations (Chang and Glover, 2010; Monti et al., 2014) based on the biologically plausible assumption of slowly varying temporal correlations with gradual changes in connectivity. While sliding window methods are arguably the most widely used, these approaches may be limited by practical issues such as the choice of the fixed window length (Lindquist et al., 2014), although more recent methods use adaptive forgetting estimation generalization that vary the window size to address non-stationarity in the time-series (Monti et al., 2017). On the other hand, change point models and HMMs have the advantage of model parsimony by limiting the distinct number of parameters. However, the performance of these methods often depend on modeling assumptions, and temporal smoothness of connectivity estimates can not be typically ensured. More importantly, since most of these existing approaches typically rely on single-subject data, they often face challenges in terms of detecting rapid changes in connectivity and may result in inaccurate estimates due to a limited information from a single individual.

Essentially, almost the entirety of the existing dynamic connectivity literature has focused on data from single individuals, due to the fact that temporal changes in connectivity are expected to be subject-specific and may not be replicated across individuals. However, recent evidence suggests that combining information across individuals in a group provides more accurate estimates for connectivity (Hindriks et al., 2016), which adheres to the commonly used statistical principle of data aggregation using multiple samples to obtain more robust estimates. Kundu et al. (2018) proposed a sub-sampling approach to compute time varying dynamic connectivity networks using multi-subject fMRI data, which resulted in considerable gains in accuracy under limited heterogeneity across samples, compared to a single-subject analyses. Unfortunately, the variations across samples may not be restricted in many practical settings. To our knowledge, there is a scarcity of carefully calibrated approaches for pooling information across heterogeneous samples in order to accurately estimate a population of (single-subject) dynamic networks. This is perhaps not surprising, given that there are considerable challenges involved. From a methodological perspective, it is not immediately clear how to effectively borrow information across individuals in a data-adaptive manner that also respects the inherent connectivity differences between heterogeneous samples. Similarly when estimating dynamic networks with V brain regions for N individuals each having T time scans, one encounters computational challenges in terms of computing NT distinct V × V connectivity matrices, which is not straightforward for high-dimensional fMRI applications. We note that while there is some related literature on joint estimation of multiple related static networks (Danaher et al., 2014; Guo et al., 2011; Lukemire et al., 2020), it is not straightforward to apply these approaches for the estimation of population of dynamic networks involving moderate to large sample sizes.

In this article, our goal is to develop a fundamentally novel hierarchical Bayesian product mixture modeling (BPMM) approach incorporating covariates (MacEachern, 1999) for estimating a population of dynamic networks corresponding to heterogeneous multi-subject fMRI data. The importance of using covariates to model known stationary networks has already been illustrated in recent literature (Zhang et al., 2018), where the networks are specified in advance. These methods suggest a strong justification for incorporating demographic, clinical, and behavioral covariates when modeling dynamic networks in order to obtain more accurate and reliable estimates (Shi and Guo, 2016). Motivated by these existing studies, the proposed BPMM framework estimates unknown dynamic networks by leveraging covariate information in order to inform the clustering mechanism under the mixture model, which is better designed to tackle heterogeneity across samples that ultimately results in more accurate network estimation. Under the proposed model, subgroups of individuals with similar dynamic connectivity profiles are identified, where the subgroup memberships are also influenced by covariate profiles and change over time in an unsupervised manner that is designed to pool information in order to estimate the dynamic networks. Another appealing feature of the proposed BPMM approach is the ability to report cluster level network summaries that are more robust to noise and heterogeneity in the data. Since the proposed approach clusters samples independently at each time scan guided by covariate information, it is clearly distinct from HMM approaches that instead cluster transient brain states across time scans. To our knowledge, the proposed approach is one of the first to estimate a population of dynamic networks incorporating covariate knowledge by integrating heterogeneous multi-subject fMRI data, which represents considerable advances.

In order to tackle the daunting task of estimating NT connectivity matrices, each of dimension V × V, the proposed approach employs dimension reduction by clustering samples under the mixture modeling framework that translates to considerable computational gains. In particular, the BPMM approach induces model parsimony by reducing the number of unique model parameters from to , where k (<< N) denotes the number of clusters2at the t-time point that is determined in an unsupervised manner. Temporal smoothness in connectivity for each network is also ensured via additional hierarchical fused lasso priors on mixture atoms in the BPMM, which results in gradual changes in connectivity that is biologically meaningful. In scenarios where sharp connectivity changes are anticipated in certain localized time windows (due to changes in experimental design in a block task experiment, or other exogeneous or endogenous factors), one may estimate these connectivity change points via a post-processing step that involves applying the total variation penalty (Vert and Bleakley, 2010) to the dynamic connectivity estimates under the proposed approach. Additional post-processing steps involving a K-means algorithm are also proposed to identify subgroups of individuals with similar dynamic connectivity patterns consolidated across time, which is particularly useful in terms of obtaining insights related to heterogeneity. Figure 1 provides a visual illustration of the proposed approach.

Fig. 1.

A schematic diagram illustrating the proposed dynamic pairwise correlation method. A mixture prior with H = 3 components is used to model dynamic correlations, where the mixture weights are modeled using covariates. The resulting networks at each time scan for each sample are allocated to one of the H clusters representing distinct network states that are represented by red, orange and blue cubes. Although the proposed method does not cluster transient states across time, the simplified representation in the Figure illustrates the similarity of brain states contained in identical colored cubes across the experimental session. Such temporal smoothness of the network is imposed via hierarchical fused lasso priors on the mixture atoms. Once, the dynamic FC is estimated, a post-processing step using K-means (Section 2.2) is applied to compute sub-groups of samples that exhibit similar dynamic connectivity patterns summarized across all time scans. The subgroups are represented by the circle, pyramid, triangle and inverted triangle shapes in the Figure and correspond to different modes of dynamic connectivity with different number of brain states represented by different patterns within each shape. The connectivity change points for each individual, as well as at a cluster level, are computed via another post-processing step that employs a group fused lasso penalty (Section 2.3). The method reports both individual and cluster-level network features.

The proposed BPMM is developed for dynamic pairwise correlations as well as dynamic precision matrices, which provide complementary interpretations of dynamic connectivity. In particular, pairwise correlations encode connections between pairs of nodes without accounting for the effects of third party nodes, whereas partial correlations report association between nodes conditional on the effects of the remaining network nodes. While our goals do not involve assessing the merits of one approach over the other (see Smith et al. (2013) for a review), the proposed development is designed to provide users with an option to implement either approach as desired and suitable for respective applications. We develop an effcient Expectation-Maximization (EM) algorithm to implement the dynamic pairwise correlation method separately for each edge, and another EM algorithm for dynamic precision matrix estimation that simultaneously involves all network nodes. We perform extensive simulations to evaluate the performance of the proposed method in contrast to existing approaches that involved a variety of dynamic network structures. The proposed methods were also used to investigate dynamic functional connectivity changes due to a high intensity, aerobic exercise ‘spin’ intervention when compared to a non-aerobic exercise, control intervention, which were administered to a heterogeneous group of sedentary adults who performed a fMRI block task experiment. Our goals are to provide connectivity insights that are complimentary to previous activation-based findings from the same study (Nocera et al., 2017), but involves analytic challenges due to the short duration of the fixation and task blocks that induce rapid connectivity changes which are usually diffcult to detect via existing methods.

The rest of the article is structured as follows. Section 2.1 develops the proposed approach for dynamic pairwise connectivity that performs edge-wise analysis (denoted as integrative dynamic pairwise connectivity with covariates or idPAC), Section 2.2 develops a post-processing clustering approach for identifying subgroups of individuals with similar dynamic connectivity patterns consolidated across time scans based on the estimated dynamic correlations under the idPAC method, and Section 2.3 outlines a post-processing strategy for estimating network change points. Section 3 extends our framework to estimate dynamic precision matrices that uses all the network nodes simultaneously for network estimation (denoted as integrative dynamic precision matrix with covariates or idPMAC). Section 4 develops a computationally efficient EM algorithm to implement the proposed idPAC and idPMAC approaches, and describes choices for tuning parameters. Section 5 reports results from extensive simulation studies, and Section 6 reports our analysis and results from the block-task fMRI experiment. Additional discussions are provided in Section 7. Throughout the article, we will use BPMM to denote the overall Bayesian product mixture modeling framework used for developing the idPAC and idPMAC approaches, as appropriate.

Methods

In this section, we propose a novel approach for estimating a population of dynamic networks using heterogeneous multi-subject fMRI data with the same number of brain volumes across all individuals. For modeling purposes, we will assume that the demeaned fMRI measurements are normally distributed with zero mean (Kundu et al., 2018) at each time scan, and that pre-whitening steps have been performed to minimize temporal correlations (see Supplementary Materials for details). We will fix some notations here. Suppose fMRI data are collected for T scans and V nodes (voxels or regions of interest) for N individuals. Denote the fMRI measurements across all the nodes at time point t as , and denote the V × T matrix of fMRI measurements for the ith individual as ( that has the tth column as . Further, denote the vector of q × 1 covariates as for the ith sample, and represent the collection of fMRI data matrices across all individuals as .

In what follows, the idPAC method for pairwise correlations (Section 2.1) and idPMAC method for partial correlations (Section 3), will involve a combination of likelihood terms and priors on the model parameters that are combined into a posterior distribution, which is used to estimate model parameters. The posterior distribution for parameter θ given data Y is defined as using Bayes theorem, where L(Y|θ) denotes the data likelihood given the parameter value θ, π(θ) represents the prior on θ under the Bayesian model, and is the marginal likelihood after integrating out all possible values of θ. Full details of the posterior distributions are provided in the Appendix.

Dynamic connectivity via pair-wise correlations

Let the unknown dynamic pairwise correlation of individual i be denoted as , and the corresponding Fisher-transformed pairwise correlations be denoted as . We propose a Bayesian hierarchical approach that models the edge-wise dynamic correlations, using data from multiple individuals. We propose the following model for edge (j, l), and jointly for t = 1, …, T and i = 1, …, n, where |⋅| denotes the L1 norm, denotes the residual variance in the likelihood term, the Fisher-transformed correlations are modeled under a mixture of Gaussians prior having H components denoted as , with the prior probability for the hth mixture component denoted as that depends on covariates, such that for all captures the (unknown) variability of the pairwise correlations under the mixture prior specification, and N(μ, Σ) denotes a multivariate Gaussian distribution with mean μ and V × V covariance matrix Σ. Under a hierarchical Bayesian specification, is estimated under the conjugate Gamma prior with shape and scale parameters respectively. The mixture prior specifies that for any given time scan t, the functional connectivity for each individual can take values revolving around any one of the H mixture atoms denoted by , that are themselves unknown and modeled under a fused lasso prior as in (1). These values are realized with respective prior probabilities that are modulated via covariates with effect sizes respectively, where with fixed as the reference group.

Modeling mixture atoms via fused lasso:

The mixture atoms are modeled under a fused lasso prior in (1) that encourages temporal smoothness of pairwise correlations by assigning small prior probabilities for large changes in the values between consecutive time scans. Although temporal smoothness in correlations is encouraged, the Bayesian approach is still equipped to accommodate sharp jumps in connectivity that may arise due to changes in experimental design or other factors. Such connectivity jumps are detected using a post-processing step (see Section 2.3) applied to the estimated dynamic connectivity under the proposed model.

Modeling mixture weights via covariates:

In order to effectively tackle heterogeneity, we incorporate supplementary covariate information when modeling the mixture weights under our mixture modeling framework in (1). By incorporating covariate information, the model is designed to achieve more accurate identification of clusters, which then naturally translates to improved estimates for dynamic FC at the level of each individual. In particular, we model via a Multinomial Logistic regression (Engel, 1988), where (with β = 0) represents the vector of unknown regression coefficients that control the contribution of the covariates to the mixture probabilities for the Hth component (h = 1, …, H − 1), in contrast to the Hth component. A large value of these regression coefficients implies increased importance of the corresponding covariate with respect to modeling a particular edge under consideration, whereas for all , indicates spurious covariates unrelated to the dynamic pairwise correlations.

The multinomial logistic regression model incorporating covariates suggests that the log-odds for each component , can be expressed as a linear combination of covariates. When two or more samples have similar covariate information, the prior specification in (1) will encourage similar mixture components to characterise the dynamic connectivity for all these samples that will result in analogous connectivity patterns. However the posterior distribution (that is used to derive parameter estimates) should be flexible enough to accurately estimate varying connectivity patterns between individuals even when they share similar covariate values, by leveraging information present in the data (as evident from extensive numerical studies in Section 5).

Role of clustering in tackling heterogeneity and pooling information:

Under model (1), each sample will be assigned to one of the H clusters at each time scan in an unsupervised manner and guided by their covariate profiles in order to model the edge-level dynamic connectivity. Due to independent clustering at each time scan, these cluster configurations change over the experimental session in a data-adaptive manner to characterize connectivity fluctuations across individuals. Such time scan specific clusters represent subgroups of individuals with similar connectivity profiles over a subset of time scans, which are learnt by pooling information across all samples within a cluster. Here, it is important to note that model (1) does not impose identical dynamic connectivity across all time scans between multiple individuals (that is biologically unrealistic), but instead encourages common connectivity patterns within subgroups of samples for a subset of time points that are learnt in a data-adaptive manner. Hence, the proposed method is designed to result in more accurate estimation compared to a single subject analysis that is not equipped to pool information across samples or a group level analysis that does not account for within sample heterogeneity. We note that although the estimation is performed separately for each edge, the connectivity estimates across all edges are consolidated to obtain connectivity change point estimates (Section 2.3) or identify subgroups with common dynamic connectivity profiles (Section 2.2).

Post-processing steps for sub-group detection

In practical neuroimaging applications, it is often of interest to detect dissimilar modes of dynamic connectivity patterns that are embodied by distinct subgroups of individuals who also differ in terms of demographic or clinical characteristics, or other factors. For example in our fMRI task study, one of the objectives is to assess variations in dynamic connectivity with respect to subgroups of samples that were assigned different interventions, and who also had varying demographic characteristics. Instead of comparing network differences between pre-specified subgroups that are likely to contain individuals with heterogeneous connectivity patterns, it is more appealing to develop a data-adaptive approach to identify subgroups that comprise individuals with homologous dynamic connectivity, and then examine connectivity variations across such subgroups and how these variations are related to intervention and other factors of interest. When estimating these subgroups, we do not require identical dynamic connectivity patterns for all individuals within subgroups, but rather expect them to have limited network differences in terms of edge strengths and connectivity change points. An inherently appealing feature of subgroup detection is that is allows one to compute cluster level change points and other aggregate network features (see Section 2.3) which are more reproducible in the presence of noise and heterogeneity, compared to a single-subject analysis. Subgroup level network summaries may be particularly beneficial in certain scenarios such as fMRI block task experiments where it may be challenging for single-subject analyses to detect rapidly evolving network features induced via quick transitions between rest and task blocks within the experimental design.

We propose an approach that consolidates the time-varying clusters of samples under the BPMM approach to detect subgroups which comprises samples with similar network-level dynamic connectivity patterns. In order to identify these subgroups, we first create a N × N similarity matrix that measures the propensity of each pair of samples to belong to the same cluster over the experimental session. This matrix is created by examining the proportion of time scans during which a pair of samples belonged to the same cluster across the experimental session, averaged across all edges. Once this similarity matrix has been computed, a K-means algorithm is applied to identify clusters of samples that exhibit similar dynamic connectivity patterns across the experimental session. The number of clusters K is determined using some goodness of fit score such as the elbow method (Thorndike, 1953), or it is fixed as the maximum number of mixture components (H) under the BPMM approach. Finally, we note that the subgroup identification step is not strictly needed under the proposed BPMM framework for dynamic network estimation, but it is an optional analysis that can be used to identify cluster-level network features in certain scenarios of interest.

Post-processing steps for connectivity change point estimation

The estimated dynamic correlations can be used to detect connectivity change points in scenarios involving sharp changes in the network during the session, such as in fMRI task experiments. Our strategy involves computing change points for each individual network (a) at the edge level that captures localized changes; and (b) at the global level that captures major disruptions in connectivity over the entire network. We compute the change points using the total variation penalty (Vert and Bleakley, 2010) that was also used in CCPD approach by Kundu et al. (2018). However the proposed methods are distinct from the two-stage CCPD approach; the latter estimates connectivity change points based on empirical time-varying connectivity measures in the first stage, and then in the second stage, computes piecewise constant networks conditional on the estimated change points that represent connectivity jumps. In contrast, proposed method pools information across samples in order to first estimate dynamic correlations that does not depend on change points and can vary continuously over time, and subsequently uses a post-processing step to compute connectivity change points without requiring piecewise constant connectivity assumptions. An appealing feature of the proposed mixture modeling framework guided by covariates is that it is more suitable for tackling divergent dynamic connectivity across samples, in contrast to empirical correlations under the CCPD approach.

Denote the vector of estimated correlations over all edges for the ith individual and at time scan t as . Then the functional connectivity change points for the ith individual may be estimated using connections across all edges via a total variation norm penalty that is defined as . In particular, the following penalized criteria is used as in Kundu et al. (2018) for detecting network level connectivity change points: where λ represents the penalty parameter and represents the piecewise constant approximation to the time series of correlations at time point t for the ith individual that also assumes the presence of an unknown number of connectivity jumps. The first term in (2) measures the error between the observed correlations and the piece-wise constant connectivity, while the second term controls the temporal smoothness of correlations for V(V − 1)/2 edges. The increment in the second term becomes negligible when the multivariate time series does not change significantly between times t and t + 1, but it takes large values corresponding to significant connectivity changes. The network change points computed via (2) represent global changes functional connectivity resulting from a subset of edges that exhibit large connectivity changes. It is important to note that not all edges are expected to exhibit changes at these estimated change points. When it is of interest to compute edge-level connectivity change points, one can simply use criteria (2) separately for each edge, so that the total variation term translates to the L1 penalty. However, it is important to note that edge-level connectivity changes represent granular fluctuations that are typically more challenging to detect in the presence of noise in fMRI.

The number of change points is determined by the penalty parameter λ, with a smaller value yielding a greater number of change points and vice-versa. Tibshirani and Wang (2008) proposed an estimate of λ based on a pre-smoothed fit of a univariate time series using a lowess estimator (Becker et al., 1988). We adapt this approach for a multivariate time series to obtain an initial estimate for λ, and then propose post-processing steps to tune this estimate in order to obtain change points, as in the CCPD approach in Kundu et al. (2018). Full details for these steps are provided in Section 3 of Supplementary Materials.

Cluster-level connectivity change point estimation:

For fMRI task experiments involving multiple subjects, subgroups of individuals are expected to share analogous dynamic connectivity patterns with limited variations across samples, as discussed in Section 2.2. The proposed total variation penalty norm in (2) is equipped to leverage information across samples within a cluster for identifying cluster level change points, which reflect aggregated dynamic connectivity changes across all samples within a cluster at the global network level. These cluster level connectivity changes are obtained by aggregating the change points obtained via (2) applied separately to each sample within the cluster, and then choosing those change points that show up repeatedly within the cluster. One can define a threshold such that all change points that appear with a high frequency (above the chosen threshold) across samples within the cluster are determined to represent cluster level change points (Kundu et al., 2018). We note that under the proposed method, it is entirely possible for individuals within a cluster to have unique connectivity changes in addition to the common cluster level change points, which reflect within sample heterogeneity. In our experience, this method typically works well in accurately recovering aggregated cluster-level connectivity changes, in certain scenarios such as block task experiments, or more generally in the presence of subgroups of individuals with similar dynamic connectivity patterns.

Extension to dynamic precision matrix estimation

We now propose a mixture model for dynamic precision matrix estimation (idPMAC) that looks at the totality of all nodes in the network, in contrast to the edge-wise analysis in Section 2.1. While the idPMAC also uses a mixture modeling framework, it is fundamentally distinct compared to the idPAC method in Section 2.1, with respect to the manner in which the mixture prior is specified and in terms of how the network edges are constructed and interpreted. The proposed approach estimates the network by computing the V × V precision matrix involving V(V − 1)/2 distinct partial correlations that are learnt by borrowing information across V nodes at each time scan. The partial correlations measure interactions between pairs of regions after removing the influence of third party nodes, which is successful in filtering out spurious correlations. Hence a zero partial correlation between two nodes implies conditional independence. The proposed idPMAC approach enables one to report graph-theoretic network summary measures that capture important patterns of network information transmission (Lukemire et al., 2020), which may not be straightforward to report using pairwise correlations (Smith et al., 2012).

Denote the V × V precision matrix over all nodes for the ith individual at the tth time point as , and note that the partial correlation between nodes k and l is given directly as (ignoring the subject-specific and time-scan specific notations). We propose a Gaussian graphical model involving product mixture priors as: for i = 1, …, N, t = 1, …, T, where , the space of symmetric positive definite matrices, E(α) denotes the Exponential distribution with scale parameter α, and denotes the vector of (V − 1) off-diagonal elements corresponding to the vth row of that are modeled using a mixture of multivariate Gaussians prior. Specifically, the dynamic connectivity at time scan t is likely to be characterised via the hth mixture component with prior probability depending on covariates, where the prior mean and precision for this unknown mixture component is given by and respectively. The idPMAC approach in (3) specifies independent mixture priors on the set of all edges related to each node and at each time scan, while ensuring that the precision matrices are symmetric and positive definite. Full details for the computational steps are presented in Section 4.

Modeling mixture atoms:

Under a hierarchical Bayesian specification, the mixture atoms or component-specific means are themselves unknown and modeled via a fused lasso prior, which encourages temporal homogeneity of partial correlations by assigning small prior probabilities for large changes in the values. In addition, systematic changes in connectivity reflected by sharp jumps may be still identified via a post-processing step in Section 2.3. The unknown prior variance on mixture components is assigned a Gamma prior and is estimated adaptively via the posterior distribution.

Modeling mixture weights via covariates:

The node level mixture weights incorporating covariates are modeled via a Multinomial Logistic regression with unknown covariate effects corresponding to time scan t that are assigned Gaussian priors, and we fix , as the reference group. The prior in (3) encourages similar clustering configurations resulting in analogous time-varying partial correlations for individuals with similar covariate profiles. However in the presence of heterogeneity, the posterior distribution under the idPMAC method is still able to identify divergent dynamic connectivity patterns even among individuals with same covariate profiles (see numerical studies in Section 5).

Role of clustering in tackling heterogeneity and pooling information:

Under model (3), each column of the precision matrix is assigned to one of the H clusters at each time scan in an unsupervised manner. Hence, the mixture modeling framework allows subsets of rows/columns of to have the same values depending on their clustering allocation at each given time scan, which is an unique feature under the idPMAC approach that is not shared by the idPAC method. This feature results in robust estimates by pooling information across nodes and samples to estimate common partial correlations, and is a necessary dimension reduction step for scenarios involving large networks. For example, all weak or absent edges can be subsumed into one cluster which yields model parsimony. In addition, divergent connectivity patterns are captured via distinct time-varying clustering configurations across individuals as derived from the posterior distribution, which accommodates heterogeneity. Hence, the clustering mechanism under the idPMAC method not only enables dimension reduction, but also provides a desirable balance between leveraging common connectivity patterns within and across networks and addressing inherent network differences across individuals.

Post-processing steps:

The post-processing steps for sub-group detection and connectivity change point estimation under the idPMAC approach can be applied in a similar manner as outlined in Sections 2.2 and 2.3. They proceed by replacing the estimated pairwise correlations with estimated dynamic partial correlations derived under the idPMAC approach in the K-means algorithm and the fused lasso criteria (2) in Sections 2.2 and 2.3 respectively.

Computational details for parameter estimation

Although one can use Markov chain Monte Carlo (MCMC) to sample the parameters from the posterior distribution, we use a maximum-a-posteriori or MAP estimators for our purposes in this article that bypasses the computational burden under a MCMC implementation. The MAP estimators are obtained by maximizing the posterior distribution for the model parameters and are derived via the Expectation-Maximization or EM algorithm. The EM algorithm is scalable to high-dimensional fMRI applications of interest that requires one to compute N × T distinct dynamic networks each involving V × V connectivity matrices. Table 1 provides a list of model parameters to be estimated via the EM steps for both the dynamic pairwise correlations and the dynamic precision matrix estimation methods.

Table 1

Summary for all model parameters under the dynamic pairwise correlations and the dynamic precision matrix methods. MC E-step refers to Monte Carlo E-step.

Notation	Description DATA	Update
Y	fMRI scanning data for all individuals	observed
ylt(i)	observed fMRI data for individual i, node l, time point t	observed
x i	covariate information for individual i	observed
σy2	Variance of ylt(i) that is empirically estimated	Pre-Fixed
Σβ	prior covariance for covariate effects β	Pre-Fixed
H	Number of components in mixture	Pre-Fixed
	Dynamic Pairwise Correlation
ρjl,ti	pairwise corr for edge (j, t) at time t (γjl,t(i) = Fisher-transformed ρjl,ti)	M-step
γh,jlt∗	mean of hth component for edge (j, t) at time point t	M-step
σγ,h2	Variance for hth mixture component	M-step
β_h,jlt	unknown regression coefficient used for modelling ξ_h,jlt(x_i)	M-step
ψh,jlt(i)	posterior probability of γjl,t(i) taking values from hth component	E-step
	Dynamic Precision Matrix
ωh,t∗	hth mixture component	M-step
σω,h2	mixture variance of hth component	M-step
βh,t	unknown regression coefficient used in Multinomial Logistic regression	M-step
ψh,vt(i)	posterior prob ωv,t(i) of taking values from component h	E-step
Ωt(i)	precision matrix for individual i at time point t	MC E-step

EM algorithm for pair-wise dynamic connectivity

EM Algorithm:

Denote the matrix containing the fMRI time series for the lth node as where represents the fMRI observations across all samples for node l and time scan t. Further, denote Δ as a latent indicator variable for the hth mixture component (that is not observed and is imputed in the proposed EM algorithm) and finally, denote by the collection of all model parameters under the specification (1) corresponding to edge . Note that under the proposed multinomial logistic regression model for incorporating covariates in (1), one has an equivalent specification under the binary latent variables distributed as where denotes a multinomial distribution with probability vector and one can marginalize out to recover the original specification in (1). The EM algorithm uses the augmented log-posterior derived in the Appendix involving the above latent mixture indicators, to compute MAP estimates for the model parameters by iteratively applying the Expectation (E) and Maximization (M) steps. The latent indicators are imputed via the E-Step by using the posterior probability of taking values from the hth mixture component, which is denoted by and updated as:

E-step :

Compute the posterior expectation for the latent cluster membership indicators as , where denotes the normal density with mean and variance .

The remaining parameters are updated via M-steps using closed form solutions except that is updated using Newton-Raphson steps. These M-steps comprise several mathematically involved derivations and are detailed in the Appendix. The E and M steps are repeated till convergence, which occurs when the absolute change in the log-posterior between successive iterations falls below a certain threshold (we use 10−4 in our implementation).

EM Algorithm for dynamic precision matrix estimation

Let us denote the collection of all the precision matrices as , and as the (V − 1)-dimensional vector of fMRI measurements at time scan t over all nodes except node v. The prior on the precision matrix can be expressed as , with the corresponding prior distributions being defined in (3). Denote by |⋅|1, the element-wise L1 norm, denote to represent the conditional variance corresponding to the fMRI measurements for the vth node given all other nodes, and let and respectively denote the diagonal and the vector of off-diagonal elements of the vth row in . Moreover use det(A) to denote the determinant of the matrix A, and write as the matrix of cross-products of the response variable, where and denote the vth diagonal element and the off-diagonal elements for the vth row respectively. Introduce latent indicator variables that follow a multinomial distribution with probability vector such that .

Denote by , the (V − 1) × (V − 1) obtained by deleting the vth row and column from . The EM algorithm uses an E step for the latent mixture indicators, as well as a Monte Carlo E step that samples from the posterior distribution in order to obtain estimates for the precision matrix. These steps are described below:

E-step for mixture component indicator :

For v = 1, …, V, use the expression: , where denotes the probability density function for the (V − 1)-dimensional normal density with mean and variance as respectively.

Monte Carlo E-step for precision matrix:

We use an E-step to update the precision matrix that computes the posterior mean by averaging MCMC samples drawn from the posterior distribution, which is equivalent to a Monte Carlo EM method (Wei and Tanner, 1990). We use this Monte Carlo approximation for the conditional expectation since it provides a computationally efficient approach to sample symmetric positive definite precision matrices via closed form posteriors. The posterior distribution for the precision off-diagonal elements are given as , where is the posterior covariance. Moreover, writing the diagonal precision matrix elements are updated via the posterior where α is pre-specified. The above steps can be alternated to sample positive definite precision matrices as in Wang (2012), and we draw several MCMC samples and average over them to approximate the conditional expectation.

The remaining parameters are updated via closed form expressions under the M step, which involve mathematically involved derivations and are detailed in the Appendix. The algorithm iterates through the E and M steps until convergence.

Tuning parameter selection

Certain tuning parameters in the BPMM need to be selected properly or pre-specified, in order to ensure optimal performance. For both dynamic pair-wise correlations and precision matrix estimation, λ is the tuning parameter used in fused lasso penalty for the mixture atoms that controls the temporal smoothness of the dynamic connectivity. We choose an optimal value for λ over a pre-specified grid of values, as the value of the tuning parameter that minimizes the BIC score. In model (1) for the dynamic pairwise correlation, the is also pre-specified as the initial mean variance over all edges and across all samples. Moreover when updating covariate effects, Σ is pre-fixed as a diagonal matrix with the diagonal terms as 100, although it is also possible to impose a hierarchical prior on Σ and update it using the posterior distribution. Extensive simulation studies revealed that the proposed approach is not sensitive to the choices of Σ as long as the variances are not chosen to be exceedingly small. Other hyper-parameters in the hierarchical Bayesian specification include α in the prior on the precision matrices (chosen as in Wang (2012)), and that results in an uninformative prior on the mixture variance.

The number of mixture components H also needs to be chosen appropriately. On the one hand, a large value of H may be used to address inherent heterogeneity, but it will also increase the running time and may generate redundant clusters that overcompensates for the variations across samples. On the other hand, a small value of H may restrict the approach to overlook connectivity variations across individuals, resulting in inaccurate estimates. One may use a data adaptive approach to select H in certain scenarios where it is reasonable to assume that the dynamic connectivity can be approximated by piecewise constant connectivity. In such cases that potentially involve block task experiments (Kundu et al., 2018), one can evaluate criteria (2) separately for each individual under different values of H, and fix the optimal choice as that which minimizes the average value of the criteria (2) across all individuals. Based on extensive empirical studies, we noticed the need for larger values for H when fitting the model for cases involving large number of nodes and samples.

Numerical experiments

Simulation set-up

Data generation:

We generate observations from Gaussian distributions with sparse and piecewise constant precision matrices that change at a finite set of change points. Moreover, the network change points are generated based on covariate information where individuals with identical covariates have partially overlapping connectivity change points. Broadly, we use the following few steps to generate the data, each of which is described in greater detail in the sequel: (i) generate a given number of change points for each subject using corresponding covariate information; (ii) conditional on the generated change points, piecewise constant networks are simulated such that the connectivity changes occur only at the given change points; (iii) conditional on the network for a given state phase, a corresponding positive definite precision matrix is generated for each time scan where non-zero off-diagonal elements represent edge strengths and zero off-diagonals represent absent edges; and (iv) the response variable for a given time point is generated from a Gaussian distribution having zero mean and the precision matrix in step (iii). Four clusters are created with 10 samples each, where the samples with each cluster have the same number of connectivity change points, common state phase specific networks and identical covariate values. However within each cluster, there are differences in locations of connectivity change points and the network edge strengths are free to vary across individuals even when they share the same network structure. All samples in the first two clusters have 3 connectivity change points each, whereas the samples in the other two clusters have 4 change points, out of a total of T = 300 time scans.

Conditional on the change points in step (i), several types of networks are constructed for each state phase in step (ii) that include: (a) Erdos Renyi network where each edge can randomly appear with a fixed probability; (b) small-world network, where the mean geodesic distance between nodes are relatively small compared with the number of nodes and which mimics several practical brain network configurations; and (c) scale-free network that resembles a hub network where the degree of network follows a power distribution. Given these networks, the corresponding precision matrix was generated in step (iii) by assigning zeros to off-diagonals for absent edges, and randomly generating edge weights from uniform [−1,1] for all important edges. To ensure the positive definiteness, the diagonal values of the precision matrix were rescaled by adding the sum of the absolute values of all elements in each row with one. Finally, the response variables were generated either (a) independently at each time point via a Gaussian graphical model, or (b) via a vector autoregressive (VAR) model where the response variables are auto-correlated across time. In both cases, sparse time-varying precision matrices having dimensions V = 40, 100, were used. The ‘VARM’ function in Matlab was used to generate temporally correlated observations under a lag-1 VAR model, where the elements in autocorrelation matrix were generated from a uniform random variable with range (−0.2, 0.2).

We generated two binary features that resulted in four distinct covariate configurations, i.e. (0,0), (0,1), (1,0), (1,1), and all samples with identical covariates were allocated to the same cluster. In addition, we also evaluated the performance of proposed method in the presence of spurious covariates that are not related to dynamic connectivity patterns. Specifically, we introduced anywhere between 1 to 8 spurious covariates for each sample (in addition to the two true covariates described earlier), which were randomly generated using uniform as well as from random normal distributions. We then investigated the performance of the proposed approach over varying number of spurious covariates. While the proposed approach is expected to work best in practical experiments involving a carefully selected set of covariates that influence dynamic connectivity patterns, our goal was also to investigate the change in performance as the number of spurious covariates increase.

Competing methods:

We perform extensive simulation studies to evaluate the performance of the proposed approach, and compare the performance with (a) change point estimation approaches such as the CCPD (Kundu et al., 2018) that can estimate single subject connectivity using multi-subject data in the presence of limited heterogeneity, and the dynamic connectivity regression (DCR) approach for single subjects proposed in Cribben et al. (2013); (b) an empirical sliding window based approach (SD) and the model-based SINGLE (Monti et al., 2014) method that uses sliding window correlations; and (c) a covariate-naive version of the proposed approach using the methods in Section 2.1 and Section 3 (denoted as BPMM-PC and BPMM-PR respectively) that employs a multinomial distribution to model the mixture weights without covariates. While methods in (a) and (c) are designed to report connectivity change points, we were also able to compute change points under the sliding window approaches in (b) by applying a post-processing step in (2) on the estimated sliding window correlations. Moreover, the data under the VAR case were prewhitened via an autoregressive integrated moving average (ARIMA) before fitting the proposed models. In particular, the ‘auto.arima’ in R was used to prewhiten the raw data, which yielded residuals that were subsequently used for analysis (more details provided in Supplementary Materials). The SINGLE approach was implemented using the python implementation in pySINGLE, while all other methods were implemented in Matlab.

Performance metrics:

We evaluate the performance of different approaches in terms of different metrics. First, we investigated the accuracy in recovering true connectivity change points at the network and edge level for each sample, using sensitivity (defined as the proportion of truly detected change points or true positives), as well as the number of falsely detected change points or false positives. In addition, the performance of the network connectivity change points at the cluster level was also evaluated by comparing the true connectivity change points for each sample within the cluster with the aggregated cluster level change points. We note that since there were variations in connectivity change points within each cluster, false positive change points are to be expected under any estimation approach; however our goal is to evaluate how well these false positives are controlled and the sensitivity in detecting true change points under different methods. In addition, we also evaluated accuracy in terms of estimating the strength of connections that is computed as a squared loss (MSE) between the estimated and the true edge-level pairwise correlations. The pairwise correlations corresponding to dynamic precision matrix approaches for computing MSE were obtained by inverting the respective precision matrices.

In order to evaluate the accuracy in dynamic network estimation, we computed the F-1 score defined as 2(Precision × Recall)/(Precision + Recall), where Precision=T P /(T P + F P) is defined as the true positive rate, and Recall=T P /(T P + F N) represents the sensitivity in estimating the edges in the network. Here, T P, F P, F N, refer to the number of true positive, false positive, and false negative edges that are obtained via binary adjacency matrices derived by thresholding the estimated absolute partial correlations. We employed reasonable thresholds (0.05) that are commonly used in literature (Kundu et al., 2018). In contrast, it was not immediately clear how to choose such thresholds for pairwise correlations given the fact that they tend to be typically larger in magnitude and have greater variability. Hence, we did not report F-1 scores corresponding to pairwise correlations, although one could do so in principle by choosing suitable thresholds to obtain binary adjacency matrices. Finally, we also evaluated the clustering performance in terms of the clustering error (CE) and Variation of Information (VI). CE (Patrikainen and Meila, 2006) is defined as the maximum overlap between the estimated clustering with the true clustering, whereas VI (Meilǎ, 2007) calculates the entropy associated with different clustering configurations.

Results

The performance in terms of recovering the true clusters of subjects is provided in Table 2, in the presence of two covariates that are both related to the true connectivity changes. It is clear from the results that incorporating covariate information results in near perfect recovery of the clusters, in contrast to the covariate-naive version of the method. For V = 100, the dynamic pairwise correlation approach seems to have a slightly higher accuracy in terms of cluster recovery compared to the dynamic precision matrix approach when data were generated from a VAR model. However when covariates are not included, the BPMM-PR method has greater clustering accuracy compared to the BPMM-PC approach, since the former is able to pool information across the whole network to inform the clustering mechanism, in contrast to an edge-byedge analysis under BPMM-PC. Table 3 reports the accuracy in recovering the true network-level change points under the proposed approaches at the level of the estimated clusters, as per discussions in Section 2.3. In this case, both idPAC and idPMAC methods are shown to have near perfect recovery of the true network connectivity change points when data were generated under GGM, and high sensitivity when data were generated under VAR. Moreover when using data from a VAR model, the idPAC method has a comparable or higher sensitivity but also higher false positives for V = 100 in terms of detecting connectivity change points at the cluster level, compared to the idPMAC method. We note that although all samples within a cluster had identical covariate information, the proposed approach was able to accommodate within cluster connectivity differences that is evident from low false positives and high sensitivity when estimating cluster level change points. Moreover as seen from Tables 4–5, the accuracy in recovering cluster level connectivity change points is considerably higher than the corresponding results at the level of individual networks. These results indicate the usefulness of aggregating information when it is reasonable to assume the existence of subgroups of individuals who share some similar facets of dynamic connectivity.

Table 2

Clustering performance under different network types. GGM implies that the Gaussian graphical model was used to generate temporally uncorrelated observations; VAR implies a vector autoregressive model that was used to generate temporally dependent observations. For the VAR case, the observations were pre-whitened before fitting the model.

	idPAC				BPMM-PC
	V=40		V=100		V=40		V=100
	CE	VI	CE	VI	CE	VI	CE	VI
GGM+Erdos-Renyi	0	0	0	0	0.64	1.93	0.62	2.19
GGM+Small-world	0	0	0	0	0.57	1.92	0.71	2.23
GGM+Scale-free	0	0	0	0	0.63	2.01	0.66	2.19
VAR+Erdos-Renyi	0	0	0	0	0.61	1.93	0.67	1.97
VAR+Small-World	0	0	0	0	0.59	1.88	0.61	1.90
VAR+Scale-Free	0	0	0	0	0.61	1.78	0.61	1.93
	idPMAC				BPMM-PR
	V=40		V=100		V=40		V=100
GGM+Erdos-Renyi	0	0	0	0	0.43	1.41	0.54	1.59
GGM+Small-world	0	0	0	0	0.41	1.41	0.51	1.68
GGM+Scale-free	0	0	0	0	0.43	1.49	0.60	1.78
VAR+Erdos-Renyi	0.08	0.25	0.04	0.17	0.54	1.51	0.66	1.88
VAR+Small-World	0	0	0.03	0.14	0.48	1.47	0.58	1.91
VAR+Scale-Free	0	0	0.04	0.11	0.49	1.42	0.63	1.75

Table 3

Cluster-based network change point estimation under the proposed approaches, assuming that all samples within a particular cluster have the same number and similar location of change points, with a limited degree of heterogeneity in functional connectivity. If this assumption holds, then the cluster level network change point estimation provides greater accuracy compared to the estimated change points at the level of individuals as reported in subsequent Tables.

	idPAC				idPMAC
	V=40		V=100		V=40		V=100
	sens	FP	sens	FP	sens	FP	sens	FP
GGM+Erdos-Renyi	1	2.15	0.99	1.58	0.97	3.94	0.99	3.18
GGM+Small-world	0.97	2.11	1	1.59	0.99	4.18	0.98	3.17
GGM+Scale-free	0.99	2.09	1	1.37	1	3.91	0.97	3.09
VAR+Erdos-Renyi	0.91	3.71	0.88	3.66	0.87	3.47	0.87	2.89
VAR+Small-world	0.84	3.44	0.8	3.09	0.82	3.45	0.81	2.98
VAR+Scale-free	0.88	3.29	0.84	3.68	0.85	3.3	0.81	3.01

Table 4

Results under the dynamic pair-wise correlation approaches for network and edge-level connectivity change-point estimation (Edge CP) accuracy and network changepoint (Network CP) estimation accuracy for V = 40, 100. GGM and VAR correspond to data generated from Gaussian graphical models and vector autoregressive models. Significantly improved metrics among the four approaches corresponding to the GGM data and separately for the VAR data, are highlighted in bold font. The standard deviations corresponding to the reported metrics are presented in separate Tables in the Supplementary Materials.

Results for V=40	Network CP		Edge CP		MSE	Network CP		Edge CP		MSE
	sens	FP	sens	FP	MSE	sens	FP	sens	FP	MSE
	BPMM-PC					idPAC
GGM+Erdos-Renyi	0.91	7.31	0.50	1.12	0.1	1	2.75	0.92	1.08	0.09
GGM+Small-world	0.92	5.99	0.47	1.03	0.12	0.98	2.77	0.92	1.01	0.08
GGM+Scale-free	0.91	7.29	0.49	1.19	0.12	1	2.81	0.92	1.1	0.09
	SD+GFL					CCPD
GGM+Erdos-Renyi	0.3	3.13	0.09	2.97	0.29	0.92	2.15	0.31	4.1	0.16
GGM+Small-world	0.29	3.31	0.09	3.08	0.27	0.92	2.18	0.29	4.17	0.21
GGM+Scale-free	0.29	3.08	0.09	2.99	0.24	0.91	2.33	0.29	4.09	0.19
	BPMM-PC					idPAC
VAR+Erdos-Renyi	0.68	6.55	0.43	1.08	0.2	0.84	5.57	0.80	1.06	0.12
VAR+Small-world	0.66	5.97	0.47	1.14	0.19	0.77	5.54	0.74	1.12	0.09
VAR+Scale-free	0.59	5.51	0.39	1.02	0.17	0.78	5.29	0.73	1.06	0.09
	SD +GFL					CCPD
VAR+Erdos-Renyi	0.41	7.72	0.13	3.06	0.26	0.55	1.12	0.18	4.33	0.21
VAR+Small-world	0.56	6.29	0.14	2.98	0.19	0.64	1.36	0.17	3.47	0.23
VAR+Scale-free	0.42	6.99	0.17	3.13	0.22	0.58	1.27	0.19	3.29	0.2
Results for V=100	Network CP		Edge CP		MSE	Network CP		Edge CP		MSE
	sens	FP	sens	FP	MSE	sens	FP	sens	FP	MSE
	BPMM-PC					idPAC
GGM +Erdos-Renyi	0.92	4.77	0.51	1.31	0.11	1	2.31	0.83	1.16	0.09
GGM +Small-world	0.91	4.69	0.49	1.33	0.1	1	2.37	0.82	1.17	0.09
GGM +Scale-free	0.91	4.71	0.50	1.31	0.11	1	2.29	0.83	1.16	0.09
	SD +GFL					CCPD
GGM +Erdos-Renyi	0.3	3.13	0.09	2.97	0.29	0.9	1.12	0.29	4.6	0.18
GGM +Small-world	0.29	3.31	0.09	3.08	0.27	0.91	1.18	0.25	4.2	0.17
GGM +Scale-free	0.29	3.08	0.09	2.99	0.27	0.91	1.02	0.27	4.4	0.17
	BPMM-PC					idPAC
VAR+Erods-Renyi	0.66	5.97	0.51	1.07	0.14	0.82	5.88	0.81	1.04	0.11
VAR+Small-world	0.59	6.03	0.41	1.02	0.14	0.75	5.44	0.74	1.05	0.12
VAR+Scale-free	0.62	5.49	0.44	0.99	0.15	0.77	5.51	0.71	1.11	0.13
	SD +GFL					CCPD
VAR+Erdos-Renyi	0.37	8.03	0.1	3.14	0.15	0.55	1.09	0.17	3.75	0.22
VAR+Small-world	0.44	7.51	0.16	2.71	0.16	0.66	1.44	0.19	3.41	0.19
VAR+Scale-free	0.36	7.72	0.18	2.88	0.18	0.59	1.31	0.17	3.44	0.19

Table 5

Results under the dynamic precision matrix estimation approaches for network and edge-level connectivity change-point estimation (Edge CP) accuracy and network changepoint (Network CP) estimation accuracy for V = 40, 100. GGM and VAR correspond to data generated from Gaussian graphical models and vector autoregressive models respectively. Significantly improved metrics among the four approaches corresponding to the GGM data and separately for the VAR data, are highlighted in bold font. The standard deviations corresponding to the reported metrics are presented in the Supplementary Materials.

Results for V=40	Network CP		Edge CP		MSE	F1	Network CP		Edge CP		MSE	F1
	sens	FP	sens	FP	MSE	F1	sens	FP	sens	FP	MSE	F1
	BPMM-PR					idPMAC
GGM+Erdos-Renyi	0.85	6.99	0.32	1.04	0.1	0.79	1	5.2	0.79	0.89	0.08	0.88
GGM+Small-world	0.88	7.14	0.33	1.16	0.08	0.77	1	5.11	0.81	0.91	0.08	0.9
GGM+Scale-free	0.87	7.36	0.33	1.19	0.08	0.71	0.97	5.6	0.77	0.92	0.07	0.89
	DCR						SINGLE
GGM+Erdos-Renyi	0.22	16.15	0.41	9.39	0.27	0.59	0.35	6.49	0.1	2.84	0.08	0.71
GGM+Small-world	0.19	11.83	0.49	9.66	0.22	0.61	0.32	6.55	0.09	2.88	0.07	0.77
GGM+Scale-free	0.21	10.92	0.49	9.058	0.23	0.62	0.33	6.01	0.09	2.94	0.07	0.69
	BPMM-PR						idPMAC
VAR+Erdos-Renyi	0.66	4.45	0.29	1.16	0.10	0.77	0.79	4.81	0.68	1.22	0.09	0.81
VAR+Small-world	0.59	5.12	0.27	1.03	0.1	0.74	0.78	4.99	0.69	1.04	0.09	0.79
VAR+Scale-free	0.61	4.77	0.31	1.04	0.12	0.77	0.76	4.64	0.71	0.99	0.09	0.82
	DCR						SINGLE
VAR+Erdos-Renyi	0.22	9.83	0.4	3.35	0.24	0.64	0.42	7.35	0.13	3.11	0.27	0.66
VAR+Small-world	0.24	10.14	0.33	3.61	0.23	0.63	0.44	7.12	0.17	3.04	0.26	0.62
VAR+Scale-free	0.21	9.98	0.32	3.61	0.22	0.59	0.38	6.77	0.21	3.36	0.23	0.6
Results for V=100	Network CP		Edge CP		MSE	F1	Network CP		Edge CP		MSE	F1
	sens	FP	sens	FP	MSE	F1	sens	FP	sens	FP	MSE	F1
	BPMM-PR						idPMAC
GGM+Erdos-Renyi	0.92	6.83	0.28	1.09	0.08	0.83	0.97	5.1	0.82	0.89	0.08	0.89
GGM+Small-world	0.91	6.98	0.31	1.19	0.09	0.81	0.97	5.44	0.81	0.99	0.07	0.87
GGM+Scale-free	0.92	7.44	0.32	1.25	0.08	0.81	0.96	5.6	0.79	0.94	0.07	0.87
	DCR						SINGLE
GGM+Erdos-Renyi	0.33	16.14	0.41	9.39	0.22	0.63	0.38	6.77	0.12	2.97	0.08	0.69
GGM+Small-world	0.31	15.88	0.4	9.66	0.27	0.59	0.35	6.48	0.12	3.02	0.08	0.71
GGM+Scale-free	0.34	16.82	0.39	10.08	0.27	0.64	0.35	7.02	0.11	2.97	0.08	0.70
	BPMM-PR						idPMAC
VAR+Erdos-Renyi	0.73	4.41	0.29	1.18	0.14	0.77	0.88	4.22	0.63	1.09	0.13	0.82
VAR+Small-world	0.56	5.22	0.22	0.91	0.11	0.78	0.72	4.87	0.61	1.09	0.1	0.81
VAR+Scale-free	0.59	5.13	0.29	1.03	0.11	0.78	0.77	4.49	0.65	1.08	0.09	0.81
	DCR						SINGLE
VAR+Erdos-Renyi	0.23	9.92	0.43	3.19	0.16	0.64	0.42	7.41	0.14	3.11	0.11	0.71
VAR+Small-world	0.31	10.23	0.37	3.37	0.19	0.67	0.47	7.66	0.13	3.28	0.12	0.69
VAR+Scale-free	0.25	10.23	0.38	3.61	0.18	0.65	0.44	7.59	0.13	3.19	0.11	0.66

Table 4 reports the performance under pair-wise correlation based approaches, i.e. idPAC, BPMM-PC, SD, and CCPD. It is clear for the results that the proposed idPAC method has a near perfect sensitivity when data were generated under GGM, and a suitably high sensitivity under the VAR model, when estimating connectivity change points. The sensitivity for network and edge change point estimation, along with the MSE in estimating the pairwise correlations are significantly improved under idPAC compared to competing approaches in Table 4. The CCPD method is shown to have the lowest false positives when estimating the network level change points, but otherwise has poor sensitivity for change point estimation and high MSE, which is potentially due to the assumption of piecewise constant connectivity. The approach based on sliding window correlations has the poorest performance across all the reported metrics, which illustrates their drawback in estimating dynamic connectivity.

Table 5 reports the performance under precision matrix based approaches, i.e. idPMAC, BPMM-PR, SINGLE, and DCR. It is evident that the proposed idPMAC method has near-perfect or high sensitivity for detecting network level change points, corresponding to data generated under GGM and VAR models respectively. It also has a suitably high sensitivity for detecting edge level connectivity change points under both cases. Similarly, the MSE for edge strength estimation and the F-1 scores for network estimation accuracy are significantly improved under the proposed method in contrast to competing approaches. Figure 2 illustrates that the F-1 score over time under the proposed dynamic precision matrix method with covariates is almost always higher across almost all time scans compared to competing methods. Moreover the DCR and SINGLE method have the least impressive performance in terms of connectivity change point estimation, which also translates to poor dynamic network estimation (low F-1 scores).

Fig. 2.

F1-score over time for one single subject under the case of dynamic partial correlation method. The vertical green lines are the true change points. Red line represents the proposed method with dynamic partial correlation (idPMAC), the cyan line represents the covariate-naive version (BPMM-PM), the blue line represents DCR, and the pink line represents SINGLE method.

Our results clearly illustrate the advantages of the proposed methods over existing approaches that are not effective in leveraging information across samples. In addition, Tables 4–5 also illustrate the gains of incorporating covariate information under the proposed idPAC and idPMAC approaches over the covariate naive BPMM counterparts. It is interesting to note that the covariate naive BPMM still fares better than existing dynamic connectivity methods that fail to pool information across samples in a systematic manner. We also note that while the presence of false positive (FP) connectivity change points are expected due to the heterogeneity across samples, the proposed approaches provide desirable control of FP even while pooling information across samples with varying networks. In fact, the FP under the proposed method are lower than all competing methods except CCPD, whose performance is otherwise less impressive in terms of significantly lower sensitivity for change point detection, and inferior network estimation as reflected by poor MSE and F-1 scores.

When comparing the relative performance between idPAC and idPMAC methods, it is evident that the former has comparable or higher sensitivity but lower false positives in terms of estimating connectivity change points at the network level, when data are generated under a GGM. When data are generated under a VAR model, the idPAC method has higher sensitivity but also higher false positives compared to idPMAC, for estimating network connectivity change points. This is also true when estimating edge-level connectivity change points. In addition, since the idPMAC method estimates all edges simultaneously, the mean squared error for estimating edge strengths is often lower compared to the idPAC method. Moreover when the number of spurious covariates is increased, both these approaches experience a drop in performance (Fig. 3), as expected. It is of note that the number of false positive change points under the dynamic pairwise correlation approach increase minimally under the scale-free and small-world networks with an increase in the number of spurious covariates. However, this robust behavior was not replicated for network change points or other metrics of interest under the dynamic pairwise correlation method. In contrast, the recovery of the true clusters is shown to be more resilient under the dynamic precision matrix approach. This is evident from the top panels in Fig. 3 that show a slower increase in the clustering error under the idPMAC method.

Fig. 3.

Performance of dynamic pairwise correlation (columns 1 and 2) and dynamic precision matrix (columns 3 and 4) methods under different number of spurious covariates represented by the X-axis. Lines with different color represent different network structure: Green (Erdos Renyi), Red (Small World), Blue (Scale Free). The top row provides the information of clustering performance (Clustering Error and Variation of Information), the middle row demonstrates the performance of network level change points estimation (sensitivity and number of False Positive estimations), and the performance of edge level change point estimation was provided in the bottle row.

The proposed approach is clearly scalable to higher dimensional networks, with the dynamic pairwise correlation method being slightly faster than the dynamic precision matrix estimation approach. We found the computation time under the proposed approaches to be slightly slower than existing dynamic connectivity methods such as DCR, which results from additional computations related to clustering and due to incorporation of covariates. Table 6 presents the computation time in minutes for all approaches implemented via Matlab on a personal desktop computer (Alienware Desktop) that had Intel(R) Core i7–4930 processor with 32GB RAM (SINGLE was implemented via Python and hence the computation time is not reported). We note that the total computation time under BPMM is expected to increase with V, T, N, which is true for most dynamic connectivity approaches.

Table 6

Computation Time (in minutes) for simulation studies involving 300 time scans and 40 samples, under all approaches implemented via Matlab version R2017a.

Method	v=20	V=40	V=100
BPMM-PC	21	80	321
BPNN-PR	25	92	348
idPAC	27	102	402
idPMAC	31	114	416
SD+GFL	3	9	44
CCPD	70	315	844
DCR	18	90	297

Finally, we also conducted sensitivity analysis of the proposed approaches with respect to the model hyperparameters (see Tables 3–4 in Supplementary Materials). We found that moderate variations in the values of hyperparameters do not result in considerable changes in performance. Further, the clustering error is seen to be more sensitive to changes in hyperparameter values, while the metrics corresponding to dynamic network estimation seem more resilient to changes in hyperparameters, which suggests a degree of robustness for the estimated dynamic network with respect to the choice of hyperparameters. Although we do expect the performance of the methods to fluctuate to a greater degree for extreme choices of hyper-parameters, this is not of immediate concern to practitioners who use recommended values for the model hyper-parameters suggested in the manuscript.

Analysis of task fMRI data

Description of the study

We analyze a block task data involving a semantic verbal fluency at Veterans Affairs Center for Visual and Neurocognitive Rehabilitation, Atlanta. In a 12-week randomized controlled trail, 33 elderly individuals (aged 60–80, 11 males, 22 females) were assigned to two intervention groups: spin aerobic exercise group (14 participants) and the non-aerobic exercise control group (19 participants). During the intervention, individuals belonging to the aerobic spin group were required to do 20–45 min of spin aerobic exercise three times a week, led by a qualified instructor. For control group, participants were asked to do the same amount of non-aerobic exercise per week, such as group balance and light muscle toning exercise. A more detailed description of the data is available in Nocera et al. (2017).

For each participant, fMRI scans were conducted with 6 blocks of semantic verbal fluency (task) conditions with 8 scans, both pre- and post-intervention. The semantic verbal fluency task involved participants looking at different categories (e.g. “colors”) at the center of video screen and they were asked to generate and speak 8 different objects associated with that category (e.g. “blue”). After task block, a rest block with 3–5 TRs would appear and participants were required to read the word “rest” out loud. A total of 74 brain scans were acquired using a 3T Siemens Trio scanner with a whole-brain, 1-shot gradient EPI scan (240 mm FOV, 3.75 × 3.75 in-plane resolution, TR=5830ms, TA=1830ms, TE=25ms, flip angle (FA)=70). Analysis of Functional NeuroImages (AFNI) software and FMRIB Software Library (FSL) were used for pre-processing, as in Nocera et al. (2017). Slice-time corrections, linear trend removal, echo planar images alignment, and motion correction were performed as a part of the pre-processing pipeline. We used 18 brain regions for analysis that were shown to be differentially activated between the two intervention groups as described in Nocera et al. (2017). These regions are listed in Table 7 and comprise more regions in the right hemisphere due to decreased activity in that hemisphere in the spin group following the intervention, as compared to the control group. We note that since these regions corresponded to group differences due to spin exercise, they can not be described as “canonical” regions associated with semantic language function, which would also comprise some additional homologous regions in the left hemisphere. Since the purpose of the study was to investigate dynamic connectivity changes between brain regions due to the intervention, an analysis based on the selected 18 regions was undertaken instead of using canonical regions.

Table 7

Summary of brain regions used for analysis. R and L are abbreviations for right and left respectively.

ROI Number	Region name	Broadmann area	MNI coordinate
1	R Cerebullum 1	NA	(5,−62, −57)
2	R Inferior Temporal Gyrus	20	(41,−27, −30)
3	R Angular Gyrus	39	(44,−56,12)
4	R Middle Frontal Gyrus	10	(23,56,−6)
5	R Middle Temporal Gyrus 1	22	(53,−12, −9)
6	L Precuneus 1	7	(−9, −74,57)
7	L Cingulate Gyrus	NA	(−9, −33,39)
8	R Precuneus	7	(6,−80,48)
9	R Cerebellum 2	NA	(35,−53, −27)
10	R Middle Temporal Gyrus 2	21	(60,−45, −6)
11	R Inferior Frontal Gyrus/precentral gyrus	44	(59,9,9)
12	R Retrosplenial Area	30	(9,−47,18)
13	R Supramarginal Gyrus	40	(41,−36,33)
14	R Pars Triangularis/MFG	45	(47,47,−9)
15	L Precuneus 2	7	(−6, −71,45)
16	L Cuneus	19	(−15, −80,27)
17	L Superior Frontal Gyrus	6	(−17, −18,69)
18	R Middle Temporal Gyrus 3	22	(60,−36,0)

Analysis outline

We performed the analysis separately for the pre-intervention and post-intervention data, under both the dynamic pairwise correlations and dynamic precision matrix estimation methods. We used age and gender as covariates for the pre-intervention dataset, while also using the type of intervention (spin or non-aerobic control) as an additional covariate for the post-intervention analysis. Our analysis is designed to: (i) investigate the clustering behavior and inspect how these clusters differ with respect to demographics and the intervention type; (ii) investigate the cluster-level network differences using network summary measures; (iii) estimate the connectivity change points and examine how well they align with the changes dictated by the block task experiment; (iv) infer nodes and edges in the network with significantly different connectivity patterns between pre- and post-intervention.

Objective (i) enables us to characterize homogeneous dynamic connectivity patterns corresponding to clusters of samples in terms of their demographic and clinical characteristics; aim (ii) will be instrumental in interpreting the cluster-level network differences that will shed light on network variations across transient network states; aim (iii) will provide insights regarding the effectiveness of the proposed approaches in terms of recovering connectivity jumps where these changes are influenced by, but often not fully aligned with, the changes in the block task experimental design (Hindriks et al., 2016; Kundu et al., 2018); and aim (iv) will inform investigators regarding dynamic connectivity differences that are associated with the type of intervention. For aim (ii), we were only able to report results under dynamic precision matrix estimation, since a graph theoretic framework is necessary to compute the network summary measures, which may not be feasible under a pairwise correlation analysis.

Cluster analysis:

As seen from Table 8, the analysis under both idPAC and idPMAC methods yielded 5 clusters consolidated over all time scans (using the K-means algorithm described in Section 2.2, although the size of the clusters were more equitable under the idPAC method. The pre-intervention analysis yielded clusters that were largely homogeneous with respect to gender. These clusters were also reasonably well-separated with respect to age under the idPAC analysis, whereas the age of the participants within clusters were more diverse under the idPMAC analysis. The post-intervention analysis yielded more heterogeneous clusters with respect to both age and gender, with only one cluster comprising all males under both the idPAC and idPMAC analyses. This suggests a realignment of the dynamic connectivity after the intervention is administered, such that individuals with similar genders and age-groups have synchronous dynamic connectivity patterns pre-intervention as identified via subgroups, but the subgroups and their composition with respect to age and gender change post-intervention. Our post-intervention analysis also suggests that the variability across clusters under the idPAC method can be largely explained via the intervention type.

Table 8

Results for analysis of block task fMRI experiments. Size refers to the number of participants in each cluster, ‘CP(Task-Rest)’ and ‘CP(Rest-Task)’ denotes the cluster level connectivity change points that were detected within +/− 2 time scan of the change in experimental design from task to fixation, and from fixation to task, respectively. ‘Spin’ refers to the percentage of individuals assigned the Spin intervention belonging to each cluster.

Method	idPAC					idPMAC
Cluster index	1	2	3	4	5	1	2	3	4	5
Cluster features	Pre-intervention					Pre-intervention
Size	8	6	8	7	4	3	5	17	6	2
% of females	0	100	0	14	100	0	100	0	100	0
Age (mean)	72.2	65.8	64.7	76.7	67.7	71.7	69	70.4	66.8	67
Age(range)	69–73	60–72	60–68	74–80	66–69	63–78	62–80	60–80	60–72	66–68
CP(Task-Rest)	6	3	4	4	4	4	5	5	4	3
CP(Rest-Task)	3	5	2	3	4	4	4	2	4	3
	Post-intervention					Post-intervention
Size	8	4	7	11	3	3	4	9	11	6
% of females	63	75	0	18	33	67	100	0	9	67
Age (mean)	67.3	65	65.1	74.5.	73.7	73.7	69.3	68.6	73	62.7
Age(range)	62–70	60–71	60–68	71–80	68–78	67–80	68–72	63–78	68–80	60–66
CP(Task-Rest)	5	6	4	3	6	3	3	5	5	5
CP(Rest-Task)	3	5	2	2	4	2	5	2	4	2
Spin(%)	0	100	100	0	100	33	0	100	9	50

Connectivity change point estimation:

Table 8 illustrates the cluster level connectivity change point estimation. We observed that under both the idPAC and idPMAC methods, the estimated change points were consistent with 4 or more (out of 6) changes in experimental design when transitioning from task to rest, except one cluster where 3 of the connectivity change points aligned with the experimental design. These patterns were consistent in both the pre- and post-intervention analysis; however the number of connectivity change points that were strongly aligned with changes in the experimental design were (on average) greater in the post-intervention analysis compared to the pre-intervention analysis. This suggests a learning effect of the task that was reflected in terms of higher concordance between the connectivity change points and the experimental design post-intervention. On the other hand, the cluster-level estimation of change points when transitioning from fixation to task was (on average) less aligned with the experimental design compared to the change points when transitioning from task to fixation, as seen in Table 8. This is somewhat expected since there were only 3–5 time scans in each fixation block, which made it extremely challenging to detect connectivity changes when transitioning from fixation to task. However, the proposed approach was still able to detect at least two, and often 3 or more connectivity change points (out of 6) aligned with the experimental design that suggests a reasonable concordance between connectivity jumps and experimental transitions from fixation to task.

In contrast, the CCPD approach detected at most one or two connectivity change points, while the DCR method was not able to detect connectivity change points at all, which makes these results appear biologically impractical given the nature of the block task experiment. Although the changes in connectivity are not expected to be fully aligned with changes in the experimental design (Hindriks et al., 2016), one expects a certain degree of synchronicity between the two. Our results indicate that this is not captured at all via existing change point methods especially when there are rapidly occurring transitions in the experimental design, which highlights their limitations. Hence, our analysis clearly illustrates the advantages of pooling information across heterogeneous samples and incorporating covariate knowledge via a mixture modeling framework, which is simply not possible using existing approaches that rely on information from single subjects as in DCR, or that use empirical methods to pool information across individuals as in CCPD.

Cluster level network differences:

In order to investigate the differences between the networks corresponding to the different clusters, we examined variations in dynamic network metrics that capture modes of information transmission in the brain. These network metrics include the characteristic path length (CPL) that measures the length of connections between nodes, and the mean clustering coeffcient (MCC) that measures the clustering tendency averaged over all network nodes. Using permutation testing, we examined p-values to evaluate which pairs of clusters exhibited significantly different network summary measures. None of the clusters had significantly different CPL values in the pre-intervention analysis, but several pairs of clusters exhibited significant CPL differences post-intervention. The CPL differences were particularly pronounced between the first and remaining clusters, as well as the last and remaining clusters in the post-intervention analysis. These two clusters also demonstrated the highest within cluster variability in CPL values amongst all clusters. Moreover, the number of pairs of clusters with significantly different MCC values increased from the pre-intervention to post-intervention analysis, with 8 out of 10 pairs of post-intervention clusters reporting significantly different MCC values compared to at least one other cluster. Hence, our results suggest greater variability in network organization between clusters in the post-intervention analysis compared to pre-intervention, which potentially reflects greater network heterogeneity after the 12 week intervention was administered.

Network differences pre- and post-intervention:

We applied paired ttest with multiplicity adjustment (using Bonferroni correction) in order to infer which edges were significantly different between pre- and post-intervention at 5% level of significance, along with identifying which network nodes contained the greatest number of differential edges. Since the magnitude of the pairwise correlations and the corresponding edge strength differences were higher, we discovered higher number of edges with differential edge strengths under the idPAC analysis. For both the idPAC and idPMAC methods, the bulk of the pre- vs post-intervention connectivity differences were concentrated in individuals in the spin group exclusively that were not present in the control group. We obtained 57 significantly different edges under the idPAC analysis, and 38 significantly different edges under the idPMAC analysis, which were exclusive to the spin group see Fig. 4. In contrast, the number of significantly different edges between the pre- and post-intervention networks under the idPAC analysis were 20 corresponding to both the spin and control groups, and 7 corresponding to the control group only. Moreover the idPMAC analysis did not produce any significant edge level differences between the pre- and post-intervention networks corresponding to both the intervention groups as well as for the control group only. Our results suggest a considerably strong realignment in dynamic connectivity after the 12-week intervention that were exclusive to the spin group, compared to negligible changes in the control group.

Fig. 4.

Circle plots for the edges that are significantly different pre- and post-intervention in spin group but not in the control group. The top and bottom panel correspond to the results under dynamic pairwise correlation and dynamic precision matrix estimation incorporating covariates, respectively. Red and blue lines correspond to lower or higher edge strengths in the pre-intervention network compared to post-intervention. RC1 and RC2 refer to the two brain regions in the right cerebellum; RMTG1-RMTG3 refer to the three brain regions in the right middle temporal gyrus; and LP1-LP2 refer to the two regions in the left precuneus. The MNI coordinates for these regions are provided in the Figure legend.

The changes between the pre- vs post intervention networks that occurred exclusively in the spin group under idPAC analysis were concentrated in the following brain regions: Right Angular Gyrus(8 edges), Left Precuneus(10 edges), Right Cerebellum(9 edges), Right Middle Temporal Gyrus(11 edges), and Right Middle Temporal Gyrus(8 edges). Similarly the following brain regions had the highest number of differential edges pre- vs post-intervention under the idPMAC analysis: Right Middle Frontal Gyrus(16 edges), Right Cerebellum(6 edges), Right Pars Triangularis/MFG(8 edges), and Right Middle Temporal Gyrus(7 edges). Two nodes, Right Cerebellum and Right Middle Temporal Gyrus had a large number of significantly differential edges under both idPAC and idPMAC analyses, while the right middle frontal gyrus had, by far, the largest number of differential edges (16) under the dynamic precision matrix analysis. In addition, we also observe that more nodes in right hemisphere of the brain have significantly differential connectivity, which is to be expected since the majority of the 18 brain regions being investigated lie in the right hemisphere.

The large number of differential connections with respect to the right cerebellum is believed to be attributable to the generation of internal models or context specific properties of an object (Moberget et al., 2014), and preferential activation during a semantic challenge (D’Mello et al., 2017). The connectivity between the right cerebellum and inferior frontal regions has been noted in earlier studies (Balsters et al., 2013), with the inferior frontal regions being responsible for ordering language and codifying the motor output for syntax (Balsters et al., 2013). Moreover, the differential connectivity in the right middle temporal gyrus is along the lines of earlier findings that illustrated the role of the left temporal gyrus as a hub for integration of sensory input into a transformation to semantic forms (Davey et al., 2016), and the corresponding connectivity differences in the right middle temporal gyrus may be attributable to a shift in laterality of involvement (Lacombe et al., 2015) due to aging. Finally, the large number of differential edges corresponding to the right middle frontal gyrus is potentially associated with semantic priming in older adults (Laufer et al., 2011). Given that this region is associated with executive function (Jolles et al., 2013; Wang et al., 2019) and is well characterized as being involved in working memory tasks, it is likely for connectivity differences to be focused on this region since the semantic task requires a continuous reference to working memory.

Discussion

In this article, we developed a novel approach that accurately estimates a population of subject-level dynamic networks by pooling information across multiple subjects in an unsupervised manner under a mixture modeling framework using covariates. The proposed approach, which is one of the first of its kind in dynamic connectivity literature, results in significant gains in dynamic network estimation accuracy, as illustrated via extensive numerical studies. The gains under the proposed method become particularly appealing compared to existing approaches in the presence of rapid transitions in connectivity as evident from our fMRI block task analysis. The proposed approach works best in fMRI task experiments involving a group of heterogeneous individuals executing the same task protocols, and in the presence of a carefully chosen set of covariates that are related to the dynamic network.

We also illustrate the robust performance of the proposed approach in the presence of a limited number of covariates that are not related to changes in connectivity, although the performance deteriorates as the number of spurious covariates increase. In the presence of a large number of features that may not be necessarily related to dynamic connectivity, one can perform a screening step to exclude unimportant predictors from the analysis. This step will involve examining the associations between each covariate and the dynamic connectivity estimates obtained from the covariate naive BPMM approach, and subsequently only retaining the covariates with significant associations for analysis using the full model. This approach is expected to work well as long as the screening step does not exclude any important covariates and manages to largely filter out spurious covariates that are unrelated to the network. In future work, we plan to extend the proposed approach to incorporate feature selection that automatically identifies significant covariates that are related to the dynamic networks, and down-weights the contribution of unimportant covariates using Bayesian shrinkage priors. We note that although our analysis included covariates that do not vary with time, the proposed BPMM approach can be easily generalized to include time-varying covariates that are collected in-scanner (such as behavioral performance), when required.

Moreover for networks with higher densities or those with large variability in the values of the non-zero elements in the precision matrix, it is possible that the proposed approach may result in sub-optimal performance under a small number of mixture components. This is due to the symmetry constraint, which may potentially impose some restrictions on the clustering of the precision matrix elements when the number of clusters is small, and hence lead to inaccurate estimates. In such cases, a larger number of clusters would be required to fit the data well. Another avenue to tackle such potential restrictions is to generalize the proposed BPMM approach in equation (3) so as to specify independent mixture priors on each element of the upper triangular precision matrix subject to the constraint which will impose the same marginal distributions on elements (i, j) and (j, i) in the precision matrix. We plan to explore such generalizations in future work.

In addition to identifying important connectivity changes, during the fMRI block task experiment, our analysis conclusively established major changes between the pre- and post-intervention networks that were exclusive to the spin group. We note that existing literature has established the role of cardiovascular fitness in regulating aging related declines in both language and motor control (McGregor et al. (2011), 2013). However, much less is known about the effect of exercise intervention on dynamic connectivity, particularly in older adults. Because connectivity is a fundamental aspect of neuronal communication required for high-level cognitive processes, it is important to understand the potential impact of aging and/or aerobic exercise interventions in aging on changes in brain connectivity.

Further, our analysis also discovered subgroups of individuals with homologous dynamic connectivity, where the heterogeneity within these subgroups with respect to intervention was higher under the idPMAC method compared to the idPAC analysis. This indicates that dynamic pairwise correlations were more accurate in classifying participants in terms of the intervention administered. It is important to note that the separation of clusters with respect to intervention reflects the distinct patterns of dynamic connectivity between the 18 brain regions specified in our study that are known to be differentially activated in spin and control groups (Nocera et al., 2017). However, if additional regions are included that may not be necessarily associated with intervention type, it is entirely possible to obtain more heterogeneous clusters that have a more equitable composition with respect to intervention group. This is due to the presence of network edges between regions that are not necessarily associated with intervention and hence behave similarly in both the spin and control groups. Future work will focus on a more general analysis involving a larger number of cannonical regions known to be associated with the semantic language function.