Sevoflurane Alters Spatiotemporal Functional Connectivity Motifs That Link Resting-State Networks during Wakefulness

Background: The spatiotemporal patterns of correlated neural activity during the transition from wakefulness to general anesthesia have not been fully characterized. Correlation analysis of blood-oxygen-level dependent (BOLD) functional magnetic resonance imaging (fMRI) allows segmentation of the brain into resting-state networks (RSNs), with functional connectivity referring to the covarying activity that suggests shared functional specialization. We quantified the persistence of these correlations following the induction of general anesthesia in healthy volunteers and assessed for a dynamic nature over time. Methods: We analyzed human fMRI data acquired at 0 and 1.2% vol sevoflurane. The covariance in the correlated activity among different brain regions was calculated over time using bounded Kalman filtering. These time series were then clustered into eight orthogonal motifs using a K-means algorithm, where the structure of correlated activity throughout the brain at any time is the weighted sum of all motifs. Results: Across time scales and under anesthesia, the reorganization of interactions between RSNs is related to the strength of dynamic connections between member pairs. The covariance of correlated activity between RSNs persists compared to that linking individual member pairs of different RSNs. Conclusions: Accounting for the spatiotemporal structure of correlated BOLD signals, anesthetic-induced loss of consciousness is mainly associated with the disruption of motifs with intermediate strength within and between members of different RSNs. In contrast, motifs with higher strength of connections, predominantly with regions-pairs from within-RSN interactions, are conserved among states of wakefulness and sevoflurane general anesthesia.


DYNAMIC CORRELATION MODEL
Let z k , k = 1, 2, ... be a vector time series of correlation values between pairs of brain regions. We assume that z k arises from an underlying state x k ∈ R N via where the function F −1 is vector-valued inverse Fisher transform which can be written as Such a function bounds the state within [−1, 1] and is the standard transformation of the cross-correlation into a normally distributed variable. We assume that x k obeys a standard, linear state space model of the form x k = A k x k−1 +w k−1 , Here, the state noise process w k ∈ R N is an N -dimensional random vector with multivariate Gaussian distribution having zero mean and covariance matrix Q w ∈ R N ×N , the noise vector v k ∈ R N is a zeromean multivariate Gaussian random vector with covariance matrix R v ∈ R N ×N . Bounded nature of correlation coefficients result in a non-Gaussian noise for the observation vector y k . The covariance Q w determines the extent to which correlation can change in successive measurements, while R v simply characterizes measurement noise. In this sense, the vector y k consists of noisy observations of pairwise correlation values. We assume that state noise and measurement noise are uncorrelated.
We seek the optimal filter for obtaining the state estimateẑ k in the sense of minimum mean-squared error Central to this problem is the calculation of the probability density function (p.d.f.) of the state vector at any given time k, conditioned on Y k = {y 0 , y 1 , . . . , y k }, the set of all the past observations. The presence of the nonlinearity (2) complicates this calculation, but only slightly since F is smooth and invertible. Thus, it follows immediately that we can obtain a surrogate measurement such that d k is linear in the state x k . Note that (1) directly yields the p.d.f. of z k from that of x k . Thus, the problem (4) reduces to the classical Kalman filter to obtain the (Gaussian) p.d.f. of x k given Y k (Kalman, 1960b,a), based on the measurements d k .
It follows directly from (3)-(5) that the p.d.f. p(y k | x k ) can be written as where

a zero-mean multivariate Gaussian with covariance matrix
R v . A standard Bayesian approach thus yields the posterior density p (x k | Y k ), which is Gaussian with covariance matrix (Σ) and mean (µ) as The Kalman update equations are thuŝ In the subsequent results, we assume an initial multivariate Gaussian prior x 0|−1 with meanx 0|−1 and covariance of P 0|−1 . We note that, for the problem (4), (8) returns the optimal estimate. After obtaining state estimatex k|k , the correlation coefficientsẑ k are approximated using fourth order Taylor expansion of (1) aroundx k|k aŝ where all the multiplication and powers of vectors are element-wise operations. Also, diag() operator applies to a square matrix and returns a vector containing the diagonal elements of the matrix.
A final obstacle in using the filter equations (8) is knowing the noise parameters that specify the statespace model. Since such information is difficult to define a priori, a principled way to update them alongside the state estimates is required. Here, we face the challenge of estimating noise parameters alongside our online filtering. To address this issue we use an EM algorithm (Ghahramani and Hinton, 1996;Shumway and Stoffer, 2010) to jointly estimate the state and model parameters from observed correlation data.
Assuming diagonal covariance matrices for model parameters, we can decouple the update procedures for each edge, resulting in a reduction of the computational burden typically associated with EM methods.
Briefly speaking, after initializing the model parameter (process and measurement noise), we perform the E-step, filtering and smoothing, and then updating the parameter in the M-step. The E and M steps alternate to convergence. In order to obviate computational issues associated with this algorithm, a maximum iteration number is set alongside the convergence criterion.

DETERMINING THE NUMBER OF MOTIFS USING AIC
Selecting the number of clusters, K, in unsupervised learning is one of the key problem. The correct choice of K is often ambiguous. There exists several heuristic methods to facilitate the selection of k such as "The Elbow Method" by (Thorndike, 1953) and a technique based on the use of hierarchical clustering by (Salimi-Khorshidi et al., 2011). Another set of approaches for determining k are information criteria methods. In this paper, we adopt the Akaike information criterion (AIC) (Akaike, 1974) approach.
AIC rewards goodness of fit, assessed by the likelihood function, but it also includes a penalty that is an increasing function of K. Figure 1 shows AIC versus K for both 0% and 1.2% vol sevoflurane. In this plot, AIC values are normalized to have maximum value of unity for each condition. It can be seen from this figure that K = 8 is the optimal choice of K for both condition.

SPATIOTEMPORAL ANALYSES OF SYNTHESIZED DATA
Here, we demonstrate the efficacy and utility of the approach in a simple example of a data set with 9 locations and 160 time windows. We assume there exists 5 spatial motifs shown in Figure 2A and each motif has the same temporal correlation trajectory with the same mean over time (Green lines shown in Figure 2B).
Before decomposing the correlation trajectories, we added noise drawn from uniform distribution inside interval (−0.1, 0.1) to all region pair correlations over time. Figure 2 shows spatiotemporal decomposition of correlation time series. Spatial and temporal motifs are shown in Figures 2A and 2B, respectively. It can be conclude from these plots that temporal information can play an important role for extracting meaningful motifs even though if all the region pairs have the same static correlation value.

ROBUSTNESS AS A FUNCTION OF INTERSUBJECT VARIABILITY
To see how consistent the spatiotemporal motifs are for different conditions (0% and 1.2%), we look at how spatial motifs at 0% and 1.2% obtained from each individual, without concatenating data from all subjects for each condition, are similar. Figure 3A demonstrates the average cardinality of each motifs over different individual. Error bars stand for SEM of the average cardinality. Figure 3B represents average similarity between motifs at 0% and 1.2% obtained from each individual for 11 seconds window size. Standard error of the mean of the average similarity is shown in Figure 3C. We can see from this figures that the similarity between first, second and eighth motifs are robust over different individual.