## Abstract

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’.

## 1. 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 1*g*). 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.

## 2. 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 1*a*) and obtain resting-state (or under other conditions, which are ideally stationary) fMRI signals at the ROIs, resulting in a multivariate time series (figure 1*b*). 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, *σ*_{i}(*t*)=1 () indicates that the *i*th ROI is active at time *t*, and *σ*_{i}(*t*)=−1 indicates that the ROI is inactive (figure 1*c*). 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 [12].

Third, we calculate the relative frequency with which each activity pattern is visited, *P*_{empirical}(** σ**) (figure 1

*d*). To

*P*

_{empirical}(

**), we fit the Boltzmann distribution given by 2.1where 2.2is the energy, and**

*σ***={**

*h**h*

_{i}} and

**J**={

*J*

_{ij}} (

*i*,

*j*=1,…,

*N*) are the parameters of the model (figure 1

*e*). We assume

*J*

_{ij}=

*J*

_{ji}and

*J*

_{ii}=0 (

*i*,

*j*=1,…,

*N*). The principle of maximum entropy imposes that we select

**and**

*h***J**such that 〈

*σ*

_{i}〉

_{empirical}=〈

*σ*

_{i}〉

_{model}and 〈

*σ*

_{i}

*σ*

_{j}〉

_{empirical}=〈

*σ*

_{i}

*σ*

_{j}〉

_{model}(

*i*,

*j*=1,…,

*N*), where 〈⋯ 〉

_{empirical}and 〈⋯ 〉

_{model}represent the mean with respect to the empirical distribution (figure 1

*d*) 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

*i*th ROI and the interaction between the

*i*th and

*j*th 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

*i*th 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
2.3where is the likelihood given by
2.4We can maximize the likelihood by a gradient ascent scheme
2.5and
2.6where 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 ** h** 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

^{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:
2.7where represents the Boltzmann distribution for a single spin, *σ*_{i}, given the other *σ*_{j}(*j*≠*i*) values fixed to *σ*_{/i}(*t*)≡(*σ*_{1}(*t*),…,*σ*_{i−1}(*t*),*σ*_{i+1}(*t*),…,*σ*_{N}(*t*)) [19]. In other words,
2.8We call the right-hand side of equation (2.7) the pseudo-likelihood. Although this is a mean-field 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 [19]. A gradient ascent updating scheme is given by
2.9and
2.10where and are the mean and correlation with respect to distribution (equation (2.8)) and are given by
2.11and
2.12respectively. 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 2.13where

*W*(

**|**

*σ***′) is a transition rate from activity pattern**

*σ***′ to activity pattern**

*σ***. As the initial condition, we impose**

*σ**P*(

**;0)=**

*σ**P*

_{empirical}(

**). By choosing 2.14where**

*σ***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 ** 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: 2.15where is the KL divergence, which quantifies the discrepancy between two distributions,

*Ω*is the set of all the 2

^{N}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

*D*

_{KL}(

*P*(

**;0) ∥**

*σ**P*(

**;**

*σ**ϵ*)) is practically used for determining

**and**

*h***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
2.16which ranges between 0 and 1 for the maximum-likelihood estimator, and is the Shannon entropy of the MEM incorporating correlations up to the *k*th 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
2.17Note 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 1*f*. 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 1

*f*, 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 2.18If 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^{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 and 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 2.19(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 1*g* as a hypothetical two-dimensional landscape. The local minimums and energy barriers in figure 1*g* correspond to those shown in the disconnectivity graph (figure 1*f*).

### (e) 0/1 versus 1/−1

We remark on two binarization schemes. In statistical physics, the pairwise MEM, or the Ising model, usually employs *σ*_{i}∈{−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 (*σ*_{i})^{2}=1 regardless of *σ*_{i} [24,25]. For neuronal spike data, 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 [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 *σ*_{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 is defined as 2.20Mathematically, the two representations are equivalent to the one-to-one relationship, , which results in 2.21and 2.22

## 3. 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 pseudo-likelihood maximization in figure 2*a*–*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 2*g*–*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. , 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 ** h**. 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,

*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

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^{N} 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 ). If 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 .**

*σ*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 maximum-likelihood estimator). For a given ℓ, we calculated *r* for each of the *t*_{max}−ℓ datasets of length ℓ, i.e. {** σ**(1),…,

**(ℓ)}, {**

*σ***(2),…,**

*σ***(ℓ+1)}, …, . The average and standard deviation of**

*σ**r*as a function of ℓ/2

^{N}are shown in figure 4

*a*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 into two halves of length , four quarters of length , eight non-overlapping samples of length and so forth. The results (figure 4*b*) were similar to those in the case of overlapping time windows (figure 4*a*).

## 4. 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, *t*_{max}, 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.

## 5. 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 () 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.

### (b) 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 *t*_{max}=9560 at 30 ROIs (12 in the DMN, 11 in the FPN and 7 in the CON).

## Authors' contributions

T.E. carried out the numerical experiments and performed the data analysis. T.W. preprocessed the fMRI data. N.M. conceived of and designed the study. T.E., T.W. and N.M. drafted the manuscript. M.O. provided algorithms and discussed theoretical aspects of the manuscript. All authors read and approved the manuscript.

## Competing interests

The authors declare that they have no competing interests.

## Funding

T.E. acknowledges the support provided through PRESTO, JST. T.W. acknowledges the support provided through the European Commission. M.O. acknowledges the support provided through JSPS KAKENHI grant no. 15H03699. N.M. acknowledges the support provided through CREST, JST and the Kawarabayashi Large Graph Project, ERATO, JST and EPSRC Institutional Sponsorship to the University of Bristol.

## Acknowledgements

Data were provided in part by the Human Connectome Project, WU-Minn Consortium (principal investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657), funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

## Footnotes

One contribution of 14 to a theme issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3726625.

- Accepted January 27, 2017.

- © 2017 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.