A new model for simultaneous dimensionality reduction and time-varying functional connectivity estimation.

An important question in neuroscience is whether or not we can interpret spontaneous variations in the pattern of correlation between brain areas, which we refer to as functional connectivity or FC, as an index of dynamic neuronal communication in fMRI. That is, can we measure time-varying FC reliably? And, if so, can FC reflect information transfer between brain regions at relatively fast-time scales? Answering these questions in practice requires dealing with the statistical challenge of having high-dimensional data and a comparatively lower number of time points or volumes. A common strategy is to use PCA to reduce the dimensionality of the data, and then apply some model, such as the hidden Markov model (HMM) or a mixture model of Gaussian distributions, to find a set of distinct FC patterns or states. The distinct spatial properties of these FC states together with the time-resolved switching between them offer a flexible description of time-varying FC. In this work, I show that in this context PCA can suffer from systematic biases and loss of sensitivity for the purposes of finding time-varying FC. To get around these issues, I propose a novel variety of the HMM, named HMM-PCA, where the states are themselves PCA decompositions. Since PCA is based on the data covariance, the state-specific PCA decompositions reflect distinct patterns of FC. I show, theoretically and empirically, that fusing dimensionality reduction and time-varying FC estimation in one single step can avoid these problems and outperform alternative approaches, facilitating the quantification of transient communication in the brain.

This is a PLOS Computational Biology Methods paper.

Introduction

When we image the brain of passive subjects with fMRI, the measured signals exhibit strong correlations even between areas that are far apart in the brain [1, 2]. These patterns of resting-state correlation, referred to as functional connectivity (FC), are interpreted as a sign that these regions are somehow engaged together in relation to one or more brain processes [3]. It is now widely recognised that FC holds important relations to mental and clinical phenotypes, is reliably subject-specific (i.e. reproducible across scanning sessions), and is also hereditary [4-6]. However, the mere existence and interpretability of within-session modulations in FC is, at least in fMRI, more controversial; see [7, 8] for arguments in both directions. An important reason for this dispute is the scarcity of time samples and the very high dimensionality of the data. In this context, quantifying modulations of FC within session is a challenging statistical problem because these changes (if they exist) are by definition spontaneous and have no obvious behavioural reference [9].

One possible strategy to quantify time-varying FC is the use of sliding-windows, where some measure of FC such as Pearson’s correlation is computed across regions for each window in the data, typically followed by the application of a clustering algorithm to extract patterns across windows [10]. Although attractive because of its simplicity, this method suffers from an important problem of statistical variability in the estimation, such that disentangling actual changes in the data from fluctuations caused by statistical noise is not straightforward [11]. Alternatives to Pearson’s correlation include phase coherence [12, 13], and the angle [14] or covariance [15] between the signal gradients; these have been shown to necessitate shorter windows but might however be more sensitive to high-frequency noise artefacts.

Dispensing with the use of windows, methods that boost the statistical power by using the entire data set in the estimation are sometimes preferred. One such method is the Hidden Markov model (HMM), which assumes that the data can be reasonably modelled using a discrete number of FC states with Markovian dynamics [16]. For example, if an HMM with twelve states was trained on 820 subjects from the Human Connectome Project (HCP) data set [17], each state would be on average estimated on 68.3h of data; compared to a typical 1min window, the statistical noise in this estimation is very small. In the HMM, the Markovian property means that the model accounts for the state dynamics using a probability matrix that encodes the probability of transitioning between each pair of states—but without modelling the previous history of state activations [18, 19]. A straightforward variant of the HMM is to model each state as a Gaussian distribution where the mean is pinned to zero in order to prioritise changes in covariance [8]. An alternative is the mixture model of Gaussian distributions [19], which has no transition probability matrix and therefore ignores the temporal structure of the data.

Unfortunately, neither the HMM nor the mixture model of Gaussian distributions are easily applicable when the data dimensionality (the number of voxels) is too high in relation with the number of time points (volumes). The two most common approaches for reducing dimensionality in this context are using an anatomical parcellation [20] and applying independent component analysis (ICA) [3]. These produce a number n of regions or components (typically a hundred or more), so that an FC matrix has n(n − 1)/2 different parameters. Often, this is still high enough for the HMM inference (or mixture model inference) to overfit and produce degenerate solutions in many data sets. For this reason, PCA is usually carried out on the ICA-derived or parcellated time series so that the HMM (or alternative method) is run on an even lower-dimensional space [21].

Whereas this two-step method works reasonably well in practice for many (but not all) fMRI data sets, in the sense that it captures within-session FC modulations at least to some extent, having PCA and the HMM estimation as two separate steps is suboptimal. This is because the PCA step specialises in maximising explained variance, and is not designed to the specific goal of quantifying within-session FC modulations. Furthermore, the use of PCA can unknowingly introduce important biases on the estimation. In this paper, I discuss these issues in detail and propose an alternative model that bypasses these problems: an HMM where each state is a PCA decomposition. I refer this model to as HMM-PCA. Critically, because the computation of PCA is based on the data covariance [22], the state-specific PCA decompositions reflect distinct patterns of FC, effectively fusing dimensionality reduction and time-varying FC estimation in one single step. Fig 1 presents a graphical illustration of both the HMM with Gaussian states over principal components (HMM-Gaussian, top), and the HMM-PCA (bottom).

Fig 1

Two different approaches for the estimation of time-varying FC on high-dimensional fMRI data.

A. PCA is first used as a dimensionality reduction step, blindly to the purpose of estimating time-varying FC; then, some state-based model (like the hidden Markov model) is run on the first principal components (PC). B. The HMM-PCA approach, where each state is a different PCA decomposition, is run directly on the high-dimensional data; given that the computation of PCA is based on the data covariance, different PCA decompositions capture different patterns of FC. See S1 Fig for representations in the form of graphical models.

Materials and methods

The problem of estimating time-varying FC in high dimensions

Let be the multivariate source signal at volume (time point) t = 1…T, so that denotes the data concatenated for all sessions and subjects –although this can also be applied at the single-subject level provided that we have a sufficient number of volumes. Here, n corresponds to anatomical parcels or ICA components, referred to generically as regions. A standard estimation of functional connectivity (FC) for subject j is an n × n matrix C containing the Pearson’s correlation coefficient for each pair of regions. Formally, the question at hand can be posed as: can we find differences between matrices C(t1) and C(t2), defined as instantaneous FC matrices at time points t1 and t2, for at least one pair of time points t1 and t2 belonging to the same scanning session?

One way to approach this problem is the use of the Hidden Markov Model (HMM), which describes the data in terms of a finite number K of states that activate or deactivate throughout the scanning time. The state time courses, reflecting these activations, are in the form of probabilities P(x = 1|) = γ, where x = 1 means that the state k is active at time point t, and γ is the estimated posterior probability of the event x = 1, given the data. The HMM is a generic model where the states can be described using any family of probability distributions. Within this general framework, we can define the states as covariance matrices Σ. In this case, the states are characterised as Wishart distributions or, equivalently, as zero-mean Gaussian distributions [8]. This means that, when the k-th state is active, the data is considered to be generated according to the distribution I refer to this approach as HMM-Gaussian. Alternatively, the mixture model of Gaussian distributions dispenses with the transition probability matrix, thus ignoring the temporal structure of the data and treating the time points (volumes) as independently distributed and exchangeable [19]. In this case, the states are also defined as Wishart or zero-mean Gaussian distributions, but the transition probability between consecutive time points is not modelled. I refer to this model as Mix-Gaussian. In either case, an instantaneous FC matrix at time point t can be defined as a linear combination of the state covariance matrices Σ using the state probabilities γ as weights.

Both the HMM and the mixture model necessitate the number of states K to be pre-specified. Although there are extensions of the models where the number of states is directly inferred from the data (referred to as infinite HMM [23] or infinite mixture model [24], respectively), their parameter estimation is computationally demanding when we have large amounts of data and, in fMRI, they often yield similar results to the “finite” models [25]. Therefore, here we prespecify K and submit the data to an inference algorithm that will return the state probabilities γ, the state parameters (Σ), and the transition probability matrix and γ.

