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][2][3][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][6][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][9][10][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][9][10][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 Dijkstralike 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.

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

(a) 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, {σ i (1), . . . , σ i (t max )}, where t max is the length of the data, σ i (t) = 1 (t = 1, . . . , t max ) indicates that the ith ROI is active at time t, and σ i (t) = −1 indicates that the ROI is inactive (figure 1c). The threshold is arbitrary, and we set it to the time average of σ i (t) for each i. The activity pattern of the entire network at time t is given by an N-dimensional vector σ ≡ (σ 1 , . . . , σ N ) ∈ {−1, 1} N , where we have suppressed t. Note that there are 2 N 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 [   Third, we calculate the relative frequency with which each activity pattern is visited, P empirical (σ ) (figure 1d). To P empirical (σ ), we fit the Boltzmann distribution given by where is the energy, and h = {h i } and J = {J ij } (i, j = 1, . . . , N) are the parameters of the model (figure 1e). We assume J ij = J ji and J ii = 0 (i, j = 1, . . . , N). The principle of maximum entropy imposes that we select h and J such that σ i empirical = σ i model and 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 i and J ij 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 i is large, the energy is smaller with σ i = 1 than with σ i = −1, such that the ith ROI tends to be active.

(b) Algorithms to estimate the pairwise maximum entropy model
In this section, we review three algorithms to estimate the parameters of the MEM, i.e. h and J.

(i) Likelihood maximization
In the maximum-likelihood method, we solve We can maximize the likelihood by a gradient ascent scheme 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 sgn( σ i ) log( σ i empirical / σ i model ) and sgn( σ i σ j ) log( σ i σ j empirical / σ i σ j model ), respectively, is also used sometimes [5,12,17,18]. Because equation (2.1) is concave in terms of h and J (which we can show by verifying that the Hessian of log L 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 N activity patterns to calculate σ i model and σ i σ j 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.

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

(iii) 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 N 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) = P empirical (σ ). By choosing In the minimum probability flow method, we look for h 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) (=P empirical (σ )), 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 D KL (P(σ ; 0) P(σ ; )) = σ ∈Ω P(σ ; 0) log[P(σ ; 0)/P(σ ; )] is the KL divergence, which quantifies the discrepancy between two distributions, Ω is the set of all the 2 N activity patterns and D = {σ ∈ Ω | P empirical (σ ) > 0} 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 t max is large such that most activity patterns appear in the data. However, when N is large or t max is small, this algorithm is efficient in terms of the computation time and memory space [20]. A gradient descent method on D KL (P(σ ; 0) P(σ ; )) is practically used for determining h and J.

(c) Accuracy indices
Fully describing an empirical distribution requires 2 N − 1 parameters, whereas the pairwise MEM only uses N + N(N − 1)/2 parameters. The pairwise MEM imposes that the first two moments of σ i 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 S k ≡ − σ ∈Ω P k (σ ) log P k (σ ) 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 ij = 0 for i, j = 1, . . . , N), gives P 1 (σ ). The pairwise MEM gives P 2 (σ ). The empirical distribution (i.e. P empirical (σ )) is identical to P N (σ ). The denominator of equation (2.16), S 1 − S N ≡ I N , 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, S 1 − S 2 ≡ I 2 , is equal to the contribution of the pairwise correlation. If I 2 /I N = 1, the pairwise correlation alone accounts for all the correlations present in the empirical data. If I 2 /I N = 0, the pairwise correlation does not deliver any information.
The second measure is defined by r = D KL (P 1 (σ ) P N (σ )) − D KL (P 2 (σ ) P N (σ )) D KL (P 1 (σ ) P N (σ )) . (2.17) 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 I 2 /I N = r [5,18].

(d) 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 E max (α ↔ α ) 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 E max (α ↔ α ) to travel between α and α . Because a large energy value corresponds to a low frequency of the activity pattern, a large E max (α ↔ α ) 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Ē αα − E(α).
To calculateĒ αα , we employ a Dijkstra-like method as follows. Consider the hypercube composed of 2 N 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Ē αα = E(α) andĒ αα = E(α ) for all α that are neighbours of α. These values are finalized and will not be changed.Ē αα for the other 2 N − 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 ).  1} (i = 1, . . . , N) rather thanσ i ∈ {0, 1}. 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 (σ i ) 2 = 1 regardless of σ i [24,25]. For neuronal spike data,σ i ∈ {0, 1} is often used [22,26,27], whereas σ i ∈ {−1, 1} is also commonly used [17,18,28,29]. For fMRI data, our previous work employedσ i ∈ {0, 1} [12][13][14][15]. The use ofσ i ∈ {0, 1} in representing neuronal spike trains has a rationale in being able to express the instantaneous firing rate in a simple form N i=1 σ i 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 σ i ∈ {−1, 1} because σ i empirical and σ i model , the logarithm of whose ratio is used in the algorithm, may have opposite signs.
The energy in the case ofσ i ∈ {0, 1} is defined as Mathematically, the two representations are equivalent to the one-to-one relationship, 2σ i − 1 = σ i (i = 1, . . . , N), which results inh

Results (a) 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, N ROI = 12), fronto-parietal network (FPN, N ROI = 11) and cingulo-opercular network (CON, N ROI = 7), using the ROIs whose coordinates were identified previously [31]. We had t max = 9560 volumes in total.
The estimated parameter values are compared between likelihood maximization and pseudolikelihood 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 2 8 activity patterns appeared at least once, i.e. D = Ω, which made the right-hand side of equation (2.15) zero. The     substantially smaller than the values for the likelihood or pseudo-likelihood maximization. By contrast, I 2 /I N exceeded unity because S 1 > S N > S 2 for the minimum probability flow method.
(b) 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.
(c) Effects of the data length  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, P empirical (σ ), would not be reliable because it is evaluated only based on a few visits to σ (divided by t max ). If t max is much larger and σ is visited proportionally many times, then we would be able to estimate P empirical (σ ) more accurately. This exercise led us to hypothesize that the accuracy scales as a function of t max /2 N . To examine this point, we carried out likelihood maximization on the fMRI data of varying length (t max /20 ≤ ≤ t max ) and calculated r (which coincides with I 2 /I N for the maximumlikelihood estimator). For a given , we calculated r for each of the t max − datasets of length , i.e. {σ (1), . . . , σ ( )}, {σ (2), . . . , σ ( + 1)}, . . ., {σ (t max − + 1), . . . , σ (t max )}. The average and standard deviation of r as a function of /2 N 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.
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 {σ (1), . . . , σ (t max )} into two halves of length = t max /2, four quarters of length = t max /4, eight non-overlapping samples of length = t max /8 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][13][14][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. 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 t max = 75-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 t max = 750-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][2][3][4]. Altered dynamics in the brain at various spatial and time scales may result in deformation of energy landscapes as compared with healthy controls.

Material and methods (a) 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 (≈15 min run −1 ) 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.