A pairwise maximum entropy model accurately describes resting-state human brain networks

Abstract

The resting-state human brain networks underlie fundamental cognitive functions and consist of complex interactions among brain regions. However, the level of complexity of the resting-state networks has not been quantified, which has prevented comprehensive descriptions of the brain activity as an integrative system. Here, we address this issue by demonstrating that a pairwise maximum entropy model, which takes into account region-specific activity rates and pairwise interactions, can be robustly and accurately fitted to resting-state human brain activities obtained by functional magnetic resonance imaging. Furthermore, to validate the approximation of the resting-state networks by the pairwise maximum entropy model, we show that the functional interactions estimated by the pairwise maximum entropy model reflect anatomical connexions more accurately than the conventional functional connectivity method. These findings indicate that a relatively simple statistical model not only captures the structure of the resting-state networks but also provides a possible method to derive physiological information about various large-scale brain networks.

ARTICLE

Received 15 May 2012 | Accepted 14 Dec 2012 | Published 22 Jan 2013

Takamitsu Watanabe1, Satoshi Hirose1, Hiroyuki Wada2, Yoshio Imai2, Toru Machida2, Ichiro Shirouzu2, Seiki Konishi1, Yasushi Miyashita1 & Naoki Masuda3

The resting-state human brain networks underlie fundamental cognitive functions and consist of complex interactions among brain regions. However, the level of complexity of the resting-state networks has not been quantied, which has prevented comprehensive descriptions of the brain activity as an integrative system. Here, we address this issue by demonstrating that a pairwise maximum entropy model, which takes into account region-specic activity rates and pairwise interactions, can be robustly and accurately tted to resting-state human brain activities obtained by functional magnetic resonance imaging. Furthermore, to validate the approximation of the resting-state networks by the pairwise maximum entropy model, we show that the functional interactions estimated by the pairwise maximum entropy model reect anatomical connexions more accurately than the conventional functional connectivity method. These ndings indicate that a relatively simple statistical model not only captures the structure of the resting-state networks but also provides a possible method to derive physiological information about various large-scale brain networks.

DOI: 10.1038/ncomms2388

A pairwise maximum entropy model accurately describes resting-state human brain networks

1 Department of Physiology, The University of Tokyo School of Medicine, Tokyo 113-0033, Japan. 2 Department of Radiology, NTT Medical Center Tokyo, Tokyo 141-8625, Japan. 3 Department of Mathematical Informatics, The University of Tokyo, Tokyo 113-8656, Japan. Correspondence and requests for materials should be addressed to N.M. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

& 2013 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms2388

During rest, the human brain is not idle but shows a large amount of spontaneously uctuating activity that is highly correlated between multiple brain regions1,2.

Previous studies used positron emission tomography or functional magnetic resonance imaging (fMRI) to reveal that brain regions interact with each other during rest and constitute several resting-state networks (RSNs)36. The RSNs, including the default mode network (DMN)5,6 and the fronto-parietal network (FPN)7,8, are highly reproducible across different healthy individuals9 and are considered to underlie fundamental and intrinsic functions such as self-referential cognitive processes10,11, maintenance of memory12 and attentional cognitive processes8. These functions are considered to originate from complex functional interactions among the brain regions belonging to the RSNs13. However, the level of complexity that is observed in the activities of an ensemble of brain regions and in the structure of these functional interactions has not been quantied for the RSNs. As a result, it remains challenging to comprehensively understand this ensemble brain activity as an integrative large-scale neural system14,15.

For microscopic neural tissues, maximum entropy models (MEMs) have been successful in estimating the complexity of neural activity patterns. For example, the pairwise MEM method (that is, a second-order model) seeks to t a relatively simple binary-state model containing up to second-order interaction terms (that is, ring rates and pairwise interaction) to empirical data of spatiotemporal spike trains. If the pairwise MEM is accurately tted to the data, activity patterns of neurons can be described with a combination of averaged activity at each recording site and pairwise correlation. If it does not, we cannot ignore higher-order interactions such as triplet interactions originating from an external common input, and accurate descriptions of the activity patterns of neurons would require more complicated (that is, higher-order) MEM models. In fact, the pairwise MEM accurately describes ring patterns in the retinal tissues of primates recorded electrophysiologically in vitro1621, ring patterns and local eld potentials (LFPs) in human cortical tissues in vitro22 and large-scale ring patterns in the visual cortex of monkeys and cats in vivo23,24. These ndings suggest the possibility that the human brain activity patterns during rest are not so complex and can be accurately described by pairwise MEMs.

If the pairwise MEM accurately explains large-scale brain activity patterns during rest, the MEM method will give us much richer information about functional interactions in the RSNs than the widely used functional connectivity (FC) analysis, which is based on Pearsons correlation coefcient between a pair of brain regions. Although FC often implies broader classes of methods for inferring connectivity in the brain, we use this term to refer to the methods on the basis of Pearsons correlation coefcient in this paper. The FC-based analysis (based on Pearsons correlation coefcient) has successfully revealed that specic functional connexions are enhanced during specic cognitive processes25,26, and the FC map serves as a promising indicator of several diseases27,28, stages of development29,30 and levels of intelligence31. However, because FC is calculated as Pearsons correlation coefcient between activities of pairs of regions3,4, the FC-based method is founded on the implicit assumption that the pairwise interactions are independent of each other32. If different pairwise interactions inuence each other, the observed correlation between regions A and B may be a natural consequence of, for example, the correlation between regions A and C and that between regions B and C. In fact, a study employing monkeys shows that, even when a pair of brain regions does not have a direct anatomical connexion, the FC between them is often large if the regions receive common input from a

third region33. Other studies also show that a signicant FC value between a pair of regions does not distinguish a direct (that is, monosynaptic) connexion from an indirect (that is, polysynaptic) connexion34,35. The FC-based methods may discard possibly rich information that would be revealed if we take into account that functional interactions inuence each other8.

Unlike the functional interactions estimated by the FC-based method, functional interactions estimated by the pairwise MEM method are not a simple collection of pairwise interactions that are determined independently of each other. The method infers organization of functional interactions based on global activity patterns (that is, activities of more than two sites that are considered simultaneously), not on the assumption that activity patterns of region pairs are independent of those of other pairs. Therefore, if the pairwise MEM can accurately describe the RSNs, the MEM method is expected to provide a better method for assessing integrative network structure than the FC-based method.

In the present study, we rst quantify the complexity inherent in activities of an ensemble of brain regions during rest by tting the MEM to resting-state fMRI data in the RSNs. Then, to validate the approximation of the RSNs by the pairwise MEM, we demonstrate that the interactions estimated by the pairwise MEM are more physiologically informative than those estimated by FC in the sense that anatomical connexions are more accurately predicted. The pairwise MEM even outperforms other methods that take into account global activity patterns beyond a collection of independently estimated pairwise interactions, including the partial correlation method36,37 and the mutual information (MI) method3739.