Because n is often large in comparison to T, PCA is typically used as an intermediate dimensionality reduction step. This way, HMM-Gaussian (or Mix-Gaussian) is estimated on = d, where represents a PCA decomposition and p is the number of principal components (PCs). The estimated FC matrices are thus low-dimensional, . Note that, across the entire PCA-reduced data set, , the p PCs are by construction orthogonal (i.e. the correlation between the columns of is zero). However, there might be periods in the data during which the time series are temporarily not orthogonal, meaning that the time series are temporarily correlated (or negatively correlated) across brain regions. Any transient departure from (the time-averaged) orthogonality will be encoded by the HMM parameters γ and Σ, and can be considered an FC modulation.

HMM-PCA: A new model for estimating time-varying FC

I now introduce the HMM-PCA mathematically. Various of the elements of this model are analogous to the probabilistic mixture model of PCA analysers introduced by Tipping and Bishop [26] (here referred to as Mix-PCA), which, similarly to Mix-Gaussian, does not account for the temporal structure of the fMRI data. The model is also closely related to that of Alvárez and Henao [27]. In brief, the main read-outs of HMM-PCA and Mix-PCA are (i) a set of K states, each characterised by a PCA decomposition; and (ii) the corresponding state time courses γ, which encode the probability of each state k to be active at each time point t. In the case of HMM-PCA, a matrix with transition probabilities between states is also estimated.

Both Mix-PCA and the HMM-PCA are based on probabilistic PCA [28], which formulates classic PCA within a probabilistic framework. Probabilistic PCA assumes the following distribution: where the covariance is given by = ′ + σ2, and is the identity matrix. The density function is therefore where is assumed to be Gaussian distributed, and || ⋅ || is the Euclidean norm. Since for our purposes we are not interested in modelling transient changes in amplitude, and we wish to concentrate the model’s explanatory power on FC modulations as much as possible, I introduce the modification = 0 to the model proposed by [26], therefore not modelling changes in amplitude explicitly.

To model the data, the Mix-PCA model uses, as states, K different PCA projections, , and their corresponding noise variance estimations , which are estimated from the data together with the state occupancies. The state covariance matrices are denoted as . The prior probability for each state to have generated unseen data is given by π, which also needs to be estimated. Therefore a second set of latent variables is required, such that x = 1 if the k-th component (or state) is responsible for having generated the observed data at time point t, and x = 0 otherwise. I define = [1…], and refer to the posterior probabilities P(x = 1|, π) = γ as the state time courses. Note that is conditionally independent of all data time points except , therefore ignoring the data temporal structure. See S1 Fig for a graphical representation of the model conditional independences of HMM-PCA and the other variants. The posterior probabilities γ can be estimated, for example, using the expectation-maximisation (EM) algorithm [26, 28], or which I provide some details below.

