Energy landscape analysis of neuroimaging data.

Computational neuroscience models have been used for understanding neural dynamics in the brain and how they may be altered when physiological or other conditions change. We review and develop a data-driven approach to neuroimaging data called the energy landscape analysis. The methods are rooted in statistical physics theory, in particular the Ising model, also known as the (pairwise) maximum entropy model and Boltzmann machine. The methods have been applied to fitting electrophysiological data in neuroscience for a decade, but their use in neuroimaging data is still in its infancy. We first review the methods and discuss some algorithms and technical aspects. Then, we apply the methods to functional magnetic resonance imaging data recorded from healthy individuals to inspect the relationship between the accuracy of fitting, the size of the brain system to be analysed and the data length.This article is part of the themed issue 'Mathematical methods in medicine: neuroscience, cardiology and pathology'.

Introduction

Altered cognitive function due to various neuropsychiatric disorders seems to result from aberrant neural dynamics in the affected brain [1-4]. Alterations in brain dynamics may also occur in the absence of disorders, in situations such as typical ageing, traumatic experiences, emotional responses and tasks. Functional magnetic resonance imaging (fMRI) provides information on the neural dynamics in the brain with a reasonable spatial resolution in a non-invasive manner. There are various analysis methods that can be used to extract the dynamics in neuroimaging data including fMRI signals, such as sliding-window functional connectivity analysis, dynamic causal modelling, oscillation analysis and biophysical modelling. In this study, we seek the potential of a different approach: energy landscape analysis.

This method is rooted in statistical physics. The main idea is to map the brain dynamics to the movement of a ‘ball’ constrained on an energy landscape inferred from neural data. A ball tends to go downhill and remains near the bottom of a basin in a landscape, whereas it sometimes goes uphill due to random fluctuations that cause it to wander around and possibly transit to another basin (figure 1g). Using the Ising model (equivalent to the Boltzmann machine and the pairwise maximum entropy model (MEM); see [5-7] for reviews in neuroscience), we can explicitly construct an energy landscape from multivariate time-series data including fMRI signals recorded at a specified set of regions of interest (ROIs). The pairwise MEM, or, equivalently, the Ising model, has been used to emulate fMRI signals [6,8-11]. More recently, we have used the pairwise MEM for fMRI data during rest [12] and sleep [13] and then developed an energy landscape analysis method and applied it to participants during rest [14] and during a bistable visual perception task [15]. In contrast with the aforementioned previous studies [6,8-11], our approach is data driven, with the parameters of the Ising model being inferred from the given data. In this paper, we review the methods and some technical details. In passing, we introduce new techniques (i.e. different inference algorithms and a Dijkstra-like method). We apply the methods to publicly shared resting-state fMRI data recorded from healthy human participants to validate the new approaches and also to examine the relationship between the accuracy of fit, the size of the brain system (i.e. number of ROIs) and the length of the fMRI data.

Figure 1.

Procedures of the energy landscape analysis for fMRI data. (a) ROIs are specified. (b) fMRI signals are measured at the ROIs. (c) The fMRI signal at each ROI and each time point is binarized into 1 (active) or −1 (inactive). (d) The normalized frequency is computed for each binarized activity pattern. (e) The pairwise MEM model (i.e. Boltzmann distribution) is fitted to the empirical distribution of the 2 activity patterns (equation (2.1)). The energy value is also obtained for each activity pattern (equation (2.2)). (f) Relationships between activity patterns that are energy local minimums are summarized into a disconnectivity graph. (g) Schematic of the energy landscape. Each local minimum corresponds to the bottom of a basin. The borders between attractive basins of different local minimums are shown by the dotted curves. Any activity pattern belongs to the basin of a local minimum. Brain dynamics can be interpreted as the motion of a ‘ball’ constrained on the energy landscape. (Online version in colour.)

Material and methods

The pipeline of the energy landscape analysis based on the pairwise MEM is illustrated in figure 1.

Pairwise maximum entropy model

First, we specify a brain network of interest (figure 1a) and obtain resting-state (or under other conditions, which are ideally stationary) fMRI signals at the ROIs, resulting in a multivariate time series (figure 1b). We denote the number of ROIs by N.