ResultsAcquisition and processing of data. We applied the pairwise MEM to fMRI signals obtained from the brain regions belonging to the DMN5,6 and FPN7,8, two representative RSNs. We separately applied the pairwise MEM to the DMN and FPN, which are relatively independent of each other in the sense that they are thought to have different functions15,40. The coordinates of the brain regions were anatomically dened based on a meta-analysis of resting-state fMRI data obtained from 183 healthy participants41,42. The DMN and FPN consisted of 12 and 11 regions, respectively (Table 1). We obtained B45 h of resting-state fMRI data from six healthy young adult participants (B7.5 h per participant). Using the data, we extracted the low-frequency signals (0.010.1 Hz) from all the brain regions included in the DMN or FPN (Fig. 1a), binarized the signals (Fig. 1b) and tted the MEM to the binarized signals (Fig. 1c).

Accurate tting of the pairwise MEM to resting-state fMRI signals. In Fig. 2a,b, each data point represents the estimated and empirical probabilities that a binarized activity pattern (that is, state vector in Fig. 1c) occurs. If the points are on the diagonal (lines in Fig. 2a,b), the model would perfectly explain the empirical activity patterns. For comparison, the results for the MEM estimated under the restriction that different regions do not interact (that is, independent MEM) are also presented in the gures. Figure 2a,b suggests that, in both the DMN and FPN, the pairwise MEM explains the empirical frequency of the activity patterns much more accurately than the independent MEM. In both the DMN and FPN, the probability of the pattern estimated by the pairwise MEM was generally closer to the empirical probability than that estimated by the independent MEM for both rare and frequent patterns (Fig. 2c,d). The pairwise MEM results in a high accuracy (85% for the DMN and 94% for the FPN). An accuracy of 85%, for example, indicates that the application of the

2 NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2013 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms2388 ARTICLE

Table 1 | Coordinates of brain regions in the RSNs.

Network Region MNI Coordinates

X Y Z DMN Anterior medial PFC 1 55 26

Ventro-medial PFC 3 40 0

Left SFG 14 36 59

Right SFG 17 35 58 Left ITG 62 33 20

Right ITG 66 17 19

Left parahippocampal gyrus 22 26 21

Right parahippocampal gyrus 25 26 18

PCC 2 29 39

Left lateral parietal 47 71 35

Right lateral parietal 54 61 36

Posterior cingulate 3 53 6

FPN Left DLPFC 48 21 38

Right DLPFC 43 21 38 Left MFG 41 1 39

Right MFG 41 1 39 Middle CC 0 31 31

Left IPL 52 54 36

Right IPL 52 51 43

Left IPS 31 63 42

Right IPS 30 65 39

Left precuneus 9 76 36

Right precuneus 10 73 39

CC, cingulate cortex; DLPFC, dorso-lateral prefrontal cortex; DMN, default mode network; FPN, fronto-parietal network; IPL, inferior parietal lobe; IPS, inferior parietal sulcus; ITG, inferior temporal gyrus; MFG, middle frontal gyrus; MNI, Montreal Neurologic Institute; PCC, posterior cingulate cortex; PFC, prefrontal cortex; RSN, resting-state network; SFG, superior frontal gyrus.

Region 1

Region 2

Region N

Merge data obtained from all the subjects and t the pairwise MEM

Probabilityof network state

P(V i)

Model

Empirical

Resting-state BOLD signal

0 0

0

1 0

0

1 1

0

1 1

1

Time

Binarization

Network state Vi

Region 2

Region N

Region 1 Region 2

Region N

Region 1

Region 1 Region 2

Region N

1 1

0

0

1

1 0

1

0 0

1

1

1 1

1

0

1 2 3

Time

T

Interaction matrix Jij

Figure 1 | Fitting the pairwise MEM to fMRI data. (a,b) Resting-state BOLD signals are extracted from anatomically dened regions (a) and binarized (b). (c) If the pairwise MEM is accurately tted to the binarized data, the estimated probability distribution of the network state will resemble the empirical distribution (top). The pairwise MEM also provides an interaction matrix (bottom).

pairwise MEM reduces the distance between the estimated and empirical distributions of patterns by 85% as compared with the independent MEM (see equation (2) in the Methods). These results show that the pairwise MEM accurately describes the fMRI signals in the RSNs.

Robustness of accurate tting of the pairwise MEM. We then examined the robustness of the accurate tting of the pairwise MEM from three perspectives. First, we evaluated the accuracy of t for various values of the threshold for binarization, because the pairwise MEM requires binarization of the originally continuous fMRI signals (Fig. 1b). The distributions of the original continuous signals for the two RSNs are shown in Figure 3a. We then set various threshold values for binarization within the range of the fMRI signal, and tested the dependence of the so-called estimation reliability on the threshold. The estimation reliability represents the precision with which the parameters in the pairwise MEM are estimated (see equation (3) in the Methods). The pairwise MEM estimated the parameter values with high reliability (499%) when the threshold lay approximately between

0.15 and 0.15 (Fig. 3b). Consistent with this result, the tting of the pairwise MEM to the data was accurate for threshold values in this range (480% for the DMN and 490% for the FPN; Fig. 3c).

The threshold of 0.15, for example, implies that the accuracy of

t remains high even if the probabilities of the two binarized states of a brain region are quite uneven (that is, activated with probability 0.78 and not activated with probability 0.22). These results suggest that the accuracy of t is robust as long as the threshold is located between 0.15 to 0.15. In the following, we

set the threshold to 0.1 because this value maximized the accuracy of t of the pairwise MEM (Fig. 3c).

Second, we conrmed that the high accuracy of our results was not due to a fallacy of the pairwise MEM when it is applied to too

small populations. When the accuracy of the tting linearly decreases with the number of regions N, one theory suggests that the accurate tting may be an artifact of a small N43. In this case, the validity of the pairwise MEM cannot be directly extrapolated to populations with a large N43. For our data, the accuracy decreases sub-linearly with N at least for the DMN (Fig. 3d). Furthermore, the theory indicates that the pairwise MEM can be informative about empirical data only when it is accurately tted to data with N signicantly larger than Nc 1/ vdt, where v and dt

are the averaged activation rate and the length of the time bin, respectively43. In the present analysis, v 0.045 s 1 and

v 0.041 s 1 for the DMN and FPN, respectively, and

dt 9.045 s. Therefore, NcE2.5 and 2.7 for the DMN and FPN,

respectively. Figure 3d suggests that the pairwise MEM was accurate for the two networks up to N 12 and 11, which are