In the HMM-PCA case, I instead have the state latent variables modelled as an order-1 Markovian process, with prior probabilities . Here, Θ is the transition probability matrix, which models the average probabilities of transitioning from one state to another—and must be estimated as well. This way, the posterior probabilities (namely, the state time courses) are defined as P(x = 1|x(, , Θ) = γ.

To use the EM algorithm to solve the HMM-PCA problem (and assuming for simplicity of notation that we have one single, continuous time series), the expected log-likelihood can be formulated as where 〈⋅〉 denotes expectation and (, σ2) comprise these variables for k = 1…K, and where

Similarly to the Mix-PCA model, the EM updates for and are coupled. Using the same estimation of the intermediate variable for both, the new parameter estimates for these variables become: where is the state-specific (sample) covariance matrix in the original space, such that the aggregated sample covariance matrix can be expressed as a weighted average of . For the full EM algorithm, it only remains the estimation of the state activation probabilities γ, the transition probability matrix Θ, and the initial probabilities for the first time point of the scans.

As is standard practice, the γ parameters are estimated using the forward-backward equations, given the likelihood for each HMM state (here, PCA decomposition) at t [18]. The update rules for Θ and are also equivalent to any other HMM (given the expected log-likelihood in 5), and can be found elsewhere [19, 29].

Results

Limitations of the two-step approaches

Previous work has shown that the HMM, when ran on PCA time series, can produce useful representations of the data [21]. However, this approach suffers from two limitations: (i) a loss of statistical efficiency in detecting time-varying FC when there is time-varying-FC-relevant information in the discarded PCs, and (ii) a bias towards the lower-order PCs that were included in the model (that is, those explaining less variance). These limitations, which I discuss next, equally apply to other probabilistic models and clustering methods such as Mix-Gaussian when applied on PCA time series. Note that the following discussion merely states two mathematical facts intrinsic to the use of PCA. How much these biases actually affect real data will likely depend on many factors including data acquisition, preprocessing, and experimental paradigm.

Loss of sensitivity

I first show the loss of sensitivity in detecting time-varying FC when the discarded PCA components contain time-varying-FC-related information. Given that within-session modulations in FC are bound to be subtle [7, 8], it is quite possible that such modulations will indeed occur in lower-order PCs. Since this is not straightforward to show without access to the ground truth, I used a simulation where, as it happens in fMRI, the temporal modulations of covariance are not very large. This simulation illustrates that HMM-PCA may outperforms competing models if there is time-varying FC in lower-order PCs. Secondarily, these simulations stress the importance of acknowledging the temporal nature of the data, which is ignored by Mix-PCA.

I generated data from two different simulation schemes. In both cases, the data were generated from a low-dimensional space and projected to the dimension of the observed data (n = 10), where some observational 10-dimensional white noise was added. (That is, the data is low-rank up to the observational noise). The nature of this projection varies according to two different states, which transitions are organised as an order-1 Markovian process with transition probabilities, if k1 = k2 and otherwise; that is, the probability of remaining in the same state is 25 times higher than that of transitioning. I sampled 10 sessions of 1000 data points each. Note that I generated the data from a Markovian process to make it simple, but the fact that the generative process aligns with the Markovian assumption of the HMM does not result in any loss of generalisation for the main point of the simulations.

In Scenario 1, I separately sampled each of the latent dimensions from a zero-mean Gaussian distribution, so that these are approximately orthogonal. I denote the dimension of the latent space as p0, and the generated latent data as 0. Two different separate cases are considered: in the first, I have p0 = 2, where I sampled two latent variables with standard deviations 2.0 and 1.5; in the second, I have p0 = 3, with respective standard deviations of 2.0, 1.5 and 1.0. I then sampled three random (standard Gaussian-distributed) projection matrices , and generated the observed data as = 0(A + A1) + ϵ when state 1 is active, and = Y0(A + A2) + ϵ when state 2 is active. Observational noise ϵ is set to have zero mean and standard deviation equal to 0.001. The ground-truth state time courses (telling when each state is active) were generated using Markov chain Monte Carlo sampling as mentioned above.

In Scenario 2, I based the sampling on actual fMRI data from the HCP (see above). In particular, I randomly chose n = 10 regions from the data and computed the data covariance matrix , which I then eigendecomposed into its principal components: = ′, where is a diagonal matrix of eigenvalues and contains the corresponding eigenvectors. I used these to create two covariance matrices 1 and 2, which were then used to sample the data for each of the two states given a multivariate Gaussian distribution with zero mean. As with the other scenario, I considered two cases: p0 = 2 and p0 = 3. In either case, the first eigenvector of both 1 and 2 was set to be the first eigenvector of , thus corresponding to the time-invariant component of FC. In the p0 = 2 case, the second eigenvector of 1 and 2 were each assigned a different permutation of the second eigenvector of ; the rest of eigenvectors were set to zero. In the p0 = 3 case, the second and third eigenvectors of 1 and 2 were assigned different permutations of the second and third eigenvector of , and the rest of eigenvectors were set to zero. These latter non-zero eigenvectors, therefore, represent the time-varying components of the FC. Again, additive observational noise ϵ was set to have zero mean and standard deviation equal to 0.001. The transition probability matrix was designed as before.

I repeated the simulations 50 times, and estimated HMM-PCA, Mix-PCA and HMM-Gaussian models, which were all set to use p = 2 components, and, for simplicity, the right number of states K = 2. Therefore, there is loss of information only in the p0 = 3 cases. For illustration, Fig 2A shows one specific instance for Simulation 1 and p = p0 = 2, where the HMM-Gaussian approach misses some of the swiftest state changes, and Mix-PCA, on the other hand, appears to be noisier due to the lack of consideration for the temporal structure of the data. Fig 2B shows the complete results. Here, accuracy corresponds to the proportion of time where the corrected state was guessed. Each dot corresponds to one instance of the simulations, representing how the HMM-PCA compares to the Mix-PCA (blue) or to the HMM-Gaussian (red). Therefore, the points that lie to the left of the diagonal line correspond to simulations where the HMM-PCA performed better than the other models, and the points that lie to the right represent simulations where the HMM-PCA performed worse. Given that K = 2, and because the order of the states is non-identifiable, this measure of accuracy ranges between 0.5 and 1.0. Models with accuracy of around 0.5 correspond to degenerate solutions, where one of the two states was obliterated by the inference. A summary of the average accuracy is shown in the text boxes. As expected, HMM-PCA has almost perfect accuracy for the p0 = 2 cases, and deteriorates to a moderate extent when p0 = 3. Most importantly, HMM-PCA also outperforms HMM-Gaussian in most simulations, confirming that approaching the problem in one single step leads to superior solutions when there are time-varying FC information in low-order PCs; i.e. when time-varying FC modulations are subtle, as it is the case in most fMRI data sets. HMM-PCA also performs better than the Mix-PCA, highlighting the importance of accounting for the temporal structure of the data. (Note that in real data the temporal structure is stronger than in these simulations, so this difference will be even larger).

Fig 2

HMM-PCA outperforms the HMM-Gaussian and Mix-PCA approaches on synthetic data.

A. Example of how the different models recover the ground-truth state time courses. B. Comparative accuracy between HMM-PCA (Y-axis), HMM-Gaussian (X-axis, red) and Mix-PCA (X-axis, blue). Each dot represents one repetition of the simulations. Accuracy is measured in terms of how well each method recovered the ground-truth state time courses. Permutation-based statistical testing revealed that HMM-PCA was always significantly more accurate than the other approaches (p<0.001).

Bias towards low-order PCA components

The previous section discussed the suboptimality of the HMM-Gaussian solutions when there is time-varying FC information in lower-order PCs. I now discussed an intrinsic bias that will manifest regardless of how the time-varying FC information distributes across the PCs, and that will occur above and beyond the natural loss of information of PCA. This is that the application of PCA systematically alters the HMM or mixture model inference with regard to the original estimation (i.e. the one obtained without using PCA), biasing it towards low-order PCs and introducing a factor of arbitrariness in the inference. Critically, this issue occurs even when we keep all PCs and retain 100% of the variance.

When states are described as Gaussian distributions, the HMM (or mixture model) inference is scale-invariant. That is, as far as we use the same random seed in the initialisation of the inference, the estimation of the state time courses will not be affected if we multiply the time series of any given region by any random scalar. Mathematically, this can be expressed as where G(; ϒ) represents the HMM inference process given some specification of hyperparameters ϒ (e.g. the number of states). That is, the inference remains unaltered after rescaling the regions’ time series by any vector , regardless of the specific values of such vector. Intuitively, the reason is because the state-specific covariance matrices Σ can adjust their diagonal (i.e. their variance) to compensate for this global scaling with no effects in the inference. (This is as far as the prior distribution of the covariance matrices acknowledges this scaling; if not, there might be small changes in the inference but rarely substantial provided that we have enough data).

However, the HMM (or mixture model) inference is not rotation-invariant, and, in particular, it is not invariant to a PCA rotation:

When we apply PCA on the data, = , the columns of are ordered according to its variance, such that the first column of (the first eigenvector) has the highest variance and the last column has the lowest variance. Because of the scale-invariance property of the HMM inference, however, these variances will be ignored. Intuitively, this means that the low-order PCs (which explain less variance in the original data) are given in principle the same weight in the inference than the high-order PCs (which explain more variance in the original data). In different words, a PCA decomposition has information in both the eigenvectors and the eigenvalues, but the eigenvalues are ignored in the HMM estimation, therefore changing the HMM estimation with regard to what would be obtained in the original data. In practice, that results in a distortion of the estimates with respect to the original data, which will become more drastic as we include more and more low-order PCs.

For example, let us consider one given subject from the HCP data set [17], eight randomly-chosen brain parcels from the (100 regions) Schaefer parcellation [30], and a fixed initialisation of the algorithm (i.e. initialising the inference with exactly the same starting state time courses, so there is no randomness anywhere throughout the inference). Fig 3A shows the estimated state time courses for the original data (top), the state time courses estimated after scaling each channel randomly (middle), and the state time courses obtained from a PCA decomposition where we kept all components so that there is no loss of information (bottom). As observed, PCA-rotating the data changes the estimation, whereas scaling does not.

Fig 3

PCA introduces an estimation bias on the HMM or mixture model.

A. The HMM inference is scale-invariant, but not PCA-rotation-invariant. State time courses produced by the HMM inference for one HCP subject on the original parcellation space (top), after applying a random scaling of the data (middle), and after PCA rotation with no loss of variance (bottom). Each colour represents a different state, so that the coloured areas indicate the probability of activation for the states across the session. The similarity between the different runs, expressed as Pearson’s correlation coefficients, are expressed on the right. B. The extent of this bias (PCA distortion) is logarithmically related to how concentrated is the variance on the first PCs (PCA concentration); that is, the more correlated are the regions on the original data, the stronger will be the bias introduced by PCA. C. Random manipulations of the data eigenvalues are correctly reflected as changes in the HMM-PCA estimations; HMM-Gaussian is not able to capture the changes.

Furthermore, this bias is related with the number of regions and how concentrated is the variance on the first PCs. This is shown in Fig 3B. The amount of PCA distortion (Y-axis) is here quantified as one minus the correlation between the state time courses obtained on the original data vs. those obtained on the PCA projection. A measure of PCA concentration (X-axis) is given by the average cumulative explained variance across PCs; for example, if the areas were perfectly correlated then the first PC would explain all the variance (1.0) and the cumulative explained variance of all PC would be 1.0, in which case the average –i.e. the PCA concentration– would be exactly 1.0; in the opposite case, if the regions were orthogonal (uncorrelated), then the cumulative explained variance of the j-th PC becomes j/n, and the average becomes exactly 0.5. In each run, I sampled a number (between 5 and 100) of regions from the Schaefer parcellation and run the HMM on both the original and the PCA-projected data (with no loss of variance). Fig 3B shows that there is a logarithmic relation between PCA concentration and PCA distortion across HMM runs, suggesting that the more correlated the regions are, the stronger is the bias introduced by PCA.

It follows the question of whether HMM-PCA has this problem. It can be shown theoretically that HMM-PCA does not have this issue, because, as discussed, the problem has to do with having state-specific error variance parameters, while in our HMM-PCA formulation the error covariance matrix (and therefore error variance) is common to all states. As we have done empirically, we can manipulate the data by randomising their eigenvalues (i.e. by multiplying the ordinary PCA weights by a random number). Since this is changing the nature of the data, the HMM estimation, if correct, should also change to reflect the manipulation. If it does not change, that would signal a problem. Fig 3C shows the correlation between the HMM estimation (i.e. between the state time courses) before and after performing this manipulation in the data, for both HMM-Gaussian and HMM-PCA. As shown, HMM-PCA reflects that the data have changed by changing its estimates, while HMM-Gaussian is mostly unaffected. Although this can only be considered as indirect evidence, it adds further evidence to the issue discussed in this section.

In summary, although PCA is often an acceptable approximation in practice, it can also arbitrarily distort the time-varying FC estimates (with respect to a non-PCA estimation) towards the lower-order PCs. This effect will be more pronounced when the regions are more correlated –i.e. when the proportion of variance explained by the different PCA components is less equally balanced.

Empirical results

Next, I demonstrate the comparative performance of the models on both real data from the HCP and synthetic data constructed by using aspects of real HCP data –see [3] and references therein for details about the HCP data acquisition and preprocessing.

Simulated data experiments

In the previous section, I have used examples to illustrate the limitations of PCA when used in combination with the HMM or similar models. I now explore how the different models behave in higher-dimensional, more realistic data. Although these data are synthetic, I used aspects of real HCP data to perform the simulations.

Specifically, using the 100-regions Schaefer parcellation [30] and taking a standardised time series for each parcel, I computed the average FC across all subjects as a n × n global covariance matrix and eigendecomposed this FC matrix into its principal components = ′. Then, for each state, I randomly chose a set of principal components so that the total explained variance of this set does not go over a certain threshold ϵ. Next, I generated state-specific covariance matrices by permuting the values within each chosen eigenvector. I followed this procedure to generate six different states, each with a different covariance matrix . This equals 100 × 99/2 parameters per ground-truth state, for a total of 29700 parameters. This was done under two different conditions: for a smaller threshold ϵ = 0.1 and for a larger threshold ϵ = 0.2; i.e., in the ϵ = 0.2 case I permuted more eigenvectors than in the ϵ = 0.1 case, making the states more different to each other, and, therefore, making the subsequent HMM estimation easier to be performed accurately. Once I had the six ground-truth states, I generated ground-truth state time courses according to a transition probability matrix where the diagonal is 25 times higher than the off-diagonal, so that it is 25 times more probable to remain in the current state than to switch to a different one. Then, I sampled 100 subjects worth of data with 1000 time points each, where the sampling was done according to a Gaussian distribution with covariance given by the active state at each time point. Finally, I estimated HMM-Gauss, HMM-PCA and Mix-PCA models on this data set for a grid of number of states K = 4, 5, 6, 7, 8 and number of principal components p = 5, 10, 20, 30, 40. I repeated the entire process 10 times, calculating, for each model, data set and combination of parameters, (i) the cross-validated likelihood, and (ii) the accuracy with regard to the ground truth. The cross-validated likelihood, used as a quantitative way to assess the models and perform model selection, is sometimes preferred to other methods based on penalised likelihood when the model assumptions are too far from the true generating distribution of the data [31], as it is the case with brain data. Accuracy in this case was defined as how well the estimated state time courses could predict (in a least-squares sense) the ground-truth state time courses.

Fig 4 shows the results of the analysis. The top panels reflect the cross-validated likelihood. In all models, the cross-validated likelihood is able to estimate that at least 6 states are required, but it can hardly distinguish between 6 and higher values of K, suggesting that for 7 and 8 states the amount of model overfitting is negligible. For HMM-PCA and Mix-PCA, the cross-validated likelihood favours solutions with larger numbers of PCs. For HMM-Gaussian (which input data is ), we cannot straightforwardly compare different numbers of PCs, and neither can it be compared to HMM-PCA (which input data is ), since the likelihood is computed using different data for each choice of p (number of PCs). The middle panels show the error in predicting the ground state time courses, and the bottom panels show that same information as a function of the number of parameters. HMM-Gaussian and HMM-PCA have similar accuracies despite the fact that the ground-truth generating distribution is Gaussian. HMM-PCA outperformed Mix-PCA, specially in the harder ϵ = 0.1 case (right panels), highlighting the importance of the HMM temporal regularisation.

Fig 4

Empirical comparison of HMM-Gaussian, HMM-PCA and Mix-PCA on simulated data.

Two different conditions were tested: one with larger (A) and the other with smaller differences between states (B). The models were inferred for a range of states K and principal components p, and the entire process was repeated 10 times (averages are presented); see main text for details. The top panels show the cross-validated likelihood for each combination of parameters and model. The middle panels show the error of each solution (averaged across repetitions of the experiment) in predicting the ground-truth state time courses. The bottom panels show the errors as a function of model complexity for each model.

Overall, this section showed that, in simulated data mimicking some aspects of real data, HMM-PCA competes well with, or outperforms, alternative models. These simulations also show that the theoretical benefits of HMM-PCA presented above will apply to a lesser or greater degree depending on the characteristics of the data set under study and the models setting.

Real data experiments

Focusing on the HMM-based solutions, I next compared HMM-PCA and HMM-Gaussian on real resting-state fMRI data using 820 subjects from the HCP data set, where each subject underwent four 15-min sessions (TR = 750ms) in the scanner. I used a data-driven parcellation obtained through spatial independent component analysis (ICA) with 50 components; again, see [3] for details about preprocessing and the computation of the ICA time series. The time series of these ICA components were then standardised separately for each session, and then submitted to a stochastic variety of the HMM inference that is specially designed to deal with big volumes of data [16]. Since the estimation of the HMM parameters may return (for the same model and data) slightly different results for each run of the inference [32], I ran the inference five times per model, with K = 12 states and p = 24 principal components each. In brief, in what follows I show that the main FC patterns were not anatomically very different between the two approaches, but that these differences were behaviourally informative. This can be considered as indirect evidence of the superior sensitivity of HMM-PCA.

Asking whether the HMM-PCA states represent meaningful patterns of FC is not straightforward here because there is no ground-truth available. Since both HMM-Gaussian and HMM-PCA build on PCA, it is however expected that both approaches should be able to capture the main trends in the data to a relatively comparable extent. Given that HMM-Gaussian was shown to produce meaningful estimations in previous work [8, 21], proving that this is the case would situate HMM-PCA on first base. As an example, Fig 5A presents connectedness maps for two given states, where connectedness (or degree) is defined as how much each region correlates with the rest of the brain. The maps were centred across states, such that, if a region exhibits a positive value within a given state, then that region is more correlated to the rest brain’s voxels within this state than on average. One of the states is closely associated to the default mode network [2], and the other to the sensorimotor and visual systems. For these two states, both methods capture largely similar anatomical features. Based on the correlation between their FC patterns (specifically, by transforming the states’ covariance matrices into correlation matrices, taking the Fisher transformation, and then correlating the off-diagonal elements of these matrices between each pair of states), I then used the Hungarian algorithm [33] to match each HMM-PCA state to a HMM-Gaussian state. On the left, Fig 5B shows the resulting HMM-PCA vs HMM-Gaussian correlations for each pair of states, where the diagonal corresponds to states that were matched, and the off-diagonal to any other pair of states. On the right, Fig 5B shows the distribution of between-state correlations for states that were matched to each other (blue) and, to provide context, the correlations between states that were not matched (red). This indicates that the main time-varying FC patterns are relatively well preserved between the two methods.

Fig 5

Comparison of HMM-Gaussian and HMM-PCA on resting-state data from 820 Human Connectome Project subjects.

A. Two example states per model, default mode network and sensorimotor-visual, where the maps reflect the degree, i.e. the total amount of connectivity between each voxel and the rest of the brain. B. The states are relatively comparable across the two models. Left: Correlation matrix between the HMM-PCA and the HMM-Gaussian states (in terms of Pearson’s correlation between the off-diagonal elements of the states’ FC matrices) across 5 runs of the inference. Right: Distribution of the diagonal elements (red) and off-diagonal (blue) elements, where the red elements reflect the similarity between corresponding states (i.e. an HMM-PCA state and an HMM-Gaussian state that correspond to each other), and the blue elements reflect the similarity between non-corresponding states. C. HMM-PCA produces solutions that are more explanatory of behaviour. Using the fractional occupancy across runs (i.e. the amount of time spent on each state), the HMM-PCA models are better able to predict behavioural traits across a range of demographical, intelligence, personality and affective-related variables; each dot represent the explained variance (in terms of Pearson’s correlation) for the cross-validated prediction of one behavioural trait, and the bars represent the average across traits. D. HMM-PCA produces solutions that are more robust than HMM-Gaussian across repetitions of the inference; robustness is here measured as the capacity to predict (using cross-validation) the state time courses of one HMM estimate using the state time courses of another HMM estimate.

Then, I sought to investigate whether the theoretical benefits of HMM-PCA have indeed a practical impact on real data. First, for each run of the inference, I extracted the fractional occupancies, defined as the percentage of time spent on each state for each session. I then used these 12 × 5 fractional occupancy values as features to predict a collection of behavioural traits. In particular, I chose 63 traits across different domains including demographical, affective, personality- and intelligence-related [5], and performed cross-validated predictions respecting the family structure of the HCP data [34]. As shown in Fig 5C, the cross-validated predictions were found to be significantly more accurate for HMM-PCA than for HMM-Gaussian (p-value = 0.001, permutation testing). In accordance with previous work [35, 36] the accuracies in predicting HCP traits using FC features are relatively modest, with the average explained variance being lowered by traits that are particularly hard to predict (e.g. the personality traits). The average is, however, significantly higher than zero. Note the grey lines connecting each trait’s prediction between the two models reflecting that the two models tend to perform better on the same traits. Furthermore, as shown in Fig 5D, the HMM-PCA solutions were also more robust across runs of the HMM inference (p-value = 0.001, permutation testing); robustness in this case was measured as the capacity to predict, through cross-validation, the state time courses of a given HMM estimate using the state time courses of another HMM estimate (for a total of 20 pairs of estimates, having estimated the models five times).

These results were obtained running the HMM on a single hyperparameter configuration (K = 12 states and p = 24 principal components), so it is possible that the edge exhibited by HMM-PCA over HMM-Gaussian does not generalise to other configurations. To test this, I estimated the models on a grid of hyperparameters, K = 4, 5, 6, 7, 8 and p = 5, 10, 20, 30, 40, and reran the predictions on the behavioural traits for each pair of hyperparameters. Fig 6 shows the results of comparing the HMM-Gaussian vs the HMM-PCA predictions, where I performed (permutation-based) statistical testing in both directions: is HMM-PCA better than HMM-Gaussian at predicting behaviour? and, is HMM-Gaussian better than HMM-PCA at predicting behaviour? The left panels show a p-value for each combination of hyperparameters, and the right panels show histograms of p-values across all configurations. Finally, I used the nonparametric combination algorithm [32, 37], to combine the different tests across all combinations of parameters into a single, aggregated p-value. This procedure showed that HMM-PCA was significantly better at predicting behaviour than HMM-Gaussian for the considered grid of hyperparameters.

Fig 6

Behavioural predictions of HMM-Gaussian and HMM-PCA across a range of parameters.

Similarly to Fig 5, the fractional occupancy was used to predict a set of behavioural traits. Permutation testing was used to test whether, for each combination of parameters, the HMM-PCA solutions were better than the HMM-Gaussian solutions (top), and whether the HMM-Gaussian solutions were better than the HMM-PCA solutions (bottom). On the left side: the corresponding p-values, where colder colours mean statistically significant; at the center: histograms of p-values across all combinations of hyperparameters; on the right side: the result of combining all the tests into one single, aggregated test, using the non-parametric combination algorithm –as used in [32].

In summary, these results suggest, albeit in an indirect way, that HMM-PCA might be better able to describe time-varying FC in high-dimensional fMRI data, and that the theoretical limitations discussed above can in fact have an impact in practice.

Discussion

By fusing dimensionality reduction and time-varying FC in one single model, the presented HMM-PCA approach can bypass some important limitations of the common two-step procedure, where dimensionality reduction and time-varying FC estimation are performed in sequence. This was shown on simulations and on real fMRI data, where the HMM-PCA states were significantly better able to predict a number of behavioural traits. Importantly, the cause of the reported biases does not inherently lie on the high dimensionality of the data, but on the use of PCA as a way to reduce such dimensionality. Therefore, although for practical reasons the examples and simulations used throughout the paper are not high-dimensional by normal standards [38], the issues addressed by these examples apply generally.

An important question about the application of these models to real fMRI data is how to decide which model is objectively better. In Fig 3, for example, I showed that the use of PCA changes the estimation from what it would be obtained without using PCA. This is a mathematical fact, but does it mean that the PCA-related prediction is worse? This question depends on what is meant by the goodness of the model. In many practical applications, goodness would be given by some combination of accuracy (here, data likelihood) and model complexity [19]. Here, I used the cross-validated likelihood to assess the models, which can be more appropriate than methods based on the penalised likelihood (like the free energy) when the assumptions of the compared models do not meet the realities of the true generating distribution [31]. Most importantly, although the HMM has modelling assumptions, these do not imply a statement about the underlying biology. For example, using eight states does not necessitate or imply the assumption that there are eight biological states in the human brain (or eight attractors in the system). For this reason, we cannot straightforwardly say that a model is more biologically plausible because it has a higher cross-validated likelihood. Still, a valid theoretical argument can be made: in the asymptotic limit of having infinite data and assuming an ideal parameter inference (i.e. such that overfitting is not a factor), any transformation (like PCA) that changes the estimation from the one obtained from the original (non-PCA) data is undesirable. From Fig 3, it follows that a PCA transformation can exert such detrimental bias, because PCA alters the estimation independently of the amount of data. The extent to which this actually happens in real fMRI data will depend on the data and the preprocessing. For instance, the bias highlighted in Fig 3 would be less severe in data where global signal regression [39] has been performed; otherwise, a single component capturing the global signal could explain a large amount of variance, therefore increasing the PCA concentration measure and therefore the distortion. A more complete description of which data sets would be more or less affected by these issues would possibly necessitate the use of simulations that are more biophysically realistic [40].

Another important question is what these models are actually expressing. Because PCA is computed on the full data covariance, it is theoretically possible for the states transitions to be driven only by changes in the pattern of relative variance across channels with little or no contribution of the between-regions covariance. Here, however, it could be observed that the covariance between regions (above and beyond the variance) is clearly different between states. Otherwise we would not see any meaningful patterns after normalising by the variance (i.e. after transforming the covariance matrices into correlation matrices). The connectivity maps shown in Fig 5A, as well as the between-state FC correlations shown in Fig 5B, support that there are within-session FC modulations driving the estimation. This is in agreement with our previous work on the HCP data, where we showed (i) that HMM states have unique information in the off-diagonal elements of the covariance matrix, and (ii) that HMM models based solely on the variance (i.e. with no covariance) did not predict behaviour as accurately as those that model the full covariance [5, 8].

A limitation of these models is that they might not be suitable to be run in very high-dimensional data, such as the original surface-space of the HCP data (around 90k vertices/voxels). On top an excessive computational cost, without having appropriate additional mechanisms to control the high-dimensionality of the parameter space, overfitting might occur due to the large number of parameters in each state. Indeed, if the data dimensionality n is very high, HMM-PCA (with n × p parameters per state) might potentially overfit more than HMM-Gaussian (with p × (p − 1)/2 parameters per state), overshadowing the advantages discussed in this paper. Altogether, the presented model and the considered alternatives are, at present, better suited to be run in intermediate spaces, such as those produced by ICA or an anatomical parcellation, than on raw whole-brain space. Of course, each of these choices entail its own biases. An efficient application of prior distributions on the state-specific PCA weights , for example promoting sparsity [19, 41], is a promising avenue to help in this direction. Other limitations, such as the incapacity to model long-term temporal dependencies or the fact that the model does not allow for overlapping states, are however not specific of HMM-PCA, but more generally of the HMM framework.

Conclusion

In this paper, I have addressed the question of how to estimate patterns of time-varying FC in fMRI data. I have shown that the standard approach of sequentially applying PCA and then feeding the resulting PCs to an HMM or mixture model, although useful in practice, may suffer from biases and loss of sensitivity. On these grounds, I have introduced a new variety of the HMM, namely HMM-PCA, where each state is a probabilistic PCA decomposition. Critically, the HMM states not only express a (linearly optimal) dimensionality-reduction of the data, but also encode a correlation pattern between regions.

Graphical representations of the four models.

(TIFF)

Click here for additional data file.

6 Jan 2021

Dear Dr Vidaurre,

Thank you very much for submitting your manuscript "The statistical challenge of finding spontaneous changes in functional connectivity in high-dimensional fMRI data" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers.

Your paper was overall very well received. Some aspects need to be better clarified and inserted in the state of the art. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Daniele Marinazzo

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Summary:

In this paper, the author proposes an HMM-PCA model to explain time-varying FC in resting state data. The idea is very nice and potentially very useful, but the methodology and results do not adequately support the claims. The paper needs to be improved.

General feedback:

To motivate this approach, I think the author needs to clearly show the generative model (in form of a graphical model) and also report the model evidence, when performing model inversion. As for the inference technique, if one commits to maximum likelihood, then a grid search over the number of HMM hidden states and the number of PCA latent states is essential. Alternatively, if one turns to Bayesian PCA (reference 19), then Bishop's hyperparameter approach will determine the effective dimensionality of PCA, but the optimal cardinality of HMM’s latent space would still rely on grid search. Either way, the model evidence (accuracy minus complexity) will determine whether this new model is superior to HMM+PCA (in 2 steps) or Gaussian mixture model, in the present application. Ideally, this superiority needs to be demonstrated on high-dimensional (simulated) data first, because the claims pertain to the curse of dimensionality. In short, face validation of a new model is conventionally established based on model comparison (in terms of model evidence) and accuracy of the posteriors (with respect to the ground truth). If there is any such thing as bias/distortion or loss of sensitivity (as the author suggests) due to the two-step procedure, these will also show up in the results, alongside the more important model evidence.

Detailed comments:

- Speaking of the issue of high dimensional data, please note that 50-100 regions/ICs is not really high dimensional when the time series are concatenated across subjects. Moreover, your simulations are low-dimensional (with 10 regions/ICs). So you are not actually tackling the high-dimensional scenario (of many voxels and few time points) neither in the simulations nor in empirical data, to show the reader the gravity of this problem and how you have solved it.

- page 2, 2nd paragraph: this criticism to the sliding window analysis (that mere fluctuations may be mistaken for non-stationarity in the connectivity dynamics) applies to the hidden Markov model as well, unless model evidence rules out the existence of only one underlying connectivity state.

- page 3, 1st paragraph: Please rephrase. Using ambiguous statements such as "works reasonably well in practice for many (but not all) fMRI data" without bringing in quantitative support or a credible reference is not appropriate. Or, for instance, saying that a method "is not designed for the final goal of the analysis" without specifying that final goal, is not precise enough for a scientific text.

- page 3, 2nd paragraph: I can't see why the FC matrices are being interpreted as instantaneous, when in practice they are the covariance matrices of the Gaussian distributed observations (of each state).

- page 3, 4th paragraph: principal components are by design orthogonal. I don't see your point about transient departure from orthogonality. Please elaborate if this is actually what you meant.

- page 3, section 2.2: when a new probabilistic model is proposed, it is common to show the (generative) model structure as a graphical model, clearly showing the conditional (in)dependencies, and the model parameters and hyperparameters.

- page 4, 1st paragraph: It seems that you are using probabilistic PCA, not Bayesian PCA (with prior over W), as per your reference 19 (Bishop, Bayesian PCA).

- page 4, below Eq.1: I think discounting the mean needs further motivation. Does it have computational or practical advantages? More importantly, do you get higher model evidence by fixing all the means at a prior value of zero?

- page 4, middle of page: "Note that the estimation of u_tk only depends on the data at time point t and on pi_k, ignoring the data temporal structure". It is true that temporal order is not accounted for in mixture models; however the responsibilities depend on all the data points through the model parameters. So each u_tk does not depend *only on time point t and on pi_k*. Please rephrase to avoid misunderstanding.

- page 4, middle of page: "In the HMM-PCA case, we instead have the state latent variables xk modelled as an order-1 Markovian process...". Please draw the generative model to show how the HMM states map to the PPCA states, which are then combined with W to generate the data dt. Also please make sure the reader realizes the difference between these two types of latent states; i.e. the (categorical) hidden states of HMM, as opposed to the (Gaussian) latent states of PPCA. It is useful to stress their dimensionalities as well (K versus p0).

- page 4, Eq.3: Complete data log-likelihood (LHS) is also a function of the model parameters.

- page 4, Eq.4: this is not the familiar form of Bishop's PPCA. If it comes from another source, please provide the reference. Otherwise please put your derivation in the appendix.

- page 6: This whole bias/distortion story of PCA (assuming it is serious) has not been demonstrated well. Note that in Fig 2B most of the computed distortion is below 2%. Plus, how often does the first PC become so dominant in a FC pattern to take up 86% of the variance? I think to motivate your approach, you can compare model evidence of HMM+PCA (as a two-step procedure) to your proposed HMM-PCA. Even Fig 2A is not convincing, because there is no ground truth here, so we can not say PCA has caused any distortion. HMM+PCA is simply suggesting a different number (and sequence) of hidden states for the HMM, which may actually be more plausible than the results on the original (unrotated) data, in terms of model evidence. Note that you have not shown model comparison results for the number of hidden states, neither on the original data, nor on the PCA-rotated data. So I think this exemplar plot is not sufficiently supporting the distortion/bias claim.

- page 6, section 2.3.2: the same problem (as above) exists when discussing loss of sensitivity. Although this claim comes with simulated data, still the dimensionality is low (10 regions) and there is no model comparison result for different number of HMM hidden states and different number of PCA components. Accuracy of the state sequence through time is not adequate support for superiority of one method over another, when the correct number of latent states have been provided and no model evidence is reported. And even these accuracy values have not been statistically compared. Please reconsider these simulations.

- page 8, section 3: Same problem. Maximum fractional occupancy and behavioral relevance can not replace model evidence. And please note that the behavioral relevance is barely different from zero. Perhaps if you compute the effect size alongside the p-value it becomes clearer that this may not be a very strong argument in favor of your model.

This model has great potential and it deserves to be presented and motivated accordingly. Best of luck.

Reviewer #2: This article concerns the statistical difficulties in estimating time-varying functional connectivity in the brain from fMRI images where spatial dimensionality is high but the temporal windows for the time-resolved networks are small. The authors show that conventional time-averaged PCA methods can exhibit bias towards the least relevant components because all are given equal weight. They propose a hybrid method which uses (i) PCA to reduce the data at each time sample, and (ii) HMM to identify changes in the state of the FC network. The proposed method is tested on simulated and real fMRI data.

The paper is very well written overall. I especially appreciated the Introduction as I am not familiar with these statistical methods. Consequently my review is rather general in nature. I garnered that the proposed method improved the statistical result by using a single-step procedure to simultaneously reduce the data and estimate the functional connectivity at each state. However it was not clear to me how the algorithm worked. Specifically, I did not understand where the K different state-specific covariance matrices originated. Do they come out of the analysis?  Is the number K something that must be chosen beforehand? Perhaps the author can clarify these basic issues for readers who are not familiar with Mix-PCA.

The layout of the manuscript could be improved. It is not clear where the Methods end and the Results begin. Methods sections 2.3 onwards would be better placed under Results.

The results themselves are well-presented and I have no issues there.

Section 2.3.1 provides a good explanation of why the lower-order components in time-averaged PCA introduce biases. The third paragraph "When we apply PCA..." summarizes it particularly well.

Section 2.3.2 provides nice results for the better sensitivity of the HMM-PCA method.

Section 3 compares the methods on real fMRI data and finds better results for the proposed method. The author admits (page 9 top) that is difficult to interpret the results in the absence of a ground truth in real fMRI data. Might it be possible to use a generative model of fMRI to test the algorithms in a scenario where the ground truth is known? Perhaps the author can comment on this in the Discussion.

Upon finishing the paper, I felt that it could have done more to make the proposed method available to the reader as a software package. At very least I encourage the author to make his existing code available as a demonstration. Otherwise the reader is obliged to transform the mathematical equations into their own software, which can be risky and time-consuming.

Minor points

----------

Section 2.1, second paragraph. "... as zero mean Gaussian distributions."

Gaussian distributions of what? I didn't understand this statement.

Section 2.1, third paragraph. "... transient departures from orthogonality."

I was confused what you meant by this. Do you mean that the principle components of one state differ from those of the next state? Can you please clarify?

Section 2.2, first paragraph.

(1) Typo: "Various [of the] elements [of] this model are ...."

(2) How is K determined? Please elaborate.

Section 2.2, second paragraph. "the Mix-PCA model uses K different PCA projections"

Again, I wasn't clear to me where these K different projections came from.

Section 2.2, second paragraph. "... using the EM algorithm"

What does "EM" mean? Perhaps I missed it but I don't think it was ever explained.

Section 2.2, third paragraph.  "...using the forward-backward equations ...."

What are the 'forward-backward' equations? Are these a standard part of the Mix-PCA method?

Section 3, Title: 'Results on real fMRI data'

This title is a little misleading as it suggest that the previous results (section2) are not based on real fMRI data. Whereas Fig 2 is based on real fMRI data, albeit manipulated data. Likewise for Fig 3, Scenario 2.

Reviewer #3: In this paper the author has suggested a new approach in the HMM family that tries to combine PCA with HMM in a more unified way. The aim of this approach is to remedy the issue of high dimensionality of fMRI data. I find the method quite interesting and can think of many instances where such an approach would be quite useful. My major issue with this paper is that it seems a little rushed (for the lack of better term) in both introduction and discussion. I think if the authors add some more information to both introduction and discussion, the paper would be improved greatly.

Specific comments:

1. First and foremost, I do not find the title of the paper representing the core of this paper. This paper does not discuss the challenge of finding spontaneous changes in functional connectivity” a lot. Instead, a new approach is proposed that might help solve some of those “challenges”. When I read the title, I thought this paper is a review/commentary paper.

2. I also find the introduction overly short. The authors have only given sliding window Pearson correlation (SWPC) as a different category of methods (the rest are all HMM). First although SWPC is one the most know method in this field, there are some more recent connectivity estimation that have been proposed recently. For example, shared trajectory by (Faghiri et al. 2020) or instantaneous phase synchrony (Pedersen et al. 2018) and even Multiplication of Temporal Derivatives (Shine et al. 2015). As these are more instantaneous estimators (compared to SWPC where the window size makes sure the estimand are not instantaneous) I find them more relevant to the current work and therefore should probably be cited too. Second, I don’t think these estimators should be compared to HMM at all. As HMM tries to directly estimates the FC states whereas these estimators aim to estimate connectivity and an additional step is usually required to estimate FC states from these connectivity time series. For example, one can use kmeans to estimate the FC states from connectivity time series (Allen et al. 2014). Right now if someone reads the introduction it seems you are comparing SWPC with HMM family which is not quite accurate.

3. In the “Loss of sensitivity” section the author states “for example, this would correspond to a data set with sampling rate or TR≈1.1s, and 15min worth of data per session.” I do not understand how the authors has calculated the TR. Based on what I understand no sampling/resampling has been done to get a sampling rate value. This simulation can be from any other TR value too. And why is this information is given at all? The authors do not use it any other place to derive any kind of result.

4. Is it possible that the stability of the other HMM variants (gaussian HMM) is lower compared to the HMM PCA? i.e. can gaussian HMM predict better if it is run more than 5 times?

5. I do not see any mention of pre-processing of the fMRI data. Even if the authors have not done this at least a short summary would be needed (or a citation of where to find that information).

6. I find the discussion of this paper quite underwhelming. I think the author need to expand this part too. for example, a little more on the limitation of the proposed method? I see that the authors mention that their method is sensitive to the temporal ordering of the time series (which is true) but it has a limit. If I am understanding it correct it only care about temporal ordering for one time point. This can be considered a limitation as I believe this method would be unable to see slower changes. This point to another limitation of this study. The proposed method only find hard changes in the state and not the more fuzzy changes (at least the simulation is designed like that).

References:

• Faghiri, A., Iraji, A., Damaraju, E., Belger, A., Ford, J., Mathalon, D., ... & Turner, J. (2020). Weighted average of shared trajectory: A new estimator for dynamic functional connectivity efficiently estimates both rapid and slow changes over time. Journal of neuroscience methods, 334, 108600.

• Pedersen, M., Omidvarnia, A., Zalesky, A., & Jackson, G. D. (2018). On the relationship between instantaneous phase synchrony and correlation-based sliding windows for time-resolved fMRI connectivity analysis. Neuroimage, 181, 85-94.

• Shine, J. M., Koyejo, O., Bell, P. T., Gorgolewski, K. J., Gilat, M., & Poldrack, R. A. (2015). Estimation of dynamic functional connectivity using Multiplication of Temporal Derivatives. NeuroImage, 122, 399-407.

• Allen, E. A., Damaraju, E., Plis, S. M., Erhardt, E. B., Eichele, T., & Calhoun, V. D. (2014). Tracking whole-brain connectivity dynamics in the resting state. Cerebral cortex, 24(3), 663-676.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: Yes: Ashkan Faghiri

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions, please see

11 Feb 2021

Submitted filename: ResponseLetter.docx

Click here for additional data file.

28 Feb 2021

Dear Dr Vidaurre,

Thank you very much for submitting your manuscript "Fusing dimensionality-reduction and time-varying functional connectivity estimation in a single model" for consideration at PLOS Computational Biology.

While we appreciate the way you addressed some issues, we feel that some important issues, as identified by reviewer 1, are still present, and to some extent more evident after this round of revision.

Some of these concerns are fundamental, and as such they could be unsurmontable. Still I think that it would be good to go through them.

This field is rapidly advancing, and we appreciate that you are providing a robust framework with modern tools. OIn the other hand it should be clear which aspects are novel with respect to previuous work (even beyond neuroimaging data analysis), and how the results can be interpreted.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Daniele Marinazzo

Deputy Editor

PLOS Computational Biology

Daniele Marinazzo

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: I appreciate that the manuscript has been improved with this revision. Performing more simulations, showing graphical models, computing model evidence, plus a number of other clarifications have all been valuable additions. I agree that HMM-PPCA is a potentially useful addition to the field of dynamic FC analysis.

However, there are still fundamental issues that we disagree about. In the following, I will elaborate on these points.

• Introduction/motivation

To start with, please remember to cite the HMM-mixture-of-Bayesian-PCA model of Mauricio Alvarez and Ricardo Henao (2007) as prior work. Actually your HMM-PPCA model seems like a special case of their formulation.

From an organization point of view, the motivation section (called limitations of two-step approach) should go under Introduction, and not Results. But more importantly, I don’t think the motivation content is appropriate.

Specifically, I still can not sympathize with the distortion story of PCA, used to motivate your model. You are saying that two-step PCA+HMM produces results different from HMM on raw data, when many PCs are included. Lets forget for now that you showed us this distortion is about 1% for 84% explained variance. I repeat what I mentioned in the previous review: I don’t understand why anyone would expect PCA+HMM results to be identical to HMM. PCA is a data transformation/reduction technique. Of course it will discard certain details to facilitate discovery of the underlying structure in the data. This is not distortion, it’s the natural compromise of data reduction.

If you still see this as an important problem, have you at least shown that HMM-PPCA solves it? Have you shown that HMM-PPCA results are very similar to HMM on raw data? I didn’t find any such results. When you stress an issue as the motivation of your work, the reader is waiting to see how you have resolved this issue.

I understand you need to motivate the model. Why not use all the motivations in the mixPPCA/mixFA literature? For example: the benefits of local linear sub-models, having fewer parameters to estimate than a full-covariance Gaussian mixture model, not having to impose axis aligned covariance matrices to reduce the number of parameters, all the advantages of probabilistic PCA over classic PCA, the temporal structure modeled with HMM, etc. I am sure you know the following references, so you can easily get inspirations from their introductions: Bishop and Tipping (1999), Bishop and Winn (2000), Ghahramani and Beal (1999).

If you have been thinking that global PCA prior to HMM is not as efficient as local PCA (using mix-PPCA or your extension as HMM-PPCA) in capturing and preserving local structure, I totally agree. But this is not what you have explained for the reader.

o Simulation Design

Simulations should be designed to showcase the important aspects of a proposed model. In the new high-dimensional simulations of the manuscript, the (effective) dimensionality of the data, per component/state, is not specified. This is important because you need to show that HMM-PPCA can recover both the number of HMM latent states and the intrinsic dimensionality of each state-specific PPCA. Actually your lower-dimensional simulations had this aspect encoded in the number of retained eigenvectors of the covariance, but the new simulations miss this point. When we don’t know the dimensionality of the principal subspace in the simulations, we can not judge the model inversion and comparison results.

o Coupled update equations

(Expected) complete data log-likelihood (Eq. 4) should be written for the whole HMM-PPCA model (not per state), in which case the EM update equations of PPCA and HMM become coupled. Intuitively, this is because the likelihood term (used to update the HMM state) now depends on PPCA posteriors as well. This is an important point to show, as you are offering a unified model. Check out the following references for similar derivations of coupled update equations:

- Bishop and Tipping's (1999) mix-PPCA (Appendix C, Eq. 73)

- Bishop and Winn’s (2000) variational treatment of mix-BPCA (Section 3.1)

- Alvarez and Henao’s (2007) derivations for the HMM-mix-BPCA model (Appendix B)

o Latent space identification

My concerns here can be best summarized as a question: How does the author intend to identify the optimal number of latent states for HMM and PPCA on empirical data? There is no answer in the manuscript.

There are standard approaches in the literature for identifying the number of latent states in such models: grid search, cross-validation, and Bayesian treatment (with hyper-priors over the loading matrix). Check out for example: (Bishop and Winn, 2000; Bishop and Tipping, 1999; Alvarez and Henao, 2007; Ghahramani and Beal, 1999) to see how the latent space dimensionality is specified and then recovered successfully. While the fully Bayesian treatment has undeniable advantages, all three solutions are valid. My point is that this problem has been solved in the literature, and your application is no exception.

Also note that the dimensionality of the principal subspace can differ across states/components in the HMM-PPCA model. In your application this means that different functional states can express their intrinsic complexities through the dimensionality of their respective principal subspaces.

So, on the one hand, optimizing the latent space is essential for model fitting on a given dataset. On the other hand, this step is important for between-group analysis, e.g. between a normal and patient group. Once the model has been fitted to each dataset separately, one may find out that the functional components from one group fit into lower dimensional subspaces, or may need fewer HMM states to be explained, than the other group, which can be nicely interpreted. If latent variable dimensionalities have not been optimized on each dataset, such between-group analysis would not be feasible.

I can see that you tried the grid search approach, and computed free energy as a proxy for model evidence. This is a valuable step in the right direction, because free energy computation can guide the choice of model parameters and also facilitate between-model comparison (e.g. HMM-PPCA vs. mix-PPCA). Apparently you think free energy has failed you in model comparison. I discuss this next.

• Model comparison with free energy

Please note that what I call free energy (F) is the negative of free energy in physics, or ELBO in machine learning. There could be a number of reasons why free energy is not guiding your model comparison to optimize the latent space dimensionality and to compare competing models of different structure:

1) I suspect that you have not computed free energy for the whole hierarchical model. I had a glance at your code, and F seems to pertain to the HMM part of the model only. I can’t see the complexity terms neither for the PPCA nor for the transition matrix. If that is indeed the case, it explains why F is not penalizing higher complexity and keeps encouraging more and more dimensions/states (similar to over-fitting of maximum likelihood schemes). You may want to have a look at the Appendix of Friston et al. (2015) to see how the free energy of a hierarchy is written in terms of the expected free energy of the lower level minus the complexity of the higher level.

