==== Front ArXiv ArXiv arxiv ArXiv 2331-8422 Cornell University arXiv:2305.19573v1 2305.19573 1 preprint Article Reliability of energy landscape analysis of resting-state functional MRI data Khanra Pitambar 1 Nakuci Johan 2 Muldoon Sarah 14 Watanabe Takamitsu 3 Masuda Naoki 14* 1 Department of Mathematics, University at Buffalo, State University of New York, Buffalo, USA. 2 School of Psychology, Georgia Institute of Technology, Atlanta, USA 3 International Research Centre for Neurointelligence, The University of Tokyo, Japan 4 Computational and Data-Enabled Science and Engineering Program, University at Buffalo, State University of New York, Buffalo, USA * corresponding author (naokimas@buffalo.edu) 31 5 2023 arXiv:2305.19573v1https://creativecommons.org/licenses/by/4.0/ This work is licensed under a Creative Commons Attribution 4.0 International License, which allows reusers to distribute, remix, adapt, and build upon the material in any medium or format, so long as attribution is given to the creator. The license allows for commercial use. nihpp-2305.19573v1.pdf Energy landscape analysis is a data-driven method to analyze multidimensional time series, including functional magnetic resonance imaging (fMRI) data. It has been shown to be a useful characterization of fMRI data in health and disease. It fits an Ising model to the data and captures the dynamics of the data as movement of a noisy ball constrained on the energy landscape derived from the estimated Ising model. In the present study, we examine test-retest reliability of the energy landscape analysis. To this end, we construct a permutation test that assesses whether or not indices characterizing the energy landscape are more consistent across different sets of scanning sessions from the same participant (i.e., within-participant reliability) than across different sets of sessions from different participants (i.e., between-participant reliability). We show that the energy landscape analysis has significantly higher within-participant than between-participant test-retest reliability with respect to four commonly used indices. We also show that a variational Bayesian method, which enables us to estimate energy landscapes tailored to each participant, displays comparable test-retest reliability to that using the conventional likelihood maximization method. The proposed methodology paves the way to perform individual-level energy landscape analysis for given data sets with a statistically controlled reliability. Maximum entropy model Ising model functional magnetic resonance imaging Bayesian approximation Permutation test Fingerprinting ==== Body pmc1. Introduction Brain activity is dynamic and nonlinear in nature. Such nonlinear brain dynamics are considered to underly many functions of the brain such as cognition, action, and learning [1–3], and mathematical modeling is widely accepted as a useful tool for simulating such brain dynamics on different scales [1–7]. There are also many methods for analyzing empirical data of neural dynamics, including dynamic causal modeling [8, 9], functional network analysis [10, 11], its dynamic variants [12–14], and hidden Markov models [15–17]. Population-level inferences are a common practice for analyzing brain activity in empirical data. However, both the structure and dynamics of the brain vary even among healthy individuals, let alone among individuals belonging to a disease group due to the heterogeneity of the disease. Therefore, although population-level inferences increase the data size and often help us to reach statistically significant observations, they may yield inaccurate results and loss of information when the observed data are individual-specific. To avoid population-level inferences, it is necessary to establish the reliability of individual-level inferences of collective brain dynamics. Finn et al. examined the role of individual variability in functional networks measured by functional magnetic resonance imaging (fMRI) and its ability to act as a fingerprint to identify individuals [18] (see [19, 20] for earlier studies). In order for individual fingerprinting to be successful, the test-retest reliability of the functional network must be higher across sessions obtained from the same individual (i.e., within-participant reliability) than across sessions obtained from different individuals (i.e., between-participant reliability). Indeed, it was found that within-participant reliability was robust and that both resting-state and task fMRI from different sessions of the same individual could be used to perform fingerprinting [21]. Other studies also confirmed the ability of functional networks from fMRI data as fingerprints of individuals, including the development of different methods to quantify and improve fingerprinting [22–26]. The ability of functional connectivity to act as individual fingerprints has also been confirmed with electroencephalogram (EEG) [27] and magnetoencephalogram (MEG) data [28, 29]. Functional networks or its dynamic variants are not the only tools for analyzing brain dynamics or fingerprinting individuals. One way to analyze fMRI or other multidimensional time series data from the brain is to infer dynamics of discrete states. Each state may correspond to a particular functional network [13, 15, 30–33] or a spatial activation pattern [17, 34, 35], and the transition from one state to another may correspond to a regime shift in the brain. Energy landscape analysis is a method to characterize brain dynamics as a movement of a stochastic ball constrained on an energy landscape inferred from the data [36–38]. Quantifications of the estimated energy landscapes such as the height of the barrier between two local minima of the energy allow intuitive interpretations; a local minimum of the energy is a particular spatial activity pattern and defines a discrete brain state. A high barrier between two local minima implies that it is difficult for the brain dynamics to transit between the two local minima. Indices from energy landscape analysis have been shown to be associated with behavior of healthy individuals in a test of bistable visual perception task [37, 39], executive function [40], fluid intelligence [41], healthy aging [42], autism [43], Alzheimer disease [44], schizophrenia [45, 46], attention deficit hyperactivity disorder [47], and epilepsy [48]. These successful applications of energy landscape analysis are likely to owe to advantages of the method compared to other related methods such as functional network analysis and hidden Markov models. For example, with energy landscape analysis, one can borrow concepts and computational tools from statistical physics of spin systems to quantify the ease of state transition by the energy barrier [38] and complexity of the dynamics by different phases (e.g., spin-glass phase) and susceptibility indices [41]. In addition, each network state is by definition a binary activity pattern among a pre-specified set of regions of interest (ROIs) and therefore relatively easy to interpret. Despite its expanding applications, the validity of the energy landscape analysis has not been extensively studied except that one can measure the accuracy of fit of the model to the given data [38, 49–52]. A high accuracy of fit does not imply that the estimated energy landscape is a reliable fingerprint for individuals. In fact, if fMRI data are nonstationary, an energy landscape estimated for the same individual in two time windows may be substantially different from each other, whereas the accuracy of fit may be high in both time windows. Furthermore, the original energy landscape analysis method requires pooling of fMRI data from different individuals unless the number of regions of interest (ROIs) to be used is relatively small (e.g., 7) or the scanning session is extremely long. This is because the method is relatively data hungry [38]. The concept of individual fingerprinting is unclear when pooling of data is necessary. In the present study, we assess potential utility of energy landscape analysis in individual fingerprinting by investigating its test-retest reliability. Specifically, we ask how much features of the estimated energy landscapes are reproducible across different sessions from the same individual as opposed to across sessions belonging to different sets of individuals. We hypothesize that test-retest reliability is higher between sessions for the same individual than between sessions for different individuals. Code for computing energy landscapes with the conventional and Bayesian methods is available on Github [53]. 2. Methods 2.1. Midnight Scan Club data We primarily use the fMRI data in the Midnight Scan Club (MSC) data set [22]. MSC data set contains five hours of resting-state fMRI data in total recorded from each of the 10 healthy human adults across 10 consecutive nights. A resting-state fMRI scanning section lasted for 30 minutes and yielded 818 volumes. Imaging was performed with a Siemens TRIO 3T MRI scanner using an echo planar imaging (EPI) sequence (TR = 2.2 s, TE = 27 ms, flip angle = 90°, voxel size = 4 mm × 4 mm × 4 mm. 36 slices). The original paper reported that the eighth participant (i.e., MSC08) fell asleep, showed frequent and prolonged eye closures, and had systematically large head motion, resulting in much less reliable data than those obtained from the other participants [22]. We also noticed that the accuracy of fitting the energy landscape, which we will explain in section 2.4, fluctuated considerably across the different sessions for the tenth participant (i.e., MSC10), suggesting unstable quality of the MSC10’s data across sessions. Therefore, we excluded MSC08 and MSC10 from the analysis. We used SPM12 (http://www.fil.ion.ucl.ac.uk/spm) to preprocess the resting-state fMRI data as follows: we first conducted realignment, unwraping, slice-timing correction, and normalization to a standard template (ICBM 152); then, we performed regression analyses to remove the effects of head motion, white matter signals, and cerebrospinal fluid signals; finally, we conducted band-pass temporal filtering (0.01–0.1 Hz). We determined the ROIs of the whole-brain network using an atlas with 264 spherical ROIs whose coordinates were set in a previous study [54]. We then removed 50 ROIs labelled ‘uncertain’ or ‘subcortical’, which left us with 214 ROIs. The 214 ROIs were labeled either of the nine functionally different brain networks, i.e., auditory network, dorsal attention network (DAN), ventral attention network (VAN), cingulo-opercular network (CON), default mode network (DMN), fronto-parietal network (FPN), salience network (SAN), somatosensory and motor network (SMN), or visual network. We merged the DAN, VAN, and CON into an attention network (ATN) to reduce the number of observables from nine to seven, as we did in our previous study [43]. This is due to the relatively short length of the data and the fact that energy landscape analysis requires sufficiently long data sets if working with 9 observables. In fact, the DAN, VAN, and CON are considered to be responsible for similar attention-related cognitive activity [54], justifying the merge of the three systems into the ATN. We call the obtained N=7 dimensional time series of the fMRI signal the whole-brain network. We calculated the fMRI signal for each of the seven networks (i.e., ATN, auditory network, DMN, FPN, SAN, SMN, and visual network) by averaging the fMRI signal over the volumes in the sphere of radius 4 mm in the ROI and over all ROIs belonging to the network. In addition to the whole-brain network, we used a separate 30-ROI coordinate system [55] and determined the multi-ROI DMN and CON. We used a different parcellation system for the DMN and CON than the 264-ROI system used for the whole-brain network. It is because the former (i.e., 30-ROI) coordinate system provides much fewer ROIs for the DMN and CON than the 264-ROI system does, which is convenient for energy landscape analysis. The original study identified 12 and 7 ROIs for the DMN and CON, respectively [55]. To reduce the dimension of the DMN, we averaged over each pair of the symmetrically located right- and left-hemisphere ROIs in the DMN into one observable. The symmetrized DMN, which we simply call the DMN, has eight ROIs because four ROIs (i.e., amPFC, vmPFC, pCC, and retro splen) in the original DMN are almost on the midline and therefore have not undergone the averaging between the right- and left-hemisphere ROIs [42]. For the CON, we used the original seven ROIs as the observables. Note that the whole-brain network contains the DMN and CON as single observables, whereas the DMN and CON we are introducing here are themselves systems containing N=8 and N=7 observables, respectively. We denote the fMRI signal for the ith ROI at time t by xiti=1,…,N;t=1,…,tmax, where N is the number of ROIs, and tmax is the number of time points. We then removed the global signals and transformed the signals into their z-values using zit=xit-mt/st, where mt and st represent the mean and standard deviation, respectively, of xit over the N ROIs at time t;mt is the global signal [56]. The global signal in resting-state functional MRI data is considered to be dominated by physiological noise mainly originating from the respiratory, scanner-related, and motion-related artifacts. Global signal removal improves various quality-control metrics, enhances the anatomical specificity of functional-connectivity patterns, and can increase the behavioral variance [57, 58]. The same or similar global signal removal was carried out in previous energy landscape studies [41, 42]. 2.2. Human Connectome Project data For validation, we also analyzed fMRI data that were recorded from healthy human participants and shared as the S1200 data in the Human Connectome Project (HCP) [59]. In the data set, 1200 adults between 22–35 years old under-went four sessions of 15-min EPI sequence with a 3T Siemens Connectome-Skyra (TR = 0.72 s, TE = 33.1 ms, 72 slices, 2.0 mm isotropic, field of view (FOV) = 208 × 180 mm) and a T1-weighted sequence (TR = 2.4 s, TE = 2.14 ms, 0.7 mm isotropic, FOV = 224 × 224 mm). Here, we limited our analysis to those included in the 100 unrelated participant subset released by the HCP. We confirmed that all these 100 participants were among the subset of participants who completed diffusion weighted MRI as well as two resting-state fMRI scans. The resting-state fMRI data of each participant are composed of two sessions, and each session is broken down into a Left-Right (LR) and Right-Left (RL) phases. We used data from participants with at least 1150 volumes in each of the four sessions after removing volumes with motion artifacts, which left us with 87 participants. For the 87 participants, we first removed the volumes with motion artifacts. Then, we used the last 1150 volumes in each session to remove possible effects of transient. We used independent component analysis (ICA) to remove nuisance and motion signals [60]. Furthermore, any volumes with frame displacement greater than 0.2 mm [61] were excised [62] because the ICA-FIX pipeline has been found not to fully remove motion-related artifacts [63, 64]. We standardized each voxel by subtracting the temporal mean. Lastly, global signal regression of the same form as that for the MSC data (see section 2.1 was used for removing remaining noise. In each volume, we averaged the fMRI signal over all the voxels within each ROI of the AAL atlas [65]. Note that this atlas is composed of 116 ROIs. Then, we mapped each cortical ROI to either of the parcellation scheme from the Schaefer-100 atlas [66]. System assignment was based on minimizing the Euclidian distance from the centroid of an ROI in the AAL to the corresponding centroid of an ROI in the Schaefer atlas. We removed 42 ROIs labeled ‘subcortical’ or ‘cerebellar’, which left us with 74 ROIs. These 74 ROIs were labelled either of the N=7 functionally different brain networks: control network, DMN, DAN, limbic network, salience/ventral attention network, somatomotor network, and visual network, altogether defining a whole-brain network. 2.3. Fitting of the pairwise maximum entropy model To carry out energy landscape analysis, we fit the pairwise maximum entropy model (MEM), also known as the Ising model, to the preprocessed fMRI data in essentially the same manner as in previous studies [38, 67]. For each session, we first binarized zit for each ith ROI (with i∈{1,…,N}) and time t (with t∈1,…,tmax) using a threshold that we set to the time average of zit. A computational study showed that binarization did not affect important information contained in originally continuous brain signals [68]. We denote the binarized signal at the ith ROI and time t by σit, which is either +1 or −1 corresponding to whether zit is larger or smaller than the threshold, respectively. The activity pattern of the entire network at time t is described by the N-dimensional vector (1) Vt=σ1t,…,σNt∈{-1,1}N. It should be noted that there are 2N activity patterns in total, enumerated as V1,…,V2N. The empirical mean activity at the ith ROI is denoted by (2) σi≡1tmax∑t=1tmax σit. The empirical mean pairwise joint activation for the ith and jth ROIs is defined by (3) σiσj≡1tmax∑t=1tmax σitσjt. The pairwise MEM maximizes the entropy of the distribution of activity patterns under the condition that σi and σiσj (with 1≤i≤j≤N) are the same between the estimated model and the empirical data. The resulting probability distribution of activity pattern V=σ1,…,σN, denoted by P(V), obeys the Boltzmann distribution [69] given by (4) P(V)=e-E(V)∑k=12N  e-EVk, where E(V) represents the energy of activity pattern V given by (5) E(V)=-∑i=1N hiσi-12∑i=1N ∑j=1N Jijσiσj. In Eq. (5), the fitting parameter hi represents the tendency for the ith ROI to be active (i.e., σi=+1), and Jij quantifies the pairwise interaction between the ith and jth ROIs. We denote the mean activity and mean pairwise activity from the estimated model by σim and σiσjm, respectively. By definition, we obtain (6) σim=∑k=12N σiVkPVk and (7) σiσjm=∑k=12N σiVkσjVkPVk. We calculated hi and Jij by iteratively adjusting σim and σiσjm towards the empirically values, i.e., σi and σiσj, respectively, using a gradient ascent algorithm. The iteration scheme is given by (8) hinew=hiold +ϵlog⁡σiσim and (9) Jijnew =Jijold +ϵlog⁡σiσjσiσjm, where superscript new and old represent the values after and before a single updating step, respectively, and ϵ is the learning rate. We set ⁡ϵ=0.2. 2.4. Accuracy of fit We evaluated the accuracy of fit of the pairwise MEM to the given fMRI data [38, 42, 50]. The accuracy index is given by (10) rD=D1-D2D1, where (11) Dℓ=∑k=12N PNVklog2⁡PNVkPℓVk is the Kullback-Leibler divergence between the probability distribution of the activity pattern in the ℓth-order (ℓ=1,2) MEM, Pℓ(V), and the empirical probability distribution of the activity pattern, denoted by PN(V). Note that P2(V) is equivalent to P(V) given by Eqs. (4) and (5). The first-order, or independent, MEM (i.e., ℓ=1) is Eq. (4) without interaction terms, that is, Jij=0∀i,j in Eq. (5). We obtain rD=1 when the pairwise MEM perfectly fits the empirical distribution of the activity pattern, and rD=0 when the pairwise MEM does not fit the data any better than the independent MEM. To assess the dependency of rD on the number of sessions to be concatenated for the estimation of the pairwise MEM, m, the network (i.e., whole-brain, DMN, or CON), and the type of concatenation (i.e., within-participant or between-participant), we examined the multivariate linear regression model given by (12) rD=β0+β1m+β2Iwhole +β3ICON+β4Iwithin. In Eq. (12), β0 is the intercept, dummy variable Iwhole is equal to 1 for the whole-brain network and 0 for the other two networks, ICON is equal to 1 for the CON, and 0 for the other two networks, and Iwithin is equal to 1 for the within-participant comparison and 0 for the across-participant comparison. 2.5. Bayesian approximation method The pairwise MEM and the subsequent energy landscape analysis have mostly been restricted to analysis of group-level data. This is because the methods in its original form are data-hungry, requiring concatenation of fMRI signals from different individuals. The length of fMRI data, tmax, that is necessary for reliably estimating the pairwise MEM with N nodes is roughly proportional to the number of states, 2N [38]. To overcome this problem and obtain the energy landscape for each individual, we employed a recently developed variational Bayes approximation method for estimating the pairwise MEM [40, 70], which runs as follows. We denote by 𝒮n the N-dimensional time series obtained from an nth session of fMRI. Different fMRI sessions typically originate from different participants in the same group (e.g., control group). We denote the number of sessions available by D. Let 𝒮 be the concatenated data, i.e., (13) 𝒮≡∪n=1D𝒮n. The variational Bayes approximation method estimates a pairwise MEM for each 𝒮n (with n∈{1,…,D}). This method introduces a prior distribution for the set of session-specific model parameters, θn= h1,h2,…,hN,J12,J13,…,JN-1,N∈RM, where n∈{1,…,D} and M=N(N+1)/2. We give the prior distribution for (14) Θ=θ1,…,θD by (15) p(Θ∣η,α)=∏n=1D ∏M′=1M pθnM′∣𝒩ηM′,1/αM′, where px∣𝒩μ,σ2 represents the probability density of x obeying the one-dimensional normal distribution with mean and variance equal to μ and σ2, respectively. Here, η=η1,…,ηM⊤∈RM is the prior mean vector, α=α1,…,αM⊤∈R+M is the prior precision vector, and ⊤ represents the transposition. In Eq. (15), we have assumed that the signals from all the D sessions are mutually independent. Now, we derive the posterior distribution of Θ. It is intractable to derive the posterior because the normal distribution is not the conjugate prior for the Boltzmann distribution. Therefore, we use a variational approximation to the posterior [71] using the normal distribution as follows: (16) q(Θ∣𝒮,η,α)=∏n=1D ∏M′=1M pθnM′∣𝒩μnM′,1/βnM′. We write μn=μn1,…,μnM⊤∈RM and βn=βn1,…,βnM⊤∈R+M, which are the posterior mean vector and the posterior precision vector for session n∈{1,…,D}, respectively. One obtains the variational approximate solution for distribution q by optimizing the evidence lower bound (ELBO), also called the free energy [40, 70]. By maximizing the free energy with respect to q, we have the posterior mean and precision vectors in terms of the prior mean and precision vectors as follows: (17) μn=η+tmaxAη,α−1(〈σ¯n〉−〈σ¯〉η). (18) βn=α+tmaxcη, where (19) Aη,α=diag⁡(α)+tmaxCη, and diag(·) represents the diagonal matrix whose entries are given by the arguments. In Eq. (17), σ‾n≡σ1,…,σN,σ1σ2,σ1σ3,…,σN-1σN⊤ is the vector composed of the empirical mean activity and empirical pairwise joint activation; ⟨σ‾⟩η is the model mean of σ‾n≡σ1,σ2,…,σN,σ1σ2,σ1σ3,…,σN-1σN⊤ when the model parameters h1,h2,…,hN,J12,J13,…,JN-1,N are given by η;Cη≡Covη⁡σ‾n is the covariance matrix of σ‾n when the model is given by η. In Eq. (18), cη is the vector composed of the diagonal element of Cη. In other words, the ith element of cη is the variance of the ith element of σ‾n under parameters η. Now, we fix q and maximize the free energy with respect to η and α to obtain the equations for updating η and α as follows: (20) ηM′=1D∑n=1D  μnM′, (21) αM′=1D∑n=1D  μnM′-ηM′2+1βnM′-1, where M′∈{1,…,M}. Thus, we have updated the posterior distribution θnM′~𝒩μnM′,1/βnM′,n∈{1,…,D},M′∈{1,…,M} using the prior distribution θnM′~𝒩ηM′,1/αM′, and then updated the prior distribution using the new posterior distribution. We summarize the steps of the variational Bayes approximation method as follows: Initialize the hyperparameters by independently drawing each ηM′ (with M′∈{1,…,M}) from the normal distribution with mean 0 and standard deviation 0.1. We also set the first N entries of the prior precision vector αM′, corresponding to hi,i∈{1,…,N}, to 6, and set the remaining M-N entries of αM′ corresponding to Jij,1≤i