much larger than Nc. Therefore, the high accuracy of the pairwise MEM for our data is not an artifact of a small N value.

Third, we examined the dependence of the results on the participants. We split the fMRI data into two groups, that is, the data obtained from participants no. 13 and those obtained from participants no. 46. The interaction weights (that is, elements of the estimated interaction matrix, denoted by Jij in equation (1);

also see Fig. 1c) obtained from the pairwise MEM were similar between the two groups (Pearsons correlation coefcient r 0.96, Po10 3 in Fig. 3e; r 0.93, Po10 3 in Fig. 3f).

Therefore, our results are robust with respect to the participants.

In summary, the pairwise MEM accurately and robustly explains the collective activities of the brain regions in the RSNs.

Interaction matrices derived by the pairwise MEM and alternative methods. If the functional interactions estimated by the pairwise MEM are closely relevant to some physiological

NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

& 2013 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms2388

DMN

DMN

FPN

FPN

1

Estimation reliability

1 Pairwise MEM

Pairwise MEM

0.04

DMN FPN

0.9

102

0.03

0.8

P model

104

P model

1

102 104 106 108

Probability

0.02

0.7

0.01

106 108

0.6

DMN

Independent MEM

Independent MEM

FPN

0

0.5

0.8 0.4 0 0.4 0.8

0.2 0.1 0 0.1

0.2

104

102

1

104

Pempirical

102

1

fMRI Signal

Threshold for binarization

Pempirical

1

1

1

102 104 106 108

1

102 104 106 108

Accuracy of fit

0.9

0.9

P(V i)

Accuracy of fit

0.8

P(V i)

0.7

0.6

DMN

0.8

FPN

0.5

0.2 0.1 0 0.1 0.2

3

6 9 12

Vi Vi

Threshold for binarization

N

Empirical data

Pairwise MEM Independent MEM

J ij for participants no. 46

J ij for participants no. 46

0.4

DMN

0.4

FPN

DMN FPN

Figure 2 | Accurate tting of the pairwise MEM. (a,b) Comparison between the probabilities of network states for empirical data and those estimated for the pairwise and independent MEMs. The red and grey dots indicate the results for the pairwise and independent MEMs, respectively. (c,d) The probability of a network state for the MEMs, P(Vi), is plotted against the network state, (Vi). The network states are shown in descending order of empirical probability. The red and grey lines represent the results for the pairwise and independent MEMs, respectively. The green lines represent the empirical data.

0.2

0.2

0

0

0.2 0.2 0

0.2

0.2

0.4

0.2 0 0.2 0.4

Jij for participants no. 13 Jij for participants no. 13

Figure 3 | Robustness of tting. In panels (ad) solid and dashed lines represent the results for the DMN and FPN, respectively. (a) Distributions of the preprocessed resting-state fMRI signals. (b) Relationships between the estimation reliability for the pairwise MEM (equation (3)) and the threshold for binarization. (c) Relationships between the accuracy of tting for the pairwise MEM (equation (2)) and the threshold for binarization. (d) Dependence of the accuracy of the tting for the pairwise MEM on the number of regions. Symbols:average, Error bar: s.d. For example, there are 220 different sub-networks (combinations) of N 3 regions for the DMN.

(e,f) Scattergrams of the interaction weights obtained from the pairwise MEM tted to the data for participants no. 13 and those obtained from participants no. 46 show that the results do not vary by subject.

basis, the relationship would provide further validation of the approximation of the RSNs by the pairwise MEM. Therefore, we next examined the possible physiological bases of the interaction weights between pairs of regions inferred by the pairwise MEM. For comparison, we referred to large-scale anatomical connectivity determined from previous diffusion tensor imaging (DTI) of the brains of 80 healthy young human adults44. To assess the performance of the pairwise MEM, we also generated four types of interaction matrices as controls, that is, those based on the FC method, those based on the inverse Gaussian model, those based on the partial correlation method36,37 and those based on the MI method37,39 (see Methods). It should be noted that the four methods do not require binarization of the fMRI signals.

The anatomical connectivity maps and the interaction matrices estimated by the ve methods are shown in Fig. 4. For both the DMN (Fig. 4a) and FPN (Fig. 4b), the |Jij| values estimated by the pairwise MEM are mainly large for pairs of brain regions with direct anatomical connexion, whereas, in the other methods, |Jij|

is large both for pairs with direct anatomical connexion and those without such connexion. Therefore, the pairwise MEM seems to outperform the other four methods in estimating the anatomical connexions if we regard both positive and negative large values of Jij as predicting the presence of anatomical connexion. The use of both positive and negative values is consistent with the anatomical evidence, because anatomical results derived from DTI are thought to reect either excitatory or inhibitory connectivity between pairs of regions.

Prediction of anatomical connexions by the pairwise MEM and the other methods. To analyse the results more quantitatively, we

compared the accuracy of the ve methods in estimating the presence or absence of the anatomical connexions between pairs of regions. Figure 5a,b shows histograms of the values of interaction weights (that is, values of the entries of the matrices shown in Fig. 4) estimated by the ve methods. As the absolute value of the functional interaction weight seems to predict the anatomical connexions better than the raw values of them, at least for all the methods except FC, we show the histograms of the absolute values of the interaction weights. However, for the FC, we also show the histograms of the raw values of the interaction weights (that is, correlation coefcient) because they are conventionally used34,41,45.

Under the pairwise MEM, anatomically connected pairs of regions in both the DMN and FPN tended to have larger absolute values of interaction weights than pairs that were not directly connected in anatomy. The interaction weights for the anatomically connected pairs were signicantly larger than those for the anatomically unconnected pairs on average

4 NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2013 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms2388 ARTICLE

DMN

Pairwise MEM FC

Inverse gaussian model Mutual information

Partial correlation

Anatomical connexion

ant. mPFC vmPFC Lt. SFG Rt. SFG

Lt. ITG Rt. ITG Lt. parahippo Rt. parahippo

PCC Lt lateral parietal Rt. lateral parietal

Rost. cingulate

ant. mPFC

vmPFC

Lt. SFG

Rt. SFG

Lt. ITG

Rt. ITG

Lt. parahippo

Rt. parahippo

Jij

Zij

PCC

Lt lateral parietal

Rt. lateral parietal

Rost. cingulate

3 3

0.5

0.5

rij

ij1 MIij

50 0.3

0

0.6

50

0.6

FPN

Anatomical connexion

Pairwise MEM FC

Lt. DLPFC Rt. DLPFC

Lt. MFG Rt. MFG

Mid. CC

Lt. IPL Rt. IPL

Lt. IPS Rt. IPS Lt. precuneus Rt. precuneus

Lt. DLPFC

Rt. DLPFC

Lt. MFG

Rt. MFG

Mid. CC

Lt. IPL

Rt. IPL

Lt. IPS

Rt. IPS

Lt. precuneus

Rt. precuneus

Jij

Zij

3 3

0.5

0.5

Inverse gaussian model

50

Partial correlation

Mutual information

rij

ij1 MIij

50

0.6

0.6

0

0.3

Figure 4 | Anatomical connexions and interaction matrices. (a,b) Upper left panels show anatomical connexion maps. The other panels show the interaction matrices based on the pairwise MEM, FC, inverse Gaussian model, partial correlation method and MI method. The abbreviations of the regions are explained in Table 1.

(Po0.001, two-sample t-tests in both the DMN and FPN; Fig. 5c,d).

The partial correlation and MI methods yielded similar signicant differences between the anatomically connected pairs and unconnected pairs (for both models, Po0.001 in DMN,

Po0.05 in FPN, two-sample t-tests; Fig. 5c,d. Also see Fig. 5a,b). In contrast, the FC and inverse Gaussian models did not predict anatomical connexions as accurately as the pairwise MEM

(Fig. 5a,b). In the DMN, the interaction weights estimated by these two methods had larger values when region pairs were anatomically connected than unconnected (FC: Po0.05, absolute value of FC: Po0.01, inverse Gaussian model: Po0.01, two-sample t-tests). However, in the FPN, the difference was insignicant for both of the two methods (P40.3).

For further comparison the ve methods, we carried out a receiver-operating characteristic (ROC) analysis, which is

NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

& 2013 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms2388

DMN

DMN

12

12

0

No. of connexions

3

0 Inverse gaussian model

8

***

* **

8

8

0.6

Z Zij

ij

2

4

4

4

0.4

J ij

1

0.2

0 0 1 2 3

4

0.3 0 0.3

0.6

0 0 0.3

Zij

0.6

Jij

Zij

0

Pairwise MEM

0 0 40 80 120 160

FC

Pairwise MEM

FC

0.6

0.4

0.2

16

20

**

0.4

0.2

10

***

***

30

12

15

8

6

1

20

8

10

ij

r ij

MI ij

4

4

5

2

10

0 0 0.2 0.4 0.6

0.8

0 0 0.2 0.4

ij1

rij

MIij

0 Partial correlation

0 Mutual information

Inverse gaussian model

Partial correlation

Mutual information

FPN

FPN

12

3

0

0.6

Zij

Zij

No. of connexions

8

***

8

8

2

0.4

4

4

4

J ij

1

0.2

0.6

FC

0 0 1 2 3

4

0.3

12

0

0 0.3

Jij Zij

0 0 0.3 0.6

Zij

0 Pairwise

MEM

FC

Pairwise MEM

40

15

**

*

12

16

30

12

8

10

20

ij

r ij

0.4

0.2

MI ij

0.4

0.2

1

8

4

5

10

4

0 0 20 40 60 80

0 0 0.2 0.4 0.6

0.8

0 0 0.2 0.4

0.6

0 Inverse gaussian model

0 Partial correlation

0 Mutual information

ij1

rij

MIij

Inverse gaussian model

Partial correlation

Mutual information

Anatomically connected pairs

Anatomically unconnected pairs

Figure 5 | Statistical comparison of interaction matrices with anatomical connexions. In all the panels, the black and grey bars indicate the results for the anatomically connected pairs and anatomically unconnected pairs, respectively. (a,b) Histograms of the interaction weights for anatomically connected and unconnected pairs. (c,d) Comparison of the average interaction weight between anatomically connected and unconnected pairs. Error bars: s.d., ***Po0.005, **Po0.01, *Po0.05, in two-sample t-tests.

sensitive to the difference in the distribution of the functional interactions between the anatomically connected and unconnected pairs (Fig. 6a,b). In short, the curves show the relationship between the false positive (that is, anatomically unconnected pairs whose interaction weights exceed a threshold) and true positive (that is, anatomically connected pairs whose interaction weights exceed the same threshold; also called the sensitivity) when we vary the threshold for classifying the region pairs into anatomically connected and unconnected groups. When the false positive is small and the true positive is large for a threshold value, the classication is accurate such that the interaction weights determined by an estimation method are relatively consistent with the anatomy. In both the DMN and FPN, the area under the curve (AUC) for the pairwise MEM was signicantly larger than those for the other four methods (DMN: Z47.3, Po1.0 10 11, Bonferroni-corrected; FPN: Z49.4,

Po1.0 10 16, Bonferroni-corrected; Fig. 6c). There was no

signicant difference between any other pair of the methods.

Validity of binarizing fMRI signals. We binarized the fMRI signals to implement the pairwise MEM. However, fMRI signals originally represent continuous blood ows to brain regions, and the validity of the binarization is unclear. Therefore, we calculated the accuracy of t of pairwise MEMs in which the signals were classied to three discrete levels, which we call the trinary MEM, by using a standard procedure described in Methods46. When the signals recorded at the different brain regions and time points were evenly divided into three groups (Supplementary Fig. S1), the accuracy of t was the highest. However, the accuracy of t in this case was approximately equal to 0.55 (Supplementary Fig. S2) and signicantly lower than those obtained for the binary pairwise MEM for both the DMN and FPN (Po10 3 in two-sample t-tests).

Despite the lower accuracy of the trinary MEM, the trinary MEM might be better at accurately predicting anatomical connectivity than the binary MEM. Therefore, we examined the relationship between the functional interactions estimated by the

6 NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2013 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms2388 ARTICLE

DMN

DMN

FPN

FPN

1

1

Pairwise MEM FC

0.8

0.8

True positive

Inverse gaussian model

Mutual information

Partial correlation

0.6

0.6

0.4

0.4

0.2

0.2

0 0

0.2

0.4

0.6

0.8

1

0 0 0.2 0.4 0.6 0.8

Absolute value of FC

Inverse gaussian model

1

False positive

False positive

1

***

***

0.9

Pairwise MEM

FC

Absolute value of FC

0.8

AUC

0.7

Parital correlation

Mutual information

0.6

0.5

Jij Zij rij MIij

Zij

ij1 Jij Zij rij MIij

Zij ij1

Jij Zij rij

MIij

Zij ij1

Figure 6 | ROC analysis. The colour code is as follows: red, pairwise MEM; green, raw values of FC; light green, absolute values of FC; blue, inverse Gaussian model; yellow, partial correlation method; purple, MI method. (a,b) ROC curves for the classication into anatomically connected and unconnected region pairs based on the interaction weights. (c) AUCs of the ROC curves shown in panels (a,b). Error bars: s.d., ***Po0.005, in two-sample t-tests with the Bonferroni correction.