2) A second problem might be the way you are parametrizing the likelihood. I remember I mentioned in the previous review that taking out the ‘mean’ is not justified. That is because when you fix all the means to zero, it’s like trying to disentangle concentric ellipsoids. You may want to put the mean back and see whether the state cardinality shows itself.

3) A third solution is adopting a fully Bayesian treatment (with prior over PPCA loadings). At least for PPCA dimensionality discovery this approach has proven very successful in the literature.

• Replacements for model evidence

The author seems to propose fractional occupancy of states and behavioral relevance as model evaluation criteria. As I explained before, these can not replace model evidence. (As a side note: one can use other approximations to model evidence, other than F, or even use sampling to compute model evidence).

Speaking of fractional occupancy, I have seen in the literature that a balanced fractional occupancy of states is more probable in healthy brains, than pathologic ones. However, fractional occupancy is not the cost function for optimizing the parameters of your model. This sort of domain-specific information can be encoded as ‘prior’ belief, e.g. on the form of the transition matrix. As for behavioral relevance, only once the model optimization is over can one inspect relation to behavior. Comparing behavioral relevance on (structurally different) models with arbitrarily chosen parameters is not meaningful, to claim superiority of one model over another.

I think taking care of these points would go beyond a standard revision. I would encourage the author to go through the reviews, repeat some of the analyses, and resubmit a fresh version of this work. I still think the idea is great and should definitely be added to the literature soon. Best wishes.