Second, we binarize the fMRI signal at each time point (i.e. in each image volume) and each ROI by thresholding the signal. Then, for each ROI i (i=1,…,N), we obtain a sequence of binarized signals representing the brain activity, , where is the length of the data, σ(t)=1 () indicates that the ith ROI is active at time t, and σ(t)=−1 indicates that the ROI is inactive (figure 1c). The threshold is arbitrary, and we set it to the time average of σ(t) for each i. The activity pattern of the entire network at time t is given by an N-dimensional vector ≡(σ1,…,σ)∈{−1,1}, where we have suppressed t. Note that there are 2 possible activity patterns in total. Binarization is not readily justified given continuously distributed fMRI signals. However, we previously showed that the pairwise MEM with binarized signals predicted anatomical connectivity of the brain better than other functional connectivity methods that are based on non-binarized continuous fMRI signals and that ternary as opposed to binary quantization did not help to improve the results [12].

Third, we calculate the relative frequency with which each activity pattern is visited, Pempirical() (figure 1d). To Pempirical(), we fit the Boltzmann distribution given by where is the energy, and ={h} and J={J} (i,j=1,…,N) are the parameters of the model (figure 1e). We assume J=J and J=0 (i,j=1,…,N). The principle of maximum entropy imposes that we select and J such that 〈σ〉empirical=〈σ〉model and 〈σσ〉empirical=〈σσ〉model (i,j=1,…,N), where 〈⋯ 〉empirical and 〈⋯ 〉model represent the mean with respect to the empirical distribution (figure 1d) and the model distribution (equation (2.1)), respectively. By maximizing the entropy of the distribution under these constraints, we obtain the Boltzmann distribution given by equation (2.1). Some algorithms for the fitting will be explained in §2b. Equation (2.1) indicates that an activity pattern with a high energy value is not frequently visited and vice versa. Values of h and J represent the baseline activity at the ith ROI and the interaction between the ith and jth ROIs, respectively. Equation (2.2) implies that, if h is large, the energy is smaller with σ=1 than with σ=−1, such that the ith ROI tends to be active.

Algorithms to estimate the pairwise maximum entropy model

In this section, we review three algorithms to estimate the parameters of the MEM, i.e. and J.

Likelihood maximization

In the maximum-likelihood method, we solve where is the likelihood given by We can maximize the likelihood by a gradient ascent scheme and where the superscripts new and old represent the values after and before a single updating step, respectively, and ϵ (>0) is a constant. A slightly different updating scheme called the iterative scaling algorithm [16], where the right-hand side of equations (2.5) and (2.6) is replaced by and , respectively, is also used sometimes [5,12,17,18]. Because equation (2.1) is concave in terms of and J (which we can show by verifying that the Hessian of is a type of sign-flipped covariance matrix, which is negative semi-definite), the gradient ascent scheme yields the maximum-likelihood estimator. Because a single updating step involves all the 2 activity patterns to calculate 〈σ〉model and 〈σσ〉model, likelihood maximization is computationally costly for large N.

Our Matlab code to calculate the maximum-likelihood estimator for arbitrary multivariate time-series data is available in the electronic supplementary material.

Pseudo-likelihood maximization

The pseudo-likelihood maximization method approximates the likelihood function as follows: where represents the Boltzmann distribution for a single spin, σ, given the other σ(j≠i) values fixed to /(t)≡(σ1(t),…,σ(t),σ(t),…,σ(t)) [19]. In other words, We call the right-hand side of equation (2.7) the pseudo-likelihood. Although this is a mean-field approximation neglecting the influence of σ on σ (j≠i), the estimator obtained by the maximization of the pseudo-likelihood approaches the maximum-likelihood estimator as [19]. A gradient ascent updating scheme is given by and where and are the mean and correlation with respect to distribution (equation (2.8)) and are given by and respectively. It should be noted that this updating rule circumvents the calculation of 〈σ〉model and 〈σσ〉model, which the gradient ascent method to maximize the original likelihood uses and involves 2 terms.

Minimum probability flow

Different from the likelihood and pseudo-likelihood maximization, the minimum probability flow method [20] is not based on the likelihood function. Consider relaxation dynamics of a probability distribution, P(;τ), on the 2 activity patterns whose master equation is given by where W(|′) is a transition rate from activity pattern ′ to activity pattern . As the initial condition, we impose P(;0)=Pempirical(). By choosing where and ′ are neighbours if they are only different at one ROI, we obtain a standard Markov chain Monte Carlo method such that P(;τ) converges to the Boltzmann distribution given by equation (2.1).

In the minimum probability flow method, we look for and J values for which P(;τ) changes little in the relaxation dynamics at τ=0 [20]. The idea is that only a small amount of relaxation is necessary if the initial distribution, i.e. P(;0) (=Pempirical()), is sufficiently close to the equilibrium distribution, i.e. the Boltzmann distribution. The Kullback–Leibler (KL) divergence between the empirical distribution, P(;0), and a probability distribution after an infinitely small relaxation time, P(;ϵ), is approximated as: where is the KL divergence, which quantifies the discrepancy between two distributions, Ω is the set of all the 2 activity patterns and is a set of activity patterns that appear at least once in the empirical data. The minimum probability flow method minimizes the last quantity in equation (2.15), i.e. the probability flow from activity patterns that appear in the data to the patterns that do not. Therefore, the method is not effective when N is small and is large such that most activity patterns appear in the data. However, when N is large or is small, this algorithm is efficient in terms of the computation time and memory space [20]. A gradient descent method on DKL(P(;0) ∥ P(;ϵ)) is practically used for determining and J.

Accuracy indices

Fully describing an empirical distribution requires 2−1 parameters, whereas the pairwise MEM only uses N+N(N−1)/2 parameters. The pairwise MEM imposes that the first two moments of σ agree between the empirical data and the model. However, the model may be inaccurate in describing higher-order correlations in the empirical data. Most previous studies used one or both of the following two measures to quantify the accuracy with which the MEM fitted the empirical data.

The first measure is defined by which ranges between 0 and 1 for the maximum-likelihood estimator, and is the Shannon entropy of the MEM incorporating correlations up to the kth order [17,21]. The so-called independent MEM, in which we suppress any interaction between the elements (i.e. J=0 for i,j=1,…,N), gives P1(). The pairwise MEM gives P2(). The empirical distribution (i.e. Pempirical()) is identical to P(). The denominator of equation (2.16), S1−S≡I, is referred to as the multi-information, which quantifies the total contribution of the second or higher order correlation to the entropy of the empirical distribution. The numerator, S1−S2≡I2, is equal to the contribution of the pairwise correlation. If I2/I=1, the pairwise correlation alone accounts for all the correlations present in the empirical data. If I2/I=0, the pairwise correlation does not deliver any information.

The second measure is defined by Note that r also ranges between 0 and 1 for the maximum-likelihood estimator [5,12,22]. If the pairwise MEM produces a distribution closer to the empirical distribution than the independent MEM does, r is large. If the pairwise MEM and the independent MEM are similar in terms of the proximity to the empirical distribution, we obtain r≈0. For the maximum-likelihood estimator, we obtain I2/I=r [5,18].

Disconnectivity graph and energy landscape

Once we have estimated the pairwise MEM, we construct a dendrogram referred to as a disconnectivity graph [23], as shown in figure 1f. In the disconnectivity graph, a leaf (with a loose end open downwards) corresponds to an activity pattern that is a local minimum of the energy, i.e. an activity pattern whose frequency is higher than any other activity pattern in the neighbourhood of . The neighbourhood of is defined as the set of the N activity patterns that are different from only at one ROI. In the disconnectivity graph, the vertical position of the endpoint of the leaf represents its energy value, which specifies the frequency of appearance, with a lower position corresponding to a higher frequency. The branching structure of the disconnectivity graph describes the energy barrier between any pair of activity patterns that are local minimums. For example, to transit from local minimum 1 to local minimum 3 in figure 1f, the brain dynamics have to overcome the height of the energy barrier (shown by the double-headed arrow). If the barrier is high, transitions between the two activity patterns occur with a low frequency.

The disconnectivity graph is obtained by the following procedures. First, we enumerate local minimums, i.e. the activity patterns whose energy is smaller than that of all neighbours. Then, for a given pair of local minimums and ′, we consider a path connecting them, , where a path is a sequence of activity patterns starting from and ending at such that any two consecutive activity patterns on the path are neighbouring patterns. We denote by the largest energy value among the activity patterns on path . The brain dynamics on this path must climb up the hill to go through the activity pattern with energy to travel between and . Because a large energy value corresponds to a low frequency of the activity pattern, a large value implies that the frequency of switching between and along this path is low. Because various paths may connect and , we consider If we remove all the rarest activity patterns whose energy is equal to or larger than , and are disconnected (i.e. no path connecting them exists). The energy barrier for the transition from to is given by .

To calculate , we employ a Dijkstra-like method as follows. Consider the hypercube composed of 2 activity patterns. By definition, two nodes (i.e. activity patterns) are adjacent to each other (i.e. directly connected by a link) if they are neighbouring activity patterns. Each node has degree (i.e. number of neighbours) N. Then, fix a local minimum activity pattern and look for for all local minimums . We set and for all that are neighbours of . These values are finalized and will not be changed. for the other 2−N−1 local minimums are initialized to . Then, we iterate the following procedures until values for all the nodes are finalized. (i) For each finalized , update E for its all unfinalized neighbours using (ii) Find with the smallest unfinalized value and finalize it. (iii) Repeat steps (i) and (ii). If we carry out the entire procedure for each local minimum , we obtain for all pairs of local minimums.

By collecting pairs of local minimums that have the same value, we specify a set of local minimums that should be located under the same branch. This information is sufficient for drawing the dendrogram of local minimums, i.e. the disconnectivity graph.

Each local minimum has a basin of attraction in the state space, Ω. Each activity pattern, denoted by , usually belongs to one of the attractive basins, which is determined as follows. (i) Unless is a local minimum, move to the neighbouring activity pattern that has the smallest energy value. (ii) Repeat step (i) until a local minimum, denoted by , is reached. We conclude that belongs to the attractive basin of . (iii) Repeat steps (i) and (ii) for all the initial activity patterns ∈Ω.

Using the information on the local minimums and attractive basins, the dynamics of the activity pattern are illustrated as the motion of a ‘ball’ on the energy landscape, as schematically shown in figure 1g as a hypothetical two-dimensional landscape. The local minimums and energy barriers in figure 1g correspond to those shown in the disconnectivity graph (figure 1f).

0/1 versus 1/−1

We remark on two binarization schemes. In statistical physics, the pairwise MEM, or the Ising model, usually employs σ∈{−1,1} (i=1,…,N) rather than . The former convention respects the symmetry between the two spin states and is also convenient in some analytical calculations of the model that exploit the relationship (σ)2=1 regardless of σ [24,25]. For neuronal spike data, is often used [22,26,27], whereas σ∈{−1,1} is also commonly used [17,18,28,29]. For fMRI data, our previous work employed [12-15]. The use of in representing neuronal spike trains has a rationale in being able to express the instantaneous firing rate in a simple form and synchronous firing of neurons by a simple multiplication [30]. For example, three neurons simultaneously fire if and only if σ1σ2σ3=1. It should also be noted that the iterative scaling algorithm for maximizing the likelihood (§2b(i)) does not generally work for σ∈{−1,1} because 〈σ〉empirical and 〈σ〉model, the logarithm of whose ratio is used in the algorithm, may have opposite signs.

The energy in the case of is defined as Mathematically, the two representations are equivalent to the one-to-one relationship, , which results in and

Results

Accuracy of the three methods

We applied the three methods to estimate the pairwise MEM to resting-state fMRI signals recorded from two healthy adult individuals in the Human Connectome Project. We extracted ROIs from three brain systems, i.e. default mode network (DMN, NROI=12), fronto-parietal network (FPN, NROI=11) and cingulo-opercular network (CON, NROI=7), using the ROIs whose coordinates were identified previously [31]. We had tmax=9560 volumes in total.

The estimated parameter values are compared between likelihood maximization and pseudo-likelihood maximization in figure 2a–f. For all the networks, the results obtained by the pseudo-likelihood maximization are close to those obtained by the likelihood maximization, in particular for J. The results obtained by the likelihood maximization and those obtained by the minimum probability flow are compared in figure 2g–j for the DMN and FPN. We did not apply the minimum probability flow method to the CON because all of the 28 activity patterns appeared at least once, i.e. , which made the right-hand side of equation (2.15) zero. The figure indicates that the estimation by the minimum probability flow deviates from that by the likelihood maximization more than the estimation by the pseudo-likelihood maximization does, in particular for . The two measures of the accuracy indices are shown in table 1 for each network and estimation method. The two indices took the same value in the case of the likelihood maximization [5,18]. In the case of the pseudo-likelihood maximization, the two accuracy indices were slightly different from each other, and both took approximately the same values as those derived from the maximum likelihood. In the case of the minimum probability flow, r was substantially smaller than the values for the likelihood or pseudo-likelihood maximization. By contrast, I2/I exceeded unity because S1>S>S2 for the minimum probability flow method.

Figure 2.

Estimated parameter values compared between the different methods. (a–f) Comparison between the likelihood maximization and the pseudo-likelihood estimation. (g–j) Comparison between the likelihood maximization and the minimum probability flow method. (a,d,g,i) DMN, (b,e,h,j) FPN and (c,f) CON. Superscripts L, PL and MP stand for likelihood maximization, pseudo-likelihood maximization and the minimum probability flow, respectively. (Online version in colour.)

Table 1.

Accuracy of fitting for each network and estimation algorithm.

	DMN		FPN		CON
	r	I2/IN	r	I2/IN	r	I2/IN
likelihood maximization	0.6921	0.6921	0.7830	0.7830	0.9744	0.9744
pseudo-likelihood maximization	0.6921	0.6972	0.7830	0.7853	0.9745	0.9744
minimum probability flow	0.6480	0.7437	0.6124	1.2295	—	—

Disconnectivity graphs

Figure 3 shows the disconnectivity graph of the DMN, FPN and CON, calculated for the parameter values estimated by likelihood maximization. The two synchronized activity patterns, i.e. the activity patterns with all ROIs being active or inactive, were local minimums. The FPN had much more local minimums than the DMN and CON did. Although the present results are opposite to our previous results using a different dataset [14], the reason for the discrepancy is unclear.

Figure 3.

Disconnectivity graphs for (a) DMN, (b) FPN and (c) CON. The activity pattern at each local minimum is also shown. Retro splen: retro splenial cortex, latP: lateral parietal cortex, pCC: posterior cingulate cortex, parahippo: parahippocampal cortex, inf templ: inferior temporal cortex, sup frontal: superior frontal cortex, vmPFC: ventromedial prefrontal cortex, amPFC: anteromedial prefrontal cortex, precun: precuneus, IPS: intraparietal sulcus, IPL: inferior parietal lobule, mCC: mid-cingulate cortex, frontal: lateral frontal cortex, dlPFC: dorsolateral prefrontal cortex, ant thal: anterior thalamus, dACC/msFC: dorsal anterior cingulate cortex/medial superior frontal cortex, al/fO: anterior insula/frontal operculum, aPFC: anterior prefrontal cortex. (Online version in colour.)

Effects of the data length

Our experiences suggest that, as the number of ROIs, N, increases, the pairwise MEM seems to demand a large amount of data to realize a high accuracy. If we use N ROIs, there are 2 possible activity patterns. Therefore, as we increase N, it is progressively more likely that many of the activity patterns are unvisited. However, the MEM assigns a positive probability to such an unvisited pattern. Even if an activity pattern is realized by the empirical data a few times, the empirical distribution, Pempirical(), would not be reliable because it is evaluated only based on a few visits to (divided by ). If is much larger and is visited proportionally many times, then we would be able to estimate Pempirical() more accurately. This exercise led us to hypothesize that the accuracy scales as a function of .

To examine this point, we carried out likelihood maximization on the fMRI data of varying length ℓ (tmax/20≤ℓ≤tmax) and calculated r (which coincides with I2/I for the maximum-likelihood estimator). For a given ℓ, we calculated r for each of the tmax−ℓ datasets of length ℓ, i.e. {(1),…,(ℓ)}, {(2),…,(ℓ+1)}, …, . The average and standard deviation of r as a function of ℓ/2 are shown in figure 4a for the DMN, FPN and CON. As expected, the accuracy improved as ℓ increased. The results for the DMN, FPN and CON roughly collapsed on a single curve. The figure suggests that, to achieve an accuracy of 0.8 and 0.9, each activity pattern should be visited ≈5 and ≈16 times on average, respectively.

Figure 4.

Accuracy of the fit of the pairwise MEM as a function of the data length, ℓ, normalized by the number of possible activity patterns, 2. We obtained samples of length ℓ by (a) overlapping sliding time windows and (b) non-overlapping time windows. Error bars represent the standard deviation. (Online version in colour.)

Because the aforementioned sampling method used overlapping time windows to make different samples strongly depend on each other, we carried out the same test by dividing the entire time series into two halves of length , four quarters of length , eight non-overlapping samples of length and so forth. The results (figure 4b) were similar to those in the case of overlapping time windows (figure 4a).

Discussion

We explained a set of computational methods to estimate the pairwise MEM and energy landscapes from resting-state fMRI data. Novel components, as compared with our previous methods [12-15], were the pseudo-likelihood maximization, the minimum probability flow and a variant of the Dijkstra method to calculate the disconnectivity graph. We applied the methods to fMRI data collected from healthy participants and assessed the amount of data needed to secure a sufficient accuracy of fit.

The present results suggest that the current method is admittedly demanding in terms of the amount of data, although the results should be corroborated with different datasets. In the application of the pairwise MEM to neuronal spike trains, the data length does not seem to pose a severe limit if the network size, N, is comparable to those in this study. This is because one typically uses a high time resolution to ensure that there are no multiple spikes within a time window (e.g. 2 ms [28], 10 ms [22], 20 ms [17,18,27]). Then, the number of data points, tmax, is typically much larger than in typical fMRI experiments. In fMRI experiments, the interval between two measurements, called the repetition time (TR), is typically 2–4 s, and a participant in the resting state (or a particular task condition) can be typically scanned for 5–15 min. Then, we would have –450, with which we can reliably estimate the pairwise MEM model up to N≈5 (figure 4), which is small. If we pool fMRI data from 10 participants belonging to the same group to estimate one MEM, we would have –4500, accommodating N≈8. This is an important limitation of our approach. Currently we cannot apply the method to relatively large brain systems (i.e. those with a larger number of ROIs), let alone voxel-based data.

We demonstrated the methods with fMRI data obtained from healthy participants. The same methods can be applied to different conditions of human participants including the case of medical applications, the topic of the present theme issue. Various neuropsychiatric disorders have been suggested to have dynamical footprints in the brain [1-4]. Altered dynamics in the brain at various spatial and time scales may result in deformation of energy landscapes as compared with healthy controls.

Data and participants

We used resting-state fMRI data publicly shared in the Human Connectome Project (acquisition Q10 in release S900 of the WU-Minn HCP data) [32]. The data were collected using a 3T MRI (Skyra, Siemens) with an echo planar imaging (EPI) sequence (TR, 0.72 s; TE, 33.1 ms; 72 slices; 2.0 mm isotopic; field of view, 208×180 mm) and T1-weighted sequence (TR, 2.4 s; TE, 2.14 ms; 0.7 mm isotopic; field of view, 224×224 mm). The EPI data were recorded in four runs () while participants were instructed to relax while looking at a fixed cross mark on a dark screen.

We used such EPI and T1 images recorded from two adult participants (one female; 22–25 years old), because the amount of the data was sufficiently large for the current analysis.

Preprocessing and extraction of region of interest data

We preprocessed the EPI data in essentially the same manner as the conventional methods that we previously used for resting-state fMRI data [12,33,34] with SPM12 (www.fil.ucl.ac.uk/spm). Briefly, after discarding the first five images in each run, we conducted realignment, unwarping, slice timing correction, normalization to the standard template (ICBM 152) and spatial smoothing (full-width at half maximum =8 mm). Afterwards, we removed the effects of head motion, white matter signals and cerebrospinal fluid signals by a general linear model. Finally, we performed temporal band-pass filtering (0.01–0.1 Hz) and obtained resting-state whole-brain data.

We then extracted a time series of fMRI signals from each ROI. The ROIs were defined as 4 mm spheres around their centre whose coordinates were determined in a previous study [31]. The signals at each ROI were those averaged over the sphere. In total, we obtained time-series data of length tmax=9560 at 30 ROIs (12 in the DMN, 11 in the FPN and 7 in the CON).