trinary MEM and the anatomical connectivity. For the trinary MEM, the functional interactions between anatomically connected pairs were larger than those between anatomically unconnected pairs (Po0.05 in two-sample t-tests, Supplementary

Fig. S3). However, ROC analysis revealed that the AUC for the trinary MEM was signicantly smaller than that for the binary MEM (Po10 3 in two-sample t-tests, Supplementary Fig. S4).

These results suggest that the binary pairwise MEM provides more information about anatomical connectivity than the trinary pairwise MEM. Such a difference usually arises because the trinary MEM requires a larger amount of data than the binary MEM (that is, 3N versus 2N data points). However, even the trinary MEM with N 4 and N 5, for which there are relatively

few activity patterns (that is, 34 81 and 35 243 patterns),

yields the accuracy of t values o0.6. These accuracy of t values are much smaller than those for the binary MEM with comparable numbers of activity patterns: the accuracy of t was 40.9 for the binary MEM with N 8, with which there are

28 256 activity patterns (Fig. 3d). Therefore, we conclude that

the binary MEM better captures fMRI signals than the trinary MEM and presumably MEMs that allow more discrete levels.

DiscussionWe have demonstrated that the pairwise MEM accurately describes activity in the human RSNs. The model was tted to both the DMN and the FPN with high accuracy and robustness. Furthermore, functional interaction matrices derived from the pairwise MEM were similar to the anatomical connexion maps. The agreement between the estimated matrices of functional interactions and the anatomical maps is better with the pairwise MEM than with the four other methods: the FC based on Pearsons correlation coefcients, inverse Gaussian model, partial correlation method and MI method. These ndings suggest that the large-scale human brain networks during rest can be captured by a relatively simple second-order model.

Our results extend the previous ndings that the pairwise MEM accurately describes activity patterns observed in a slice or small parts of the brain1619,2224. In particular, Tang and colleagues22 showed that the pairwise MEM approximates the statistics of multiunit LFPs recorded from cultured slices of the human brain cortex. As LFPs are considered to be predictive of blood-oxygenation-level-dependent signals (that is, fMRI signals)47, it is reasonable that the present study has successfully tted the pairwise MEM to fMRI signals. Although ne-scale neural networks (o300 mm) seem to require more complex models, such as higher-order MEMs19,24, in terms of the spatial scale, we have extended the applicability of the pairwise MEM method from microscopic neural populations (B500 mm)

to macroscopic populations at the level of the entire brain (B10 cm).

In previous studies, large FC values have been suggested to imply the presence of direct anatomical connexion between regions2,34,45. In particular, a study using high-resolution fMRI conducted an ROC analysis similar to ours, and reported that the FC classies the anatomical connexions with high accuracy (AUC 0.79)34. Our present study produced similar accuracy

values for the FC (AUCs in the DMN and FPN of 0.73 and 0.67, respectively; Fig. 6c) but also showed that the pairwise MEM gives higher accuracy (AUC in the DMN and FPN of 0.89 and 0.85, respectively; Fig. 6c). These results suggest that the pairwise MEM captures the anatomical properties of large-scale brain networks to a better extent than the FC-based method does, although the pairwise MEM is a statistical model and not based on the physiology of the brain.

As second-order interaction implies networks consisting of nodes and links, the pairwise MEM may serve as a tool for investigating large-scale human brain networks. As previous studies did for networks generated by the resting-state FC (see a review14), we can in principle calculate various graph-theoretical quantities for the RSNs dened by the pairwise MEM. This approach calls for a sufcient amount of resting-state fMRI data

NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

& 2013 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms2388

to estimate an interaction matrix among a relatively large number of brain regions and efcient algorithms to estimate the MEM for large N, such as those developed in previous studies19,48,49. We can also apply the pairwise MEM to fMRI data measured during psychological tasks, though several modications of the method may be in order. In previous fMRI studies, such task-specic functional interactions were estimated by psychophysiological interaction analysis50, Granger causality analysis51 and dynamic causal modelling52. Among these methods, psychophysiological interaction analysis, like FC-based methods, is a seed-based analysis and does not allow us to collectively estimate functional interactions among the entire set of recorded sites. The use of dynamic causal modelling and Granger causality modelling still seems to be limited8,53. Therefore, the pairwise MEM may serve as an alternative method for inferring task-specic brain networks. In particular, recent theoretical developments of the MEM incorporate causality (signal ow from one region to another)54,55, which will be an indispensable component for modelling task-related data.

Methods

Participants and fMRI data acquisition. Six healthy right-handed subjects (aged 2023 years; three males) participated in the experiments after written informed consent was obtained. The MRI scanning was conducted using a 3T MRI scanner (Philips Achieva X 3T Rel. 2.6, Best, The Netherland). T1-weighted structural images were obtained for anatomical reference (resolution 0.81 0.81 1.20

mm3). Functional imaging used gradient-echo echo-planar sequences(TR 9.045 s, TE 35 ms, ip angle 901, resolution 2 2 2 mm3, 75 slices).

The entire procedure for the MRI scanning was approved by the institutional review board of The University of Tokyo School of Medicine.

During the functional imaging, the participants were instructed to passively view a xation point on the screen. One passive xation session took B5 min, and each participant underwent 90 sessions (B7.5 h) on 68 separate days. In each session, we discarded the rst ve images (9.045 s/volume 5 volumes 45.225 s) to

exclude the inuence of transient processes before the equilibrium of longitudinal magnetization. In total, 17,820 volumes (2,970 volumes per participant) of resting-state functional images were obtained.

Preprocessing of fMRI data. In line with previous studies of resting-state fMRI25,26,30, we preprocessed the fMRI data as follows. First, for each session and each participant, the functional images were realigned, their slice-timing was corrected and they were normalized to the standard template image (ICBM 152) by SPM8 (http://www.fil.ion.ucl.ac.uk/spm/

Web End =www.l.ion.ucl.ac.uk/spm/). Second, for each session, after temporal bandpass ltering (0.010.1 Hz) with Matlab scripts written in-house, the data underwent spatial smoothing (FWHM 8 mm) in SPM8. Third, the smoothed

data originating from all the sessions were combined for each participant. Finally, for each participant, the combined data were corrected for their head motion, whole-brain signals, ventricular signals, white matter signals and run effect. Therefore, the fMRI signals used in our analysis have the same unit as that of the so-called beta-value5658 except for a normalization factor.

We determined the coordinates of the brain regions belonging to the DMN and FPN based on a previous meta-analysis of resting-state fMRI studies30,42 (the DMN: 12 regions, the FPN: 11 regions, Table 1). We subjected the preprocessed data corresponding to these brain regions to the following analysis.