• References:

Alvarez, M. and Henao, R., 2007. Hidden Markov Bayesian Principal Component Análisis. Lecture Notes in Computer Science, Neural Information Processing/ICONIP.

Bishop, C.M. and Winn, J.M., 2000, June. Non-linear Bayesian image modelling. In European Conference on Computer Vision (pp. 3-17). Springer, Berlin, Heidelberg.

Tipping, M.E. and Bishop, C.M., 1999. Mixtures of probabilistic principal component analyzers. Neural computation, 11(2), pp.443-482.

Ghahramani, Z. and Beal, M.J., 1999, December. Variational Inference for Bayesian Mixtures of Factor Analysers. In NIPS (Vol. 12, pp. 449-455).

Friston, K., Zeidman, P. and Litvak, V., 2015. Empirical Bayes for DCM: a group inversion scheme. Frontiers in systems neuroscience, 9, p.164.

Reviewer #2: The authors have addressed all of my concerns.

Reviewer #3: I think the author have answered all my concerns. the paper is in a good state to be published.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: Yes: Ashkan Faghiri

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions, please see

21 Mar 2021

Submitted filename: Review2.docx

Click here for additional data file.

31 Mar 2021

Dear Dr Vidaurre,

We are pleased to inform you that your manuscript 'Dimensionality reduction and time-varying functional connectivity estimation in one single model' has been provisionally accepted for publication in PLOS Computational Biology.

When you submit the final version, please fix the typo pointed out by the reviewer.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Daniele Marinazzo

Deputy Editor

PLOS Computational Biology

Daniele Marinazzo

Deputy Editor

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: I would like to thank the author for his scientific attitude in taking the reviews. I think switching to cross-validation was a wise decision for a maximum-likelihood scheme. The other revisions have also improved the paper and made it clearer.

Just a typo: Eq (4) last line: is right, (x_0, not x_t)

I recommend this version for publication. Best wishes to the author and congratulations.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

9 Apr 2021

PCOMPBIOL-D-20-02157R2

Dimensionality reduction and time-varying functional connectivity estimation in one single model

Dear Dr Vidaurre,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Katalin Szabo

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol