A pairwise maximum entropy model accurately describes restingstate human brain networks.  
Jump to Full Text  
MedLine Citation:

PMID: 23340410 Owner: NLM Status: MEDLINE 
Abstract/OtherAbstract:

The restingstate human brain networks underlie fundamental cognitive functions and consist of complex interactions among brain regions. However, the level of complexity of the restingstate 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 regionspecific activity rates and pairwise interactions, can be robustly and accurately fitted to restingstate human brain activities obtained by functional magnetic resonance imaging. Furthermore, to validate the approximation of the restingstate 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 restingstate networks but also provides a possible method to derive physiological information about various largescale brain networks. 
Authors:

Takamitsu Watanabe; Satoshi Hirose; Hiroyuki Wada; Yoshio Imai; Toru Machida; Ichiro Shirouzu; Seiki Konishi; Yasushi Miyashita; Naoki Masuda 
Publication Detail:

Type: Journal Article; Research Support, NonU.S. Gov't 
Journal Detail:

Title: Nature communications Volume: 4 ISSN: 20411723 ISO Abbreviation: Nat Commun Publication Date: 2013 
Date Detail:

Created Date: 20130123 Completed Date: 20130620 Revised Date: 20130711 
Medline Journal Info:

Nlm Unique ID: 101528555 Medline TA: Nat Commun Country: England 
Other Details:

Languages: eng Pagination: 1370 Citation Subset: IM 
Affiliation:

Department of Physiology, The University of Tokyo School of Medicine, Tokyo 1130033, Japan. 
Export Citation:

APA/MLA Format Download EndNote Download BibTex 
MeSH Terms  
Descriptor/Qualifier:

Brain
/
physiology* Entropy* Female Humans Magnetic Resonance Imaging Male Models, Neurological Nerve Net / physiology* ROC Curve Reproducibility of Results Rest / physiology* Signal Processing, ComputerAssisted Statistics as Topic Young Adult 
Comments/Corrections 
Full Text  
Journal Information Journal ID (nlmta): Nat Commun Journal ID (isoabbrev): Nat Commun ISSN: 20411723 Publisher: Nature Pub. Group 
Article Information Download PDF Copyright © 2013, Nature Publishing Group, a division of Macmillan Publishers Limited. All Rights Reserved. openaccess: Received Day: 15 Month: 05 Year: 2012 Accepted Day: 14 Month: 12 Year: 2012 Electronic publication date: Day: 22 Month: 01 Year: 2013 Volume: 4First Page: 1370 Last Page: PubMed Id: 23340410 ID: 3660654 Publisher Item Identifier: ncomms2388 DOI: 10.1038/ncomms2388 
A pairwise maximum entropy model accurately describes restingstate human brain networks  
Takamitsu Watanabe1  
Satoshi Hirose1  
Hiroyuki Wada2  
Yoshio Imai2  
Toru Machida2  
Ichiro Shirouzu2  
Seiki Konishi1  
Yasushi Miyashita1  
Naoki Masudaa3  
1Department of Physiology, The University of Tokyo School of Medicine, Tokyo 1130033, Japan 

2Department of Radiology, NTT Medical Center Tokyo, Tokyo 1418625, Japan 

3Department of Mathematical Informatics, The University of Tokyo, Tokyo 1138656, Japan 

amasuda@mist.i.utokyo.ac.jp 
During rest, the human brain is not idle but shows a large amount of spontaneously fluctuating activity that is highly correlated between multiple brain regions^{1}^{, 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 restingstate networks (RSNs)^{3}^{, 4}^{, 5}^{, 6}. The RSNs, including the default mode network (DMN)^{5}^{, 6} and the frontoparietal network (FPN)^{7}^{, 8}, are highly reproducible across different healthy individuals^{9} and are considered to underlie fundamental and intrinsic functions such as selfreferential cognitive processes^{10}^{, 11}, maintenance of memory^{12} and attentional cognitive processes^{8}. These functions are considered to originate from complex functional interactions among the brain regions belonging to the RSNs^{13}. 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 quantified for the RSNs. As a result, it remains challenging to comprehensively understand this ensemble brain activity as an integrative largescale neural system^{14}^{, 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 secondorder model) seeks to fit a relatively simple binarystate model containing up to secondorder interaction terms (that is, firing rates and pairwise interaction) to empirical data of spatiotemporal spike trains. If the pairwise MEM is accurately fitted 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 higherorder 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, higherorder) MEM models. In fact, the pairwise MEM accurately describes firing patterns in the retinal tissues of primates recorded electrophysiologically in vitro^{16}^{, 17}^{, 18}^{, 19}^{, 20}^{, 21}, firing patterns and local field potentials (LFPs) in human cortical tissues in vitro^{22} and largescale firing patterns in the visual cortex of monkeys and cats in vivo^{23}^{, 24}. These findings 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 largescale 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 Pearson's correlation coefficient 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 Pearson's correlation coefficient in this paper. The FCbased analysis (based on Pearson's correlation coefficient) has successfully revealed that specific functional connexions are enhanced during specific cognitive processes^{25}^{, 26}, and the FC map serves as a promising indicator of several diseases^{27}^{, 28}, stages of development^{29}^{, 30} and levels of intelligence^{31}. However, because FC is calculated as Pearson's correlation coefficient between activities of pairs of regions^{3}^{, 4}, the FCbased method is founded on the implicit assumption that the pairwise interactions are independent of each other^{32}. If different pairwise interactions influence 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 region^{33}. Other studies also show that a significant FC value between a pair of regions does not distinguish a direct (that is, monosynaptic) connexion from an indirect (that is, polysynaptic) connexion^{34}^{, 35}. The FCbased methods may discard possibly rich information that would be revealed if we take into account that functional interactions influence each other^{8}.
Unlike the functional interactions estimated by the FCbased 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 FCbased method.
In the present study, we first quantify the complexity inherent in activities of an ensemble of brain regions during rest by fitting the MEM to restingstate 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 method^{36}^{, 37} and the mutual information (MI) method^{37}^{, 38}^{, 39}.
We applied the pairwise MEM to fMRI signals obtained from the brain regions belonging to the DMN^{5}^{, 6} and FPN^{7}^{, 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 functions^{15}^{, 40}. The coordinates of the brain regions were anatomically defined based on a metaanalysis of restingstate fMRI data obtained from 183 healthy participants^{41}^{, 42}. The DMN and FPN consisted of 12 and 11 regions, respectively (Table 1). We obtained ∼45 h of restingstate fMRI data from six healthy young adult participants (∼7.5 h per participant). Using the data, we extracted the lowfrequency signals (0.01–0.1 Hz) from all the brain regions included in the DMN or FPN (Fig. 1a), binarized the signals (Fig. 1b) and fitted the MEM to the binarized signals (Fig. 1c).
In Fig. 2a, 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), 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 figures. Figure 2a 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). 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 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.
We then examined the robustness of the accurate fitting of the pairwise MEM from three perspectives. First, we evaluated the accuracy of fit 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 socalled 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 (>99%) when the threshold lay approximately between −0.15 and 0.15 (Fig. 3b). Consistent with this result, the fitting of the pairwise MEM to the data was accurate for threshold values in this range (>80% for the DMN and >90% for the FPN; Fig. 3c). The threshold of −0.15, for example, implies that the accuracy of fit 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 fit 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 fit of the pairwise MEM (Fig. 3c).
Second, we confirmed 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 fitting linearly decreases with the number of regions N, one theory suggests that the accurate fitting may be an artifact of a small N^{43}. In this case, the validity of the pairwise MEM cannot be directly extrapolated to populations with a large N^{43}. For our data, the accuracy decreases sublinearly 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 fitted to data with N significantly larger than , where and δt are the averaged activation rate and the length of the time bin, respectively^{43}. In the present analysis, =0.045 s^{−1} and =0.041 s^{−1} for the DMN and FPN, respectively, and δt=9.045 s. Therefore, N_{c}≈2.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 N_{c}. 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. 1–3 and those obtained from participants no. 4–6. The interaction weights (that is, elements of the estimated interaction matrix, denoted by J_{ij} in equation (1); also see Fig. 1c) obtained from the pairwise MEM were similar between the two groups (Pearson's correlation coefficient r=0.96, P<10^{−3} in Fig. 3e; r=0.93, P<10^{−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.
If the functional interactions estimated by the pairwise MEM are closely relevant to some physiological 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 largescale anatomical connectivity determined from previous diffusion tensor imaging (DTI) of the brains of 80 healthy young human adults^{44}. 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 method^{36}^{, 37} and those based on the MI method^{37}^{, 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 five methods are shown in Fig. 4. For both the DMN (Fig. 4a) and FPN (Fig. 4b), the J_{ij} values estimated by the pairwise MEM are mainly large for pairs of brain regions with direct anatomical connexion, whereas, in the other methods, J_{ij} 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 J_{ij} 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 reflect either excitatory or inhibitory connectivity between pairs of regions.
To analyse the results more quantitatively, we compared the accuracy of the five methods in estimating the presence or absence of the anatomical connexions between pairs of regions. Figure 5a shows histograms of the values of interaction weights (that is, values of the entries of the matrices shown in Fig. 4) estimated by the five 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 coefficient) because they are conventionally used^{34}^{, 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 significantly larger than those for the anatomically unconnected pairs on average (P<0.001, twosample ttests in both the DMN and FPN; Fig. 5c).
The partial correlation and MI methods yielded similar significant differences between the anatomically connected pairs and unconnected pairs (for both models, P<0.001 in DMN, P<0.05 in FPN, twosample ttests; Fig. 5c. Also see Fig. 5a). In contrast, the FC and inverse Gaussian models did not predict anatomical connexions as accurately as the pairwise MEM (Fig. 5a). In the DMN, the interaction weights estimated by these two methods had larger values when region pairs were anatomically connected than unconnected (FC: P<0.05, absolute value of FC: P<0.01, inverse Gaussian model: P<0.01, twosample ttests). However, in the FPN, the difference was insignificant for both of the two methods (P>0.3).
For further comparison the five methods, we carried out a receiveroperating characteristic (ROC) analysis, which is sensitive to the difference in the distribution of the functional interactions between the anatomically connected and unconnected pairs (Fig. 6a). 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 classification 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 significantly larger than those for the other four methods (DMN: Z>7.3, P<1.0 × 10^{−11}, Bonferronicorrected; FPN: Z>9.4, P<1.0 × 10^{−16}, Bonferronicorrected; Fig. 6c). There was no significant difference between any other pair of the methods.
We binarized the fMRI signals to implement the pairwise MEM. However, fMRI signals originally represent continuous blood flows to brain regions, and the validity of the binarization is unclear. Therefore, we calculated the accuracy of fit of pairwise MEMs in which the signals were classified to three discrete levels, which we call the trinary MEM, by using a standard procedure described in Methods^{, 46}. When the signals recorded at the different brain regions and time points were evenly divided into three groups (Supplementary Fig. S1), the accuracy of fit was the highest. However, the accuracy of fit in this case was approximately equal to 0.55 (Supplementary Fig. S2) and significantly lower than those obtained for the binary pairwise MEM for both the DMN and FPN (P<10^{−3} in twosample ttests).
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 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 (P<0.05 in twosample ttests, Supplementary Fig. S3). However, ROC analysis revealed that the AUC for the trinary MEM was significantly smaller than that for the binary MEM (P<10^{−3} in twosample ttests, 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, 3^{N} versus 2^{N} data points). However, even the trinary MEM with N=4 and N=5, for which there are relatively few activity patterns (that is, 3^{4}=81 and 3^{5}=243 patterns), yields the accuracy of fit values <0.6. These accuracy of fit values are much smaller than those for the binary MEM with comparable numbers of activity patterns: the accuracy of fit was >0.9 for the binary MEM with N=8, with which there are 2^{8}=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.
We have demonstrated that the pairwise MEM accurately describes activity in the human RSNs. The model was fitted 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 Pearson's correlation coefficients, inverse Gaussian model, partial correlation method and MI method. These findings suggest that the largescale human brain networks during rest can be captured by a relatively simple secondorder model.
Our results extend the previous findings that the pairwise MEM accurately describes activity patterns observed in a slice or small parts of the brain^{16}^{, 17}^{, 18}^{, 19}^{, 22}^{, 23}^{, 24}. In particular, Tang and colleagues^{22} 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 bloodoxygenationleveldependent signals (that is, fMRI signals)^{47}, it is reasonable that the present study has successfully fitted the pairwise MEM to fMRI signals. Although finescale neural networks (<300 μm) seem to require more complex models, such as higherorder MEMs^{19}^{, 24}, in terms of the spatial scale, we have extended the applicability of the pairwise MEM method from microscopic neural populations (∼500 μm) to macroscopic populations at the level of the entire brain (∼10 cm).
In previous studies, large FC values have been suggested to imply the presence of direct anatomical connexion between regions^{2}^{, 34}^{, 45}. In particular, a study using highresolution fMRI conducted an ROC analysis similar to ours, and reported that the FC classifies 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 largescale brain networks to a better extent than the FCbased method does, although the pairwise MEM is a statistical model and not based on the physiology of the brain.
As secondorder interaction implies networks consisting of nodes and links, the pairwise MEM may serve as a tool for investigating largescale human brain networks. As previous studies did for networks generated by the restingstate FC (see a review^{14}), we can in principle calculate various graphtheoretical quantities for the RSNs defined by the pairwise MEM. This approach calls for a sufficient amount of restingstate fMRI data to estimate an interaction matrix among a relatively large number of brain regions and efficient algorithms to estimate the MEM for large N, such as those developed in previous studies^{19}^{, 48}^{, 49}. We can also apply the pairwise MEM to fMRI data measured during psychological tasks, though several modifications of the method may be in order. In previous fMRI studies, such taskspecific functional interactions were estimated by psychophysiological interaction analysis^{50}, Granger causality analysis^{51} and dynamic causal modelling^{52}. Among these methods, psychophysiological interaction analysis, like FCbased methods, is a seedbased 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 limited^{8}^{, 53}. Therefore, the pairwise MEM may serve as an alternative method for inferring taskspecific brain networks. In particular, recent theoretical developments of the MEM incorporate causality (signal flow from one region to another)^{54}^{, 55}, which will be an indispensable component for modelling taskrelated data.
Six healthy righthanded subjects (aged 20–23 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). T1weighted structural images were obtained for anatomical reference (resolution=0.81 × 0.81 × 1.20 mm^{3}). Functional imaging used gradientecho echoplanar sequences (TR=9.045 s, TE=35 ms, flip angle=90°, resolution=2 × 2 × 2 mm^{3}, 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 fixation point on the screen. One passive fixation session took ∼5 min, and each participant underwent 90 sessions (∼7.5 h) on 6–8 separate days. In each session, we discarded the first five images (9.045 s/volume × 5 volumes=45.225 s) to exclude the influence of transient processes before the equilibrium of longitudinal magnetization. In total, 17,820 volumes (2,970 volumes per participant) of restingstate functional images were obtained.
In line with previous studies of restingstate fMRI^{25}^{, 26}^{, 30}, we preprocessed the fMRI data as follows. First, for each session and each participant, the functional images were realigned, their slicetiming was corrected and they were normalized to the standard template image (ICBM 152) by SPM8 ( www.fil.ion.ucl.ac.uk/spm/). Second, for each session, after temporal bandpass filtering (0.01–0.1 Hz) with Matlab scripts written inhouse, 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, wholebrain 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 socalled ‘betavalue'^{56}^{, 57}^{, 58} except for a normalization factor.
We determined the coordinates of the brain regions belonging to the DMN and FPN based on a previous metaanalysis of restingstate fMRI studies^{30}^{, 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.
We fit the pairwise MEM to the restingstate fMRI data in essentially the same manner as did the previous studies of microscopic neural activities^{16}^{, 17}^{, 22}^{, 23}^{, 24}. 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 wholebrain 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 fit (Fig. 3b). The binarized activity of brain region i at discrete time t, denoted by σ_{i}^{t}, is either on (+1) or off (0). The network state (that is, pattern) at time t is given in vector form by 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 2^{N} possible network states. The empirical activation rate of region i, denoted by 〈σ_{i}〉, is given by . The empirical pairwise joint activation rate of regions i and j, denoted by 〈σ_{i}σ_{j}〉, is given by .
The pairwise MEM maximizes the entropy of the distribution of activity patterns under the restriction that 〈σ_{i}〉 and 〈σ_{i}σ_{j}〉 (1≤i≤j≤N) are preserved. It is known that such a distribution has the form , where P(V_{i}) is the probability of network state V_{i}, and
h_{i} represents the tendency of activation at region i, and J_{ij}, called the interaction weight, represents functional interaction between regions i and j. Positive and negative values of J_{ij} are interpreted as excitatory and inhibitory interactions, respectively. For an estimated probability distribution, the expected activation rate, 〈σ_{i}〉_{m}, and the expected pairwise joint activation rate, 〈σ_{i} σ_{j}〉_{m}, are given by and , respectively, where σ_{i}(V_{j}) indicates the activity (either +1 or 0) of σ_{i} under network state V_{j}.
We calculated h_{i} and J_{ij} by iteratively adjusting 〈σ_{i}〉_{m} and 〈σ_{i}σ_{j}〉_{m} toward 〈σ_{i}〉 and 〈σ_{i}σ_{j}〉, respectively, by using a gradient ascent algorithm. The iteration scheme is given by and , where α is the learning rate. We set α=0.75 to reduce computation time. We confirmed that the results were robust with respect to the choice of α ranging from 0.1 to 1.0.
As in previous studies^{16}^{, 17}^{, 22}^{, 23}^{, 24}, we evaluated the effectiveness of the pairwise MEM by using two information theoretic quantities. To define them, we also obtained the MEM without pairwise interaction between regions (that is, independent MEM) by determining h_{i} (1≤i≤N) under the condition that J_{ij}=0 (1≤i, j≤N). For the independent MEM, 〈σ_{i}σ_{j}〉_{m} is equal to 〈σ_{i}〉_{m}〈σ_{j}〉_{m}.
Using the resultant goodness of fit obtained for the two types of the MEMs, we define the accuracy index for the pairwise MEM by
where is the Kullback–Leibler 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 P_{N}. To evaluate reliability of the fitting, we calculate r_{S}=(S_{1}−S_{2})/(S_{1}−S_{N}), where is the entropy of the distribution of the network state in the k th order model. Note that S_{N} is the entropy of the empirical data. Using r_{S}, we define the estimation reliability by
The estimation reliability is equal to 1 if the values of h_{i} and J_{ij} are estimated without error^{59}.
As a large sample size is desirable for statistical fitting, we trained the pairwise MEM using the entire data and tested it against the same data. To exclude the possibility of over fitting, 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 fit by training the MEM using one subset and testing the estimated model against the other subset. (iii) We calculated the accuracy of fit 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 fit between the cross validation case (that is, (ii)) and the dependent trainingtesting case (that is, (iii)). We did not find significant differences between the two cases (cross validation versus dependent trainingtest cases: 0.67±0.0068 versus 0.68±0.0058 for the DMN, P>0.41; 0.76±0.0071 versus 0.76±0.0051 for the FPN, P>0.35, mean±s.d. in paired ttests). 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.
To examine the validity of binarizing fMRI signals, we estimated the accuracy of fit 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 σ_{i}^{t} is set to 1, 0 and −1 if the fMRI signal was larger than ɛ, between −ɛ and ɛ and smaller than −ɛ, respectively. Then, we fitted the trinary pairwise MEM in the almost same manner as that for fitting the binary pairwise MEM^{46}.
As a control of the functional interaction, we calculated the FC based on the Pearson's correlation coefficient 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 restingstate FC^{25}^{, 26}. After calculating the FC for each scanning session and each participant, we carried out Fisher's transformation^{28} and averaged the Zvaluebased FC across sessions and participants. As is consistent with a previous study^{60}, the FC matrices for the DMN and FPN were similar across participants as well as sessions for our data: Pearson's correlation coefficient between an arbitrary pair of the FC matrices of possibly different participants and sections was >0.81. The large correlation values suggest that the averaged FC is a good representation of the FC for a single participant and session.
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 , where X_{t}=[X_{1}(t), X_{2}(t),…X_{N}(t)] is the Ndimensional row vector representing the continuousvalued activity at the N regions at time t, is the row vector of the activity at the N regions averaged over sessions and participants, and is the N × N covariance matrix. In this model, we define the interaction weight between regions i and j as the (i, j) element of matrix Γ^{−1}.
As another control, we estimated functional interactions using the partial correlation method^{36}^{, 37}. Using the inverse of the covariance matrix, we defined the partial correlation between brain regions i and j, denoted by r_{ij}, by .
As another control, we estimated functional interactions using the MI method. In accordance with the previous studies^{37}^{, 38}^{, 39}, we first estimated the MI between regions i and j at frequency ω_{k} by , where mCoh_{ij} (ω_{k}) represents the multiple coherence between fMRI signals at regions i and j at frequency ω_{k}^{37}. We then obtained the MI between i and j, MI_{ij}, by averaging MI_{ij}(ω_{k}) over ω_{k} (0.01 Hz≤ω_{k}≤0.1 Hz).
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 adults^{42}. In the study^{44}, the authors obtained the DTI data in a Siemens Sonata 1.5T MRI scanner (Siemens Medical Systems, Germany) by using a twicerefocused singleshot EchoPlanar Imagingbased sequence (3 mm slice thickness with no interslice 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 classified 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 twosample ttest. 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 ttest. However, for the FC, we also submitted the raw FC values to the ttest, as was done in previous reports^{34}.
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 classification threshold. For an ROC curve, the AUC is equal to the area below the ROC curve and represents the classification accuracy of the estimation method^{34}. The deviation of the AUC was estimated based on a bootstrap method.
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.
How to cite this article: Watanabe, T. et al. A pairwise maximum entropy model accurately describes restingstate human brain networks. Nat. Commun. 4:1370 doi: 10.1038/ncomms2388 (2013).
Supplementary Figures S1S4.
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 GrantsinAid for Scientific 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 GrantinAid for Specially Promoted Research (19002010) to Y.M., a GrantinAid for Scientific Research B (22300134) to S.K. and a research grant from Takeda Science Foundation to Y.M.
References
Fox M. D., & Raichle M. E.,. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neurosci.8, 700–711 (Year: 2007) .17704812  
Sporns O.,. The human connectome: a complex network. Ann. N Y Acad. Sci.1224, 109–125 (Year: 2011) .21251014  
Friston K. J.,, Frith C. D.,, Liddle P. F., & Frackowiak R. S.,. Functional connectivity: the principalcomponent analysis of large (PET) data sets. J. Cereb. Blood Flow Metab.13, 5–14 (Year: 1993) .8417010  
Biswal B.,, Yetkin F. Z.,, Haughton V. M., & Hyde J. S.,. Functional connectivity in the motor cortex of resting human brain using echoplanar MRI. Magn. Reson. Med.34, 537–541 (Year: 1995) .8524021  
Raichle M. E.,et al.A default mode of brain function. Proc. Natl Acad. Sci. USA98, 676–682 (Year: 2001) .11209064  
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. USA100, 253–258 (Year: 2003) .12506194  
Beckmann C. F.,, DeLuca M.,, Devlin J. T., & Smith S. M.,. Investigations into restingstate connectivity using independent component analysis. Philos. Trans. R Soc. Lond. B Biol. Sci.360, 1001–1013 (Year: 2005) .16087444  
Cole D. M.,, Smith S. M., & Beckmann C. F.,. Advances and pitfalls in the analysis and interpretation of restingstate FMRI data. Front. Syst. Neurosci.4, 8 (Year: 2010) .20407579  
Damoiseaux J. S.,et al.Consistent restingstate networks across healthy subjects. Proc. Natl Acad. Sci. USA103, 13848–13853 (Year: 2006) .16945915  
Fair D. A.,et al.The maturing architecture of the brain's default network. Proc. Natl Acad. Sci. USA105, 4028–4032 (Year: 2008) .18322013  
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. USA106, 8719–8724 (Year: 2009) .19433790  
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, 4407–4420 (Year: 2011) .21430142  
Buckner R. L.,, AndrewsHanna J. R., & Schacter D. L.,. The brain's default network: anatomy, function, and relevance to disease. Ann. N Y Acad. Sci.1124, 1–38 (Year: 2008) .18400922  
Bullmore E., & Sporns O.,. Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci.10, 186–198 (Year: 2009) .19190637  
Deco G.,, Jirsa V. K., & Mcintosh A. R.,. Emerging concepts for the dynamical organization of restingstate activity in the brain. Nat. Rev. Neurosci.12, 43–56 (Year: 2011) .21170073  
Shlens J.,et al.The structure of multineuron firing patterns in primate retina. J. Neurosci.26, 8254–8266 (Year: 2006) .16899720  
Schneidman E.,, Berry M. J.,, Segev R., & Bialek W.,. Weak pairwise correlations imply strongly correlated network states in a neural population. Nature440, 1007–1012 (Year: 2006) .16625187  
Ganmor E.,, Segev R., & Schneidman E.,. The architecture of functional interaction networks in the retina. J. Neurosci.31, 3044–3054 (Year: 2011) .21414925  
Ganmor E.,, Segev R., & Schneidman E.,. Sparse loworder interaction network underlies a highly correlated and learnable neural population code. Proc. Natl Acad. Sci. USA108, 9679–9684 (Year: 2011) .21602497  
Schwartz G.,, Macke J.,, Amodei D.,, Tang H., & Berry M. J.,. Low error discrimination using a correlated population code. J. Neurophysiol.108, 1069–1088 (Year: 2012) .22539825  
Tkacik G.,, Prentice J. S.,, Balasubramanian V., & Schneidman E.,. Optimal population coding by noisy spiking neurons. Proc. Natl Acad. Sci. USA107, 14419–14424 (Year: 2010) .20660781  
Tang A.,et al.A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro. J. Neurosci.28, 505–518 (Year: 2008) .18184793  
Yu S.,, Huang D.,, Singer W., & Nikolic D.,. A small world of neuronal synchrony. Cerebral Cortex18, 2891–2901 (Year: 2008) .18400792  
Ohiorhenuan I. E.,et al.Sparse coding and highorder correlations in finescale cortical networks. Nature466, 617–621 (Year: 2010) .20601940  
Tambini A.,, Ketz N., & Davachi L.,. Enhanced brain correlations during rest are related to memory for recent experiences. Neuron65, 280–290 (Year: 2010) .20152133  
Chadick J. Z., & Gazzaley A.,. Differential coupling of visual cortex with default or frontalparietal network based on goals. Nat. Neurosci.14, 830–832 (Year: 2011) .21623362  
Stam C. J.,, Jones B. F.,, Nolte G.,, Breakspear M., & Scheltens P.,. Smallworld networks and functional connectivity in Alzheimer's disease. Cereb. Cortex17, 92–99 (Year: 2007) .16452642  
Dinstein I.,et al.Disrupted neural synchronization in toddlers with autism. Neuron70, 1218–1225 (Year: 2011) .21689606  
Supekar K.,, Musen M., & Menon V.,. Development of largescale functional brain networks in children. PLoS Biol.7, e1000157 (Year: 2009) .19621066  
Dosenbach N. U. F.,et al.Prediction of individual brain maturity using fMRI. Science329, 1358–1361 (Year: 2010) .20829489  
van den Heuvel M. P.,, Stam C. J.,, Kahn R. S., & Hulshoff Pol H. E.,. Efficiency of functional brain networks and intellectual performance. J. Neurosci.29, 7619–7624 (Year: 2009) .19515930  
Zuo X. N.,et al.Network centrality in the human functional connectome. Cereb. Cortex22, 1862–1875 (Year: 2012) .21968567  
Adachi Y.,et al.Functional connectivity between anatomically unconnected areas is shaped by collective networklevel effects in the macaque cortex. Cereb. Cortex22, 1586–1592 (Year: 2012) .21893683  
Honey C. J.,et al.Predicting human restingstate functional connectivity from structural connectivity. Proc. Natl Acad. Sci. USA106, 2035–2040 (Year: 2009) .19188601  
Lu J.,et al.Focal pontine lesions provide evidence that intrinsic functional connectivity reflects polysynaptic anatomical pathways. J. Neurosci.31, 15065–15071 (Year: 2011) .22016540  
Salvador R.,et al.Neurophysiological architecture of functional magnetic resonance images of human brain. Cereb. Cortex15, 1332–1342 (Year: 2005) .15635061  
Salvador R.,et al.A simple view of the brain through a frequencyspecific functional connectivity measure. NeuroImage39, 279–289 (Year: 2008) .17919927  
Salvador R.,, Anguera M.,, Gomar J. J.,, Bullmore E. T., & PomarolClotet E.,. Conditional mutual information maps as descriptors of net connectivity levels in the brain. Front. Neuroinform.4, 115 (Year: 2010) .21151357  
Lynall M.E.,et al.Functional connectivity and brain networks in schizophrenia. J. Neurosci.30, 9477–9487 (Year: 2010) .20631176  
Bullmore E.,et al.Generic aspects of complexity in brain imaging data and other biological systems. NeuroImage47, 1125–1134 (Year: 2009) .19460447  
Dosenbach N. U. F.,et al.A core system for the implementation of task sets. Neuron50, 799–812 (Year: 2006) .16731517  
Fair D. A.,et al.Functional brain networks develop from a ‘local to distributed' organization. Plos Comput. Biol.5, e1000381 (Year: 2009) .19412534  
Roudi Y.,, Nirenberg S., & Latham P. E.,. Pairwise maximum entropy models for studying large biological systems: when they can work and when they can't. Plos Comput. Biol.5, e1000380 (Year: 2009) .19424487  
Gong G.,et al.Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. Cereb. Cortex19, 524–536 (Year: 2009) .18567609  
Vincent J. L.,et al.Intrinsic functional architecture in the anaesthetized monkey brain. Nature447, 83–86 (Year: 2007) .17476267  
Mora T.,, Walczak A. M.,, Bialek W., & Callan C. G.,. Maximum entropy models for antibody diversity. Proc. Natl Acad. Sci. USA107, 5405–5410 (Year: 2010) .20212159  
Logothetis N. K.,, Pauls J.,, Augath M.,, Trinath T., & Oeltermann A.,. Neurophysiological investigation of the basis of the fMRI signal. Nature412, 150–157 (Year: 2001) .11449264  
Tanaka T.,. Meanfield theory of Boltzmann machine learning. Phys. Rev. E58, 2302–2310 (Year: 1998) .  
Hinton G. E.,. Training products of experts by minimizing contrastive divergence. Neural. Comput.14, 1771–1800 (Year: 2002) .12180402  
Friston K. J.,et al.Psychophysiological and modulatory interactions in neuroimaging. NeuroImage6, 218–229 (Year: 1997) .9344826  
Goebel R.,, Roebroeck A.,, Kim D.S., & Formisano E.,. Investigating directed cortical interactions in timeresolved fMRI data using vector autoregressive modeling and Granger causality mapping. Magn. Reson. Imaging21, 1251–1261 (Year: 2003) .14725933  
Friston K.,. Dynamic causal modelling. NeuroImage19, 1273–1302 (Year: 2003) .12948688  
Stephan K. E.,et al.Ten simple rules for dynamic causal modeling. NeuroImage49, 3099–3109 (Year: 2010) .19914382  
Marre O.,, Boustani, El S.,, Frégnac Y., & Destexhe A.,. Prediction of spatiotemporal patterns of neural activity from pairwise correlations. Phys. Rev. Lett.102, 138101 (Year: 2009) .19392405  
Roudi Y., & Hertz J.,. Mean field theory for nonequilibrium network reconstruction. Phys. Rev. Lett.106, 048702 (Year: 2011) .21405370  
de Pasquale F.,et al.A cortical core for dynamic integration of functional networks in the resting human brain. Neuron74, 753–764 (Year: 2012) .22632732  
Zhou H., & Desimone R.,. Featurebased attention in the frontal eye field and area V4 during visual search. Neuron70, 1205–1217 (Year: 2011) .21689605  
Badre D.,, Doll B. B.,, Long N. M., & Frank M. J.,. Rostrolateral prefrontal cortex and individual differences in uncertaintydriven exploration. Neuron73, 595–607 (Year: 2012) .22325209  
Yeh F. C.,et al.Maximum entropy approaches to living neural networks. Entropy12, 89–106 (Year: 2010) .  
van Dijk K. R. A.,et al.Intrinsic functional connectivity as a tool for human connectomics: theory, properties, and optimization. J. Neurophysiol.103, 297–321 (Year: 2010) .19889849 
Article Categories:

Previous Document: The mechanism of ultrafast structural switching in superionic copper (I) sulphide nanocrystals.
Next Document: A highmobility twodimensional electron gas at the spinel/perovskite interface of ?Al(2)O(3)/SrTiO...