Pairwise MEM. We t the pairwise MEM to the resting-state fMRI data in essentially the same manner as did the previous studies of microscopic neural activities16,17,2224. We added a constant to the preprocessed signal at each region such that the average of the preprocessed signals over T 17,820 snapshots (that is,

functional images) is equal to zero at each region. Then the signals were binarized, because the pairwise MEM deals with snapshots of binarized signals recorded simultaneously at multiple sites.

We set the bin width to 9.045 s, that is, one whole-brain image was taken in one repetition time (TR) of 9.045 s, because a long bin width is effective at decorrelating data that were recorded in adjacent time points. With a small bin width, the fMRI signals (that is, network states) in close time points would be highly correlated with each other. Therefore, the use of such a small bin width would not increase the effective number of samples.

We experimented with various threshold values for binarization, and chose a value of 0.1, because it maximized the accuracy of t (Fig. 3b,c). The binarized activity of brain region i at discrete time t, denoted by sit, is either on ( 1) or off

(0). The network state (that is, pattern) at time t is given in vector form byVt st1; st2; . . . ; stN ; where N is the number of the brain regions of interest (that

is, N 12 in the DMN and N 11 in the FPN). It should be noted that there are 2N

possible network states. The empirical activation rate of region i, denoted by /siS,

is given by hsii 1/T PTt 1 sti. The empirical pairwise joint activation rate of

regions i and j, denoted by /sisjS, is given by hsisji 1/T PTt 1 stistj.

The pairwise MEM maximizes the entropy of the distribution of activity patterns under the restriction that /siS and /sisjS (1rirjrN) are preserved. It is known that such a distribution has the form PVi e EV / P2i 1 e EV , where P(Vi) is

the probability of network state Vi, andEVi X

N

i 1

hisi 1/2 X

N

i 1

X

j 1; j6i

N

Jijs s 1

hi represents the tendency of activation at region i, and Jij, called the interaction weight, represents functional interaction between regions i and j. Positive and negative values of Jij are interpreted as excitatory and inhibitory interactions, respectively. For an estimated probability distribution, the expected activation rate, /siSm, and the expected pairwise joint activation rate, /si sjSm, are given by

hsiim P2j 1 siVjPVj and hsisjim P2k 1 siVksjVkPVk, respectively,

where si(Vj) indicates the activity (either 1 or 0) of si under network state Vj.

We calculated hi and Jij by iteratively adjusting /siSm and /sisjSm toward /siS and /sisjS, respectively, by using a gradient ascent algorithm. The iteration scheme is given by hnewi holdi a loghsii/hsiim and Jnewij Joldij a loghsisji/hsisjim, where a is the learning rate. We set a 0.75 to reduce computation time. We

conrmed that the results were robust with respect to the choice of a ranging from0.1 to 1.0.

As in previous studies16,17,2224, we evaluated the effectiveness of the pairwise

MEM by using two information theoretic quantities. To dene them, we also obtained the MEM without pairwise interaction between regions (that is, independent MEM) by determining hi (1rirN) under the condition that Jij 0

(1ri, jrN). For the independent MEM, /sisjSm is equal to /siSm/sjSm.

Using the resultant goodness of t obtained for the two types of the MEMs, we dene the accuracy index for the pairwise MEM by

rD D1 D2/D1 2

where Dk P2i 1 PNVi log2PNVi/PkVi is the KullbackLeibler divergence

between the probability distribution of the network state in the k th order model (k 1 and 2 for the independent and pairwise MEMs, respectively) and the empirical

distribution of the network state, denoted by PN. To evaluate reliability of the tting, we calculate rS (S1 S2)/(S1 SN), where Sk P2i 1 PkVi logPkVi is the

entropy of the distribution of the network state in the k th order model. Note that SN is the entropy of the empirical data. Using rS, we dene the estimation reliability by

rS/rD 3

The estimation reliability is equal to 1 if the values of hi and Jij are estimated without error59.

Cross validation. As a large sample size is desirable for statistical tting, we trained the pairwise MEM using the entire data and tested it against the same data. To exclude the possibility of over tting, we carried out cross validation as follows: (i) We divided the T observed snapshots (that is, network states) into the two subsets of equal size, that is, T/2, such that each snapshot independently belongs to either subset with probability 1/2. It should be noted that we are not concerned with the temporal structure of the data. (ii) We calculated the accuracy of t by training the MEM using one subset and testing the estimated model against the other subset.(iii) We calculated the accuracy of t by using one of the two subsets for both training and testing. (iv) We repeated procedures (i), (ii) and (iii) ten times and compared the accuracy of t between the cross validation case (that is, (ii)) and the dependent training-testing case (that is, (iii)). We did not nd signicant differences between the two cases (cross validation versus dependent training-test cases:0.670.0068 versus 0.680.0058 for the DMN, P40.41; 0.760.0071 versus0.760.0051 for the FPN, P40.35, means.d. in paired t-tests). Therefore, we concluded that the use of the same data in both training and testing did not give rise to serious problems in the present study.

Trinary pairwise MEM. To examine the validity of binarizing fMRI signals, we estimated the accuracy of t of the trinary pairwise MEM as follows: First, we discretized the fMRI signal at each brain region i and each time t into one of the three different levels by thresholding. The discretized activity sit is set to 1, 0 and

1 if the fMRI signal was larger than e, between e and e and smaller than e,

respectively. Then, we tted the trinary pairwise MEM in the almost same manner as that for tting the binary pairwise MEM46.

Conventional FC. As a control of the functional interaction, we calculated the FC based on the Pearsons correlation coefcient between the activity of every pair of regions included in the DMN or FPN. The calculation procedures were the same as those used in previous studies of the resting-state FC25,26. After calculating the FC for each scanning session and each participant, we carried out Fishers transformation28 and averaged the Z-value-based FC across sessions and participants. As is consistent with a previous study60, the FC matrices for the DMN

8 NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2013 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms2388 ARTICLE

and FPN were similar across participants as well as sessions for our data: Pearsons correlation coefcient between an arbitrary pair of the FC matrices of possibly different participants and sections was 40.81. The large correlation values suggest that the averaged FC is a good representation of the FC for a single participant and session.

Inverse Gaussian model. As another control, we estimated functional interaction using the inverse Gaussian model. In this model, the probability density function of the brain activity is assumed to obey the multivariate Gaussian distribution given by PIGMXt / e 1/2X XG

X

X , where Xt [X1(t), X2(t),yXN(t)] is the

N-dimensional row vector representing the continuous-valued activity at the N regions at time t,

X is the row vector of the activity at the N regions averaged over sessions and participants, and G 1/T PTt 1Xt

XXt

XT is the N N

covariance matrix. In this model, we dene the interaction weight between regions i and j as the (i, j) element of matrix G 1.

Partial correlation method. As another control, we estimated functional interactions using the partial correlation method36,37. Using the inverse of the covariance matrix, we dened the partial correlation between brain regions i and j, denoted by rij, by rij G 1ij/

G 1iiG 1jj

q .

MI method. As another control, we estimated functional interactions using the MI method. In accordance with the previous studies3739, we rst estimated the MI between regions i and j at frequency ok by MIijok 1/2 ln1 mCohijok,

where mCohij (ok) represents the multiple coherence between fMRI signals at regions i and j at frequency ok37. We then obtained the MI between i and j, MIij, by averaging MIij(ok) over ok (0.01 Hzrokr0.1 Hz).

Comparison with anatomical connexions. We used the anatomical connexions among the regions in the DMN and FPN that are determined by a previous DTI study employing 80 healthy young adults42. In the study44, the authors obtained the DTI data in a Siemens Sonata 1.5T MRI scanner (Siemens Medical Systems, Germany) by using a twice-refocused single-shot Echo-Planar Imaging-based sequence (3 mm slice thickness with no inter-slice gap, 40 axial slices,TR 6,400 ms, TE 88 ms, 6 diffusion directions). The DTI data were

preprocessed in SPM5 and mapped to a template of ICBM 152 in the Montreal Neurological Institution space. The anatomical connexions were determined with the use of DTI deterministic tractography.

We classied the estimated values of the interaction weights between region pairs into a group of anatomically connected pairs and a group of anatomically unconnected pairs. We compared the two groups by using the two-sample t-test. As the anatomical connexions observed in DTI do not distinguish excitatory and inhibitory interactions, we submitted absolute values of the estimated interaction weights to the t-test. However, for the FC, we also submitted the raw FC values to the t-test, as was done in previous reports34.

By conducting the ROC analysis, we obtained the accuracy with which the magnitude of the interaction weight predicts the presence and absence of the anatomical connexions. We drew ROC curves by plotting the true positive (that is, anatomically connected pairs whose interaction weights exceed a threshold) against the false positive (that is, anatomically unconnected pairs whose interaction weights exceed the same threshold) for various values of the classication threshold. For an ROC curve, the AUC is equal to the area below the ROC curve and represents the classication accuracy of the estimation method34. The deviation of the AUC was estimated based on a bootstrap method.

References

1. Fox, M. D. & Raichle, M. E. Spontaneous uctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neurosci. 8, 700711 (2007).

2. Sporns, O. The human connectome: a complex network. Ann. N Y Acad. Sci. 1224, 109125 (2011).

3. Friston, K. J., Frith, C. D., Liddle, P. F. & Frackowiak, R. S. Functional connectivity: the principal-component analysis of large (PET) data sets.J. Cereb. Blood Flow Metab. 13, 514 (1993).4. Biswal, B., Yetkin, F. Z., Haughton, V. M. & Hyde, J. S. Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537541 (1995).

5. Raichle, M. E. et al. A default mode of brain function. Proc. Natl Acad. Sci. USA 98, 676682 (2001).

6. Greicius, M. D., Krasnow, B., Reiss, A. L. & Menon, V. Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc. Natl Acad. Sci. USA 100, 253258 (2003).

7. Beckmann, C. F., DeLuca, M., Devlin, J. T. & Smith, S. M. Investigations into resting-state connectivity using independent component analysis. Philos. Trans. R Soc. Lond. B Biol. Sci. 360, 10011013 (2005).

8. Cole, D. M., Smith, S. M. & Beckmann, C. F. Advances and pitfalls in the analysis and interpretation of resting-state FMRI data. Front. Syst. Neurosci. 4, 8 (2010).

9. Damoiseaux, J. S. et al. Consistent resting-state networks across healthy subjects. Proc. Natl Acad. Sci. USA 103, 1384813853 (2006).

10. Fair, D. A. et al. The maturing architecture of the brains default network. Proc. Natl Acad. Sci. USA 105, 40284032 (2008).

11. Christoff, K., Gordon, A. M., Smallwood, J., Smith, R. & Schooler, J. W. Experience sampling during fMRI reveals default network and executive system contributions to mind wandering. Proc. Natl Acad. Sci. USA 106, 87198724 (2009).

12. Sestieri, C., Corbetta, M., Romani, G. L. & Shulman, G. L. Episodic memory retrieval, parietal cortex, and the default mode network: functional and topographic analyses. J. Neurosci. 31, 44074420 (2011).

13. Buckner, R. L., Andrews-Hanna, J. R. & Schacter, D. L. The brains default network: anatomy, function, and relevance to disease. Ann. N Y Acad. Sci. 1124, 138 (2008).

14. Bullmore, E. & Sporns, O. Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186198 (2009).

15. Deco, G., Jirsa, V. K. & Mcintosh, A. R. Emerging concepts for the dynamical organization of resting-state activity in the brain. Nat. Rev. Neurosci. 12, 4356 (2011).

16. Shlens, J. et al. The structure of multi-neuron ring patterns in primate retina.J. Neurosci. 26, 82548266 (2006).17. Schneidman, E., Berry, M. J., Segev, R. & Bialek, W. Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440, 10071012 (2006).

18. Ganmor, E., Segev, R. & Schneidman, E. The architecture of functional interaction networks in the retina. J. Neurosci. 31, 30443054 (2011).

19. Ganmor, E., Segev, R. & Schneidman, E. Sparse low-order interaction network underlies a highly correlated and learnable neural population code. Proc. Natl Acad. Sci. USA 108, 96799684 (2011).

20. Schwartz, G., Macke, J., Amodei, D., Tang, H. & Berry, M. J. Low error discrimination using a correlated population code. J. Neurophysiol. 108, 10691088 (2012).

21. Tkacik, G., Prentice, J. S., Balasubramanian, V. & Schneidman, E. Optimal population coding by noisy spiking neurons. Proc. Natl Acad. Sci. USA 107, 1441914424 (2010).

22. Tang, A. et al. A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro. J. Neurosci. 28, 505518 (2008).

23. Yu, S., Huang, D., Singer, W. & Nikolic, D. A small world of neuronal synchrony. Cerebral Cortex 18, 28912901 (2008).

24. Ohiorhenuan, I. E. et al. Sparse coding and high-order correlations in ne-scale cortical networks. Nature 466, 617621 (2010).

25. Tambini, A., Ketz, N. & Davachi, L. Enhanced brain correlations during rest are related to memory for recent experiences. Neuron 65, 280290 (2010).

26. Chadick, J. Z. & Gazzaley, A. Differential coupling of visual cortex with default or frontal-parietal network based on goals. Nat. Neurosci. 14, 830832 (2011).

27. Stam, C. J., Jones, B. F., Nolte, G., Breakspear, M. & Scheltens, P. Small-world networks and functional connectivity in Alzheimers disease. Cereb. Cortex 17, 9299 (2007).

28. Dinstein, I. et al. Disrupted neural synchronization in toddlers with autism. Neuron 70, 12181225 (2011).

29. Supekar, K., Musen, M. & Menon, V. Development of large-scale functional brain networks in children. PLoS Biol. 7, e1000157 (2009).

30. Dosenbach, N. U. F. et al. Prediction of individual brain maturity using fMRI. Science 329, 13581361 (2010).

31. van den Heuvel, M. P., Stam, C. J., Kahn, R. S. & Hulshoff Pol, H. E. Efciency of functional brain networks and intellectual performance. J. Neurosci. 29, 76197624 (2009).

32. Zuo, X. N. et al. Network centrality in the human functional connectome. Cereb. Cortex 22, 18621875 (2012).

33. Adachi, Y. et al. Functional connectivity between anatomically unconnected areas is shaped by collective network-level effects in the macaque cortex. Cereb. Cortex 22, 15861592 (2012).

34. Honey, C. J. et al. Predicting human resting-state functional connectivity from structural connectivity. Proc. Natl Acad. Sci. USA 106, 20352040 (2009).

35. Lu, J. et al. Focal pontine lesions provide evidence that intrinsic functional connectivity reects polysynaptic anatomical pathways. J. Neurosci. 31, 1506515071 (2011).

36. Salvador, R. et al. Neurophysiological architecture of functional magnetic resonance images of human brain. Cereb. Cortex 15, 13321342 (2005).

37. Salvador, R. et al. A simple view of the brain through a frequency-specic functional connectivity measure. NeuroImage 39, 279289 (2008).

38. Salvador, R., Anguera, M., Gomar, J. J., Bullmore, E. T. & Pomarol-Clotet, E. Conditional mutual information maps as descriptors of net connectivity levels in the brain. Front. Neuroinform. 4, 115 (2010).

NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

& 2013 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms2388

39. Lynall, M.-E. et al. Functional connectivity and brain networks in schizophrenia. J. Neurosci. 30, 94779487 (2010).

40. Bullmore, E. et al. Generic aspects of complexity in brain imaging data and other biological systems. NeuroImage 47, 11251134 (2009).

41. Dosenbach, N. U. F. et al. A core system for the implementation of task sets. Neuron 50, 799812 (2006).

42. Fair, D. A. et al. Functional brain networks develop from a local to distributed organization. Plos Comput. Biol. 5, e1000381 (2009).

43. Roudi, Y., Nirenberg, S. & Latham, P. E. Pairwise maximum entropy models for studying large biological systems: when they can work and when they cant. Plos Comput. Biol. 5, e1000380 (2009).

44. Gong, G. et al. Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. Cereb. Cortex 19, 524536 (2009).

45. Vincent, J. L. et al. Intrinsic functional architecture in the anaesthetized monkey brain. Nature 447, 8386 (2007).

46. Mora, T., Walczak, A. M., Bialek, W. & Callan, C. G. Maximum entropy models for antibody diversity. Proc. Natl Acad. Sci. USA 107, 54055410 (2010).

47. Logothetis, N. K., Pauls, J., Augath, M., Trinath, T. & Oeltermann, A. Neurophysiological investigation of the basis of the fMRI signal. Nature 412, 150157 (2001).

48. Tanaka, T. Mean-eld theory of Boltzmann machine learning. Phys. Rev. E 58, 23022310 (1998).

49. Hinton, G. E. Training products of experts by minimizing contrastive divergence. Neural. Comput. 14, 17711800 (2002).

50. Friston, K. J. et al. Psychophysiological and modulatory interactions in neuroimaging. NeuroImage 6, 218229 (1997).

51. Goebel, R., Roebroeck, A., Kim, D.-S. & Formisano, E. Investigating directed cortical interactions in time-resolved fMRI data using vector autoregressive modeling and Granger causality mapping. Magn. Reson. Imaging 21, 12511261 (2003).

52. Friston, K. Dynamic causal modelling. NeuroImage 19, 12731302 (2003).53. Stephan, K. E. et al. Ten simple rules for dynamic causal modeling. NeuroImage 49, 30993109 (2010).

54. Marre, O., Boustani, El, S., Frgnac, Y. & Destexhe, A. Prediction of spatiotemporal patterns of neural activity from pairwise correlations. Phys. Rev. Lett. 102, 138101 (2009).

55. Roudi, Y. & Hertz, J. Mean eld theory for nonequilibrium network reconstruction. Phys. Rev. Lett. 106, 048702 (2011).

56. de Pasquale, F. et al. A cortical core for dynamic integration of functional networks in the resting human brain. Neuron 74, 753764 (2012).

57. Zhou, H. & Desimone, R. Feature-based attention in the frontal eye eld and area V4 during visual search. Neuron 70, 12051217 (2011).

58. Badre, D., Doll, B. B., Long, N. M. & Frank, M. J. Rostrolateral prefrontal cortex and individual differences in uncertainty-driven exploration. Neuron 73, 595607 (2012).

59. Yeh, F. C. et al. Maximum entropy approaches to living neural networks. Entropy 12, 89106 (2010).

60. van Dijk, K. R. A. et al. Intrinsic functional connectivity as a tool for human connectomics: theory, properties, and optimization. J. Neurophysiol. 103, 297321 (2010).

Acknowledgements

We thank Hideaki Shimazaki and Taro Toyoizumi for valuable discussions. We also thank Ms Suzuki for technical assistance of the MRI acquisition. This work was supported by Grants-in-Aid for Scientic Research (23681033, and Innovative Areas Systems Molecular Ethology (No. 20115009)) from MEXT Japan to N.M., a grant from the Japan Society for the Promotion of Science Research Foundation for Young Scientists (222882) to T.W., a Grant-in-Aid for Specially Promoted Research (19002010) to Y.M., a Grant-in-Aid for Scientic Research B (22300134) to S.K. and a research grant from Takeda Science Foundation to Y.M.

Author contributions

T.W. and N.M. designed the research. S.H., H.W., Y.I., T.M., I.S. and S.K. conducted imaging experiments. T.W. and N.M. analysed the data and wrote the manuscript. Y.M. discussed the results and commented on the manuscripts.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Watanabe, T. et al. A pairwise maximum entropy model accurately describes resting-state human brain networks. Nat. Commun. 4:1370 doi: 10.1038/ncomms2388 (2013).

10 NATURE COMMUNICATIONS | 4:1370 | DOI: 10.1038/ncomms2388 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2013 Macmillan Publishers Limited. All rights reserved.