Skip to main content

Main menu

  • HOME
  • CONTENT
    • Early Release
    • Featured
    • Current Issue
    • Issue Archive
    • Collections
    • Podcast
  • ALERTS
  • FOR AUTHORS
    • Information for Authors
    • Fees
    • Journal Clubs
    • eLetters
    • Submit
  • EDITORIAL BOARD
  • ABOUT
    • Overview
    • Advertise
    • For the Media
    • Rights and Permissions
    • Privacy Policy
    • Feedback
  • SUBSCRIBE

User menu

  • Log in
  • My Cart

Search

  • Advanced search
Journal of Neuroscience
  • Log in
  • My Cart
Journal of Neuroscience

Advanced Search

Submit a Manuscript
  • HOME
  • CONTENT
    • Early Release
    • Featured
    • Current Issue
    • Issue Archive
    • Collections
    • Podcast
  • ALERTS
  • FOR AUTHORS
    • Information for Authors
    • Fees
    • Journal Clubs
    • eLetters
    • Submit
  • EDITORIAL BOARD
  • ABOUT
    • Overview
    • Advertise
    • For the Media
    • Rights and Permissions
    • Privacy Policy
    • Feedback
  • SUBSCRIBE
PreviousNext
Research Articles, Systems/Circuits

Weak Higher-Order Interactions in Macroscopic Functional Networks of the Resting Brain

Xuhui Huang, Kaibin Xu, Congying Chu, Tianzi Jiang and Shan Yu
Journal of Neuroscience 25 October 2017, 37 (43) 10481-10497; DOI: https://doi.org/10.1523/JNEUROSCI.0451-17.2017
Xuhui Huang
1Brainnetome Center,
2National Laboratory of Pattern Recognition,
3Research Center for Brain-inspired Intelligence, Institute of Automation, and
4Center for Excellence in Brain Science and Intelligence Technology, Chinese Academy of Sciences, Beijing 100190, China,
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Xuhui Huang
Kaibin Xu
1Brainnetome Center,
2National Laboratory of Pattern Recognition,
4Center for Excellence in Brain Science and Intelligence Technology, Chinese Academy of Sciences, Beijing 100190, China,
5University of Chinese Academy of Sciences, Beijing 100049, China,
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Congying Chu
1Brainnetome Center,
2National Laboratory of Pattern Recognition,
4Center for Excellence in Brain Science and Intelligence Technology, Chinese Academy of Sciences, Beijing 100190, China,
5University of Chinese Academy of Sciences, Beijing 100049, China,
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
Tianzi Jiang
1Brainnetome Center,
2National Laboratory of Pattern Recognition,
4Center for Excellence in Brain Science and Intelligence Technology, Chinese Academy of Sciences, Beijing 100190, China,
6The Clinical Hospital of Chengdu Brain Science Institute, MOE Key Lab for Neuroinformation, University of Electronic Science and Technology of China, Chengdu 625014, China, and
7The Queensland Brain Institute, University of Queensland, Brisbane, Queensland 4072, Australia
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • ORCID record for Tianzi Jiang
Shan Yu
1Brainnetome Center,
2National Laboratory of Pattern Recognition,
4Center for Excellence in Brain Science and Intelligence Technology, Chinese Academy of Sciences, Beijing 100190, China,
5University of Chinese Academy of Sciences, Beijing 100049, China,
  • Find this author on Google Scholar
  • Find this author on PubMed
  • Search for this author on this site
  • Article
  • Figures & Data
  • Info & Metrics
  • eLetters
  • PDF
Loading

Abstract

Interactions among different brain regions are usually examined through functional connectivity (FC) analysis, which is exclusively based on measuring pairwise correlations in activities. However, interactions beyond the pairwise level, that is, higher-order interactions (HOIs), are vital in understanding the behavior of many complex systems. So far, whether HOIs exist among brain regions and how they can affect the brain's activities remains largely elusive. To address these issues, here, we analyzed blood oxygenation level-dependent (BOLD) signals recorded from six typical macroscopic functional networks of the brain in 100 human subjects (46 males and 54 females) during the resting state. Through examining the binarized BOLD signals, we found that HOIs within and across individual networks were both very weak regardless of the network size, topology, degree of spatial proximity, spatial scales, and whether the global signal was regressed. To investigate the potential mechanisms underlying the weak HOIs, we analyzed the dynamics of a network model and also found that HOIs were generally weak within a wide range of key parameters provided that the overall dynamic feature of the model was similar to the empirical data and it was operating close to a linear fluctuation regime. Our results suggest that weak HOI may be a general property of brain's macroscopic functional networks, which implies the dominance of pairwise interactions in shaping brain activities at such a scale and warrants the validity of widely used pairwise-based FC approaches.

SIGNIFICANCE STATEMENT To explain how activities of different brain areas are coordinated through interactions is essential to revealing the mechanisms underlying various brain functions. Traditionally, such an interaction structure is commonly studied using pairwise-based functional network analyses. It is unclear whether the interactions beyond the pairwise level (higher-order interactions or HOIs) play any role in this process. Here, we show that HOIs are generally weak in macroscopic brain networks. We also suggest a possible dynamical mechanism that may underlie this phenomenon. These results provide plausible explanation for the effectiveness of widely used pairwise-based approaches in analyzing brain networks. More importantly, it reveals a previously unknown, simple organization of the brain's macroscopic functional systems.

  • default mode network
  • frontoparietal network
  • functional connectivity
  • pairwise correlation
  • resting-state fMRI

Introduction

By using functional magnetic resonance imaging (fMRI), complex neuronal dynamics have been observed even when a subject is not engaged in specific tasks, reflecting that different brain areas interact with each other in the resting state (Biswal et al., 1995; Raichle et al., 2001; Greicius et al., 2003; Greicius et al., 2004; Beckmann et al., 2005; Fox et al., 2005; Fox and Raichle, 2007). During the past two decades, the most widely used method for studying these interactions is the analysis of functional connectivity (FC), which is defined as a pairwise correlation in BOLD activities observed from two spatially separated brain areas (Biswal et al., 1995; Greicius et al., 2003; Bullmore and Sporns, 2009). Based on FC, many functional systems of the brain, including the default mode network (DMN; Raichle et al., 2001; Greicius et al., 2003; Buckner et al., 2008) and the frontoparietal network (FPN; Daffner et al., 2003; Dosenbach et al., 2008; Zanto and Gazzaley, 2013), have been identified and studied intensively, providing important insights into the organization of macroscopic neuronal interactions in both normal and pathological conditions (Greicius et al., 2004; Greicius et al., 2007; Bullmore and Sporns, 2009; Brier et al., 2012). However, for a network with more than two nodes, there could be interactions among triplets, quadruplets, etc., of nodes termed higher-order interactions (HOIs). These HOIs can play important role for information processing [e.g., exclusive-OR (XOR) function] and cannot be revealed by pairwise-based analysis (Schneidman et al., 2003; Schneidman et al., 2006; Ganmor et al., 2011; Yu et al., 2011). Indeed, HOIs have been observed in a number of complex systems (Ohiorhenuan et al., 2010; Ganmor et al., 2011; Ohiorhenuan and Victor, 2011; Shemesh et al., 2013). More importantly, it has been reported that neuronal interactions at the mesoscopic level, reflected by spiking activities, exhibit significant HOIs (Ohiorhenuan et al., 2010; Ganmor et al., 2011). These results suggest the possibility that HOIs may also exist at the macroscopic level and may even play a role in orchestrating interactions among different brain areas. Nevertheless, so far, the HOIs among macroscopic brain networks remain largely unknown.

It is known that HOIs in a network are determined by a number of factors, including the system's intrinsic organization (e.g., whether it contains operation such as XOR) and other confounding influences such as the size (Roudi et al., 2009), topology, and node composition (Schneidman et al., 2006), as well as the inputs that drive the system (Ganmor et al., 2011; Macke et al., 2011). A recent study has claimed that pairwise interactions are good enough in characterizing the activities of two functional networks, the FPN and the DMN, under a specific condition: with the global signal (GS) regressed (Watanabe et al., 2013). However, many important questions remain. First, we do not know how these results can be extended to cover other functional systems under various conditions; for example, different network sizes, observation scales, spatial relations among nodes, and lack of GS regression (noGSR). Second, what are the dynamical mechanisms that may explain the HOI structures that can be observed from the empirical data? Answering these questions will provide important information regarding the organization of neuronal interactions at the macroscopic level. In the present study, we combine comprehensive HOI analysis based on both fMRI data recorded during the resting state and the simulation of a dynamic network model under various conditions to address these issues.

Materials and Methods

Dataset and preprocessing

The resting state fMRI (rs-fMRI) data analyzed here came from all 100 subjects (46 males and 54 females) of the Human Connectome Project (HCP) Q3 public data release (http://www.humanconnectome.org/data). During the resting-state scan, subjects underwent passive fixation in two separate sessions. In the two-fold cross-validation analyses, we separated the whole dataset into two halves, each one containing the data of one session from all subjects. Parameters for the resting data were as follows: TR = 720 ms, TE = 33.1 ms, FA = 52°, 2 × 2 × 2 mm voxels, FOV = 208 × 180 mm, and 72 oblique axial slices (Feinberg et al., 2010; Moeller et al., 2010; Setsompop et al., 2012; Van Essen et al., 2012). Each session contained 1200 time points for a total of 14.4 min. In total, 228,000 frames (approximately 46 h) of rs-fMRI data were used, with the first 60 frames of each scanning discarded.

The released fMRI data went through initial preprocessing by the HCP (Glasser et al., 2013). The detailed processing step adopted by HCP included: (1) gradient distortion correction, (2) motion correction, (3) distortion correction, (4) registration to the structural scan, (5) spline resampling to Montreal Neurological Institute (MNI) 2 × 2 × 2 mm3 space using the FSL software package (Jenkinson et al., 2002; Smith et al., 2004), and (6) intensity normalization to mean of 10,000 and bias field correction. We added the following preprocessing steps for our analysis: (1) band-pass filtering the time-series (0.01 Hz < f < 0.1 Hz); (2) censoring motion related frames with framewise displacement >0.2 mm (Ciric et al., 2017), performed in the analysis comparing the effect of different methods of motion artifacts rejection; and (3) regressing nuisance signals of 36 regressors described previously (Ciric et al., 2017). The regressors included tissue mean signal T (white matter mean signal, cerebrospinal fluid mean signal, global mean signal) and six head motion sequences R, their ([T, R]) derivatives, quadratic terms, and squares of derivatives. To examine the effect of global signal on higher-order interaction, we preprocessed the data separately with and without regressing the global mean signal (GSR vs noGSR). For the pipeline of noGSR, the four components of the GS were dropped (global mean, first derivatives of global mean, square of global mean, and square of the first derivatives).

We extracted mean region-of-interest (ROI) signals of six typical functional networks from all the brain regions (50 ROIs), which included the FPN (11 ROIs), the DMN (12 ROIs), the dorsal attention network (DAN, 8 ROIs), the executive control network (CON, 5 ROIs), the salience network (SAL, 7 ROIs), and the sensory-motor network (SMN, 7 ROIs). We defined each ROI by using a sphere with radius of 2, 5, and 10 mm, centering at the peak coordinate chosen based on previous studies (Brier et al., 2012; Watanabe et al., 2013; Zhan et al., 2016). All ROIs were masked by a gray matter probability map thresholded at 10% provided by the SPM12 software package (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/), with corresponding peak coordinates listed in Table 1.

View this table:
  • View inline
  • View popup
Table 1.

Coordinates of ROIs of six functional networks analyzed in the present study chosen according to Brier et al., 2012, Watanabe et al., 2013, and Zhan et al., 2016

To estimate HOIs in finer scale and spatial closer networks comparing with the 6 typical networks, we constructed 9-node networks, with all of their nodes (2 mm radius ROI for each one) completely contained within individual 10-mm-radius ROIs of 6 typical functional networks. For such networks, their nodes were centered at 8 vertices and geometric center of a Embedded Image × 10 mm3 cuboid that shared the same center with the 10-mm-radius ROI. In total, due to the constraint of the gray matter mask, 23 such finer-scale networks out of all 50 10-mm-radius ROIs could be constructed.

Data analysis

Outline of network analysis.

To reduce the dimensionality of the fMRI data to analyze HOIs properly, we dichotomized the continuous signals into binary time series. Specifically, we converted the continuous BOLD signals into discrete, binary time series of activity patterns by using different thresholds and different time bins (threshold ranges from 0.25 to 2 times the SD for individual ROIs; time bin width equals 1, 2, 4, 6, and 12 times of TR). The empirical pattern probability distribution was fitted with an independent model (Schneidman et al., 2006), a pairwise model (the Ising model; Schneidman et al., 2006; Tang et al., 2008; Yu et al., 2008), and a model with thresholding-induced HOIs [Macke et al., 2009; Yu et al., 2011; the dichotomized Gaussian (DG) model]. The strengths of different-order interactions were estimated directly from the distribution of experimental data and corresponding Ising and DG models (Yu et al., 2011). The fitting accuracy was quantified by Jensen–Shannon (JS) divergence DJS (Schneidman et al., 2006; Yu et al., 2011; see below for details).

Three interaction models used in the present study.

The independent model (Schneidman et al., 2006) is a maximal entropy model containing no interactions. Its pattern probability distribution is given by P(σ1, σ2, …, σn) = Embedded Image exp(∑ihiσi) where (σ1, σ2, …, σn) denotes the state pattern of a network. The state variable σi is the binary state (+1 or −1) of the region i. The intrinsic preference hi is the parameter representing the likelihood of individual nodes to be active, which is determined by the empirically observed activity rate < σi >. The partition function Z is a normalization factor.

The theoretical basis for understanding interaction structure in neural networks is information geometry, in which the complex interactions are decomposed into different orders (Nakahara and Amari, 2002). The second-order interactions were studied using the Ising model (Schneidman et al., 2006; Tang et al., 2008; Yu et al., 2008), which is a maximal entropy model containing pairwise interactions only and its expression is P(σ1, σ2, …, σn) = Embedded Image exp(∑ihiσi + ∑i<jJijσiσj). The intrinsic preference hi is the strength of the external field and the exchange interaction Jij is the strength of the second-order interaction between the regions i and j. The value of both hi and Jij were determined by the observed activity rate < σi > and pairwise correlation < σiσj >, respectively.

The DG model (Macke et al., 2009; Yu et al., 2011) is a two-layer model containing a bottom layer of continuous, multivariate Gaussian variables u and the top layer obtained by converting these Gaussian variables into binary variables by thresholding (i.e., dichotomization). The state variable u = (u1, u2, …, un) ∼ N(μ, δ) where n is the dimensionality and N(μ, δ) is a standard multivariate Gaussian distribution with the mean of μ and the covariance of δ. To make the binary variables fit the mean r and covariance ∑ of the experimentally observed data, the choice of μi and δij need to satisfy the relationship of μi = Φ−1(ri) and ∑ij = Φ2(μi, μj, δij) − Φ1(μi)Φ1(μj), respectively. The functions of φ1 and φ2 are the cumulative probability of one- and two-dimensional standard Gaussian distribution, respectively, and φ−1 is the inverse function of φ1. We applied the algorithm in Macke et al., 2009 to obtain the parameters μ and δ in the multivariate Gaussian distribution. Based on the distribution, the pattern probabilities P(σ1, σ2, …, σn) were calculated by integrating the multivariate Gaussian distribution from −∞ to the thresholding value or from the thresholding value to +∞.

Calculating strengths of different order interactions.

Strengths of different order interactions can be measured by the method used previously (Ganmor et al., 2011) when the pattern probability distribution P(σ1, σ2, …, σn) was known. Generally, the probability can be calculated as follows: P(σ1, σ2, …, σn)= Embedded Imageexp(∑i αi σi + ∑i < j βij σi σj + ∑i < j < k γijk σi σj σk + …) where α, β, and γ are strengths of the first-, second-, and third-order interactions (3OIs), respectively. Specifically, α1 = log P(100…) + log Z where P(100…) is the probability of the pattern in which the first node is active while all other nodes are quiescent, and Z = 1/P(000…); β12 = log P(1100…) + log Z − ∑i=12αi; and γ123 = log P(11100…) + log Z − ∑i=13αi − β12 − β13 − β23.

Time bin and thresholding methods for analyzing empirical data.

Different time bins (1, 2, 4, 6, and 12 times of TR) were used here to analyze the empirical data. In all presented results for the empirical data, the time bin was chosen to equal one time of TR. For cases with larger time bins, two methods were used to determinate the value for individual bins: the peak-value method and the mean-value method (the maximum value or the mean value among the data within a given time bin). Three different thresholding methods were tested here: (1) the absolute value method, which uses a predefined value as the threshold; (2) the global method, which uses the SD of all ROIs in a network to calculate the threshold; and (3) the scaled method, which uses the SD for individual ROIs to calculate the threshold. For the global and scaled methods, the threshold used was the summation of the mean and certain times of SD, as defined above. The global method was applied for all results of the empirical data shown here. Similar results can be obtained for the different time bin and thresholding methods (e.g., see Fig. 5C,D for the results obtained using the peak-value and mean-value thresholding methods). We note that the results presented in the main text were based on positive thresholding, that is,the dataset of positive deflections in the BOLD signals. The results based on various levels of negative thresholds were the same as that of positive thresholds.

Quantification of model performance.

To measure the accuracy of models' prediction to empirically observed pattern probabilities, we used the JS divergence (Schneidman et al., 2006; Yu et al., 2011). The JS divergence between probability distributions P and Q is defined by DJS(P, Q) = 0.5 × DKL(P‖M) + 0.5 × DKL(Q‖M), where M = (P + Q)/2 and DKL is the Kullback–Leibler divergence, defined as DKL(P‖Q) = ∑ipi log(pi/qi). We note that the JS divergence can quantify the distance between two distributions in a statistical sense, but it cannot address the functional implication of such distance. To quantify the improvement in performance of a model over the other, we define the relative accuracy index (RAI) as |DJSmodel1 − DJSmodel2|/DJSmodel1, which is adapted from the accuracy index used previously (Watanabe et al., 2013). In addition, to measure the higher-order information beyond the model, we quantified the entropy difference between various models and the empirical data as Imodel(H) = Hmodel − Hdata, where H is the entropy calculated by H = − ∑Pilog(Pi).

To determine whether a parametric model can provide accurate enough prediction for the data, we used a two-fold cross-validation method. The original dataset was separated into two equally sized sets. The Ising and the DG model fitted to the first set were used to predict the probabilities for the second set. Their prediction performance was compared with the “Exp” method, which uses the observed pattern probabilities of the first set to predict directly that of the second one. Previous study has shown that the performance of Exp method was highly sensitive to the amount of training data (Yu et al., 2011), so the two-fold cross-validation, which can provide the largest possible training dataset compared with other cross-validation methods, was used here to test the performance of the Ising and DG models in a more stringent way.

Examining statistical significance.

To determine whether the improvements of the DG model over the Ising and the Exp models in approximating both the activities of all six functional networks and internetwork activities were statistically significant, we separated the entire dataset randomly into two halves to perform the two-fold cross-validation test 100 times and, in each time, the performance of the three methods (Ising, DG, and Exp) was quantified by DJS, providing the distributions to assess the mean and variability of their performance. Statistical significance was then examined using the permutation test (Nichols and Holmes, 2002). Specifically, DJS of any two methods in all 100 times of two-fold cross-validation was first pooled together and then permuted and resampled (without replacement) 10,000 times, providing the permutation distribution of difference in DJS under the null hypothesis that the two models' performance was the same. This null hypothesis was rejected only if the actual difference between the mean DJS of these two models (averaged across 100 times of two-fold cross-validation) was larger than 95% of the cases in the permutation distribution. In addition, conventional statistical tests were conducted when enough independent samples could be obtained in two-fold cross-validations, such as in the cases for three-node subnetworks, finer-scale and spatially close networks, and 10-node networks with nodes randomly chosen from all 50 ROIs. In those cases, to examine the statistical significance, we applied the Kruskal–Wallis test, followed by the Mann–Whitney U test for post hoc analysis when it was appropriate. The significance level was set to be p < 0.05 for both the Kruskal–Wallis and Mann–Whitney U tests.

Method for estimating the significant HOIs on subnetwork state probabilities.

To estimate whether empirically observed 3OIs of specific subnetworks were significantly different from the DG model's prediction, we used the corresponding DG model to generate 1000 synthetic datasets, each with the same sample size of the experimental data. An empirically observed 3OI was considered significant if its value was outside the 95% confidence interval of 3OIs calculated from this synthetic dataset.

Data generated by three-neuron interaction model.

The data in Figure 1, B and C, were generated according to the probability distribution P as follows: Embedded Image The strengths of the second- and third-order interactions in the system were changed by adjusting the parameters βij and γijk, respectively. Mean correlation observed in the system was calculated by averaging Pearson correlation coefficients (PCCs) of all pairs of nodes.

Network model for generating simulated BOLD signals

Overview of the network model.

The network model consists of a number of “brain regions.” Each “region” is a collection of interacting neurons, with 80% excitatory and 20% inhibitory ones (Abeles, 1991; Brunel and Wang, 2001; Wang, 2002; Braitenberg and Schüz, 2013; Deco et al., 2014). The long-range connections between different regions were set according to the predefined connections and were exclusively excitatory. Within each region, neuronal activities were affected by the balance between the excitation and the inhibition, which was approximated by a slow-variable dominant dynamical equation using the mean-field method (Wong and Wang, 2006; Deco et al., 2013). In this system, the most important parameter that can influence the whole network's behavior is the strength of the long-range connections bridging different areas. For the sake of simplicity, we assumed the same strength for all interareal connections (termed global coupling). To simulate the fMRI signals, we then used the Balloon–Windkessel hemodynamic model that converts the fast neuronal dynamics unfolding at the neural network level to slower BOLD signals (Friston et al., 2003; Deco et al., 2013).

Mean-field neural networks.

In mean-field neural networks, the dynamics of each node, or brain region, could be described by the following nonlinear stochastic differential equations (Wong and Wang, 2006; Deco et al., 2013): Embedded Image Embedded Image Embedded Image where Si(t), H(xi), and zi(t) denote the average synaptic activity variable (the gating variable), the mean firing rate and the total inputs for the local region i, respectively. The connection matrix C represents the underlying structure of the network, with Cij = 1 if there is a connection from region j to I and Cij = 0 otherwise. The global coupling strength G controls the interarea interactions. The local coupling strength w contains the effects of both excitation and inhibition within individual areas. The transmission delay τij represents the time needed for activities traveling between areas i and j. The external input to the model, σξi(t), has two components: ξi(t) = ϵξcom(t) + (1 − ϵ)ξuni(t) where ξcom(t) and ξuni(t) are different standard Gaussian noise, whereas ξcom(t) is the common input to all areas and ξuni(t) is unique for individual areas. The parameter ϵ controls the ratio of common inputs. Parameters used for the activation function H were as follows: a = 270(VnC)−1, b = 108 Hz, and d = 0.154 s. The kinetic parameter γ = 0.641/1000 (the factor 1000 was used to convert the unit of time from second to millisecond) and the time constant τs = 100 ms. The synaptic conductance J = 0.2609 nA. The local coupling strength w was set to 0.85. The background input I0 = 0.3 nA. σvi(t) denotes the synaptic noise in the region i and ξi(t) is the standard Gaussian noise with the amplitude σ = 0.25. The transmission delay τij is usually set to be zero except for the control analysis, in which it was randomly chosen from the range of [0, 50] ms.

The Balloon–Windkessel hemodynamic model.

The Balloon–Windkessel hemodynamic model was used to translate the simulated neural synaptic signals into BOLD signals (Friston et al., 2003; Deco et al., 2013). The BOLD signal of a cortical region I was computed by the mean synaptic activity Si (Eq. 2) within this region. In short, for the ith region, the synaptic activity Si induces an increase in the vasodilatory signal xi (Eq.5). The vasodilatory signal xi is also affected by a leakage term (κixi in Eq. 5) and the feedback from the blood inflow fi, which is in turn proportional to xi itself (Eq. 6). The change in the blood volume vi is the linear function of fi and a leakage term (Eq. 7). The deoxyhemoglobin content qi is a nonlinear function of fi, vi, the resting oxygen extraction fraction ρ and itself (Eq. 8). Embedded Image Embedded Image Embedded Image Embedded Image Finally, the BOLD signal was calculated by the following: Embedded Image where V0 = 0.02, k1 = 7ρ, k2 = 2, and k3 = 2ρ − 0.2. The parameters used in the present study were set according to a previous study (Friston et al., 2003). It should be noted that, in this model, the BOLD signal of a region is the function of the synaptic activity within itself and is not affected by synaptic activities of other regions.

Network simulations with different topological structures.

Most of the simulation results were based on the structure of the FPN, with the exception of the case in which the cluster coefficient of networks was changed. To study the effects of inhibitory connections, we started with the FPN topology and randomly changed the directed excitatory connections to inhibitory ones with increasing probability, thereby gradually adjusting the ratio of inhibitory connections. To examine the effects of different cluster coefficients, we started with regular 11-node networks with mean degree < k > = 4, in which all nodes were arranged in a ring configuration and individual nodes were connected to two nearest neighbors from each side. Then, we rewired the directed connections randomly with an increasing probability (Watts and Strogatz, 1998). For each level of the cluster coefficient examined, six networks were generated.

Down-sampled data.

In our simulation, both the synaptic activities and the BOLD signals were down-sampled with an interval of TRs = 240 ms for the HOI analysis and the total time for each simulation was 60 h. Similar results were obtained from the data with different time bin widths; for example: 3, 6, 12, 24, and 36 times TRs, by using the global thresholding operation based on the peak value in the larger time bins.

Criteria of the network model constrained by the empirical data.

To ensure that the simulation can capture important features of the empirical rs-fMRI signals, two criteria were used: (1) the distribution of the BOLD signal amplitude should be unimodal and nearly symmetric; that is, a Gaussian-like distribution, which was the case for the empirical data; and (2) the overall level of pairwise correlations observed in simulations should cover the values exhibited in the empirical data.

Results

We first illustrated how HOIs can be important in determining a system's behavior and how they can introduce spurious pairwise correlations. We then analyzed the effects of HOIs in recordings of macroscopic brain dynamics based on fMRI data collected in 100 human subjects during the resting state. Finally, we constructed a two-layer model to simulate the generation of BOLD signals in distributed brain networks, examining the possible mechanisms underlying the behavior of HOIs that we observed in the empirical data.

Why are HOIs important?

Let us assume a small network composed of three nodes and each one could be either active or quiescent at a given moment. In Figure 1A, we depict one of such system (System 1) having three nodes that behave independently. In this case, there were quite a few instances in which three nodes were all active together by chance, representing a (trivial) synchronous state of the network. In another system (System 2), the probability of any nodes to be active (the first-order statistics), as well as the pairwise correlation between all pairs (the second-order statistics) were exactly the same as in System 1; that is, the correlations were all zero. Despite having the same first- and second-order statistics, System 2 exhibited no synchronous activities due to the injected third-order interaction (equivalent to the XOR operation) that actively prevented them. This simple example clearly illustrates that HOIs could be very important in shaping the systems' overall behavior and such an influence could escape our attention completely had the system been examined only at the pairwise level.

Figure 1.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 1.

Illustrating the effects of HOIs on system behavior. A, Activities of 2 3-node systems, with the same active probability (0.5) for each node and zero pairwise correlation between each pair of nodes. Synchronous active states (marked by red) observed in System 1 was absent in System 2 due to a specific form of 3OI (XOR). B, Mean correlation of a 3-node system is plotted as a function of 3OIs in the system, with various active probabilities (color coded), while all second-order interactions (2OIs) were fixed to be zero. C, Mean correlation of a 3-node system is plotted as a function of 3OIs in the system, with various strengths of 2OIs (color coded), whereas all active probabilities were fixed to be 0.04.

It was not just functionally important interactions such as the XOR operation in the above example that could be undetected at the pairwise level, HOIs can also “project” to the pairwise level, resulting in spurious pairwise correlations. In Figure 1, B and C, we show the results for two three-node systems, in which we kept the pairwise interactions stable while the 3OIs were increased systematically (see Materials and Methods for details). Despite the unchanged true pairwise interactions, the measured correlations between any two nodes were affected profoundly by the 3OIs. Importantly, there is no simple way to filter out such influence because the effects are largely determined by how active the individual nodes are (Fig. 1B) and how strong the true pairwise interactions are (Fig. 1C). In a large system composed of many nodes, there are many different HOIs; for example, fourth order, fifth order, etc., which can also “project” their effects to the pairwise level, adding yet another level of complexity.

In a system with abundant HOIs, the measured pairwise correlations are very difficult to interpret and the functional network extracted based on such measures should not be considered as a faithful representation of the true, underlying interaction structure.

HOIs in macroscopic functional networks

To determine whether HOIs are present in actual brain dynamics, we analyzed activities of six typical brain networks, including the FPN (containing 11 ROIs; Daffner et al., 2003; Watanabe et al., 2013; Zanto and Gazzaley, 2013), the DMN (containing 12 ROIs; Raichle et al., 2001; Greicius et al., 2003; Watanabe et al., 2013), the DAN (containing 8 ROIs; Fox et al., 2006; Brier et al., 2012; Zhan et al., 2016), the CON (containing 5 ROIs; Seeley et al., 2007; Brier et al., 2012; Zhan et al., 2016), the SAL (containing 7 ROIs; Zhou et al., 2010; Zielinski et al., 2010; Brier et al., 2012; Zhan et al., 2016), and the SMN (containing 7 ROIs; Fox et al., 2009; Brier et al., 2012; Zhan et al., 2016), in the rs-fMRI data of 100 human subjects collected in HCP (see Materials and Methods for details). Below, we first present the results in details for the FPN and the DMN and then show more comprehensive results for all six typical functional networks and other networks with various parameters.

Effects of HOIs in the FPN and DMN

As shown in Figure 2A, the BOLD signals recorded from all ROIs within individual networks were first converted to binary activity patterns using a thresholding operation, with 1 indicating superthreshold episodes (active) and −1 indicating subthreshold episodes (quiescent). Next, we measured the probabilities of all possible patterns of the networks, which provide important information about the system's behavior. For example, for FPN (with 11 nodes), we measured the probabilities for all 211 different activity patterns. We then tried to predict/explain all these probabilities using various parametric models, which take into account different orders of statistics measured from the empirical data. These models included: (1) the independent model, which introduces the same first-order statistics (how active a node is) as empirically observed but assumes that all nodes behave independently (Schneidman et al., 2006; Yu et al., 2011); (2) the pairwise maximal entropy model or Ising model, which introduces the same first- and second-order (pairwise correlations) statistics as empirically observed but assumes all HOIs are zero (Schneidman et al., 2006; Tang et al., 2008; Yu et al., 2011); and (3) the DG model. In addition to the observed first- and second-order statistics, the DG model also takes into account a specific structure of HOIs due to the thresholding operation used in extracting the binary activity patterns (Amari et al., 2003; Macke et al., 2009; Yu et al., 2011). To determine how well those models can predict/explain the probabilities of all patterns, we used the JS divergence (DJS; see Materials and Methods for details). Comparing different models' performances gave us a quantitative measure of the relative importance of interactions of different orders in determining the networks' behavior.

Figure 2.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 2.

Performance in approximating fMRI data with different interaction models. A, Resting-state BOLD signals were extracted from ROIs in typical networks (e.g., the FPN and the DMN) and converted to binary patterns through a thresholding operation. B, C, Performance of different models, including the independent model (“Ind”), the Ising model, and the DG model, in approximating network activities of the FPN (B) and the DMN (C) under the GSR condition. Observed pattern probabilities are plotted against model predictions. Each data point represents one pattern. Solid line indicates the equality. For visual clarity, patterns with probabilities <0.001 are not shown. D, E, Same analyses as shown in B and C but under the noGSR condition for the FPN (D) and the DMN (E). In B–E, the BOLD signals used were based on 10-mm-radius ROIs and the threshold used here was the mean plus 1 SD. The prediction accuracy was quantified by the JS divergence (DJSmodel), which is specified for different models in the figure.

In fMRI signal analysis, GSR is a widely adopted preprocessing step (Fox et al., 2009; Murphy et al., 2009; Van Dijk et al., 2012). However, the GS contains important information about neuronal dynamics and can be used as a biomarker in monitoring brain diseases (Schölvinck et al., 2010; Yang et al., 2014). Importantly, it may also introduce strong HOIs to the system and drastically change the overall interactions structure (Macke et al., 2011; Yu et al., 2011). Therefore, in our HOI analysis, we examined the GSR and noGSR conditions separately.

In Figure 2, B and C, we show results for the FPN and the DMN, respectively, under the GSR condition. Clearly, if the probabilities were predicted using the first-order statistics only (the independent model, indicated by black dots in Fig. 2B,C), then the prediction had enormous errors. If the second-order statistics were included into the prediction (the Ising model, indicated by blue asterisks in Fig. 2B,C), then the accuracy was greatly improved. Measured by the RAI, defined as |DJSmodel1 − DJSmodel2|/DJSmodel1; see Materials and Methods for details; Watanabe et al., 2013), the prediction was improved by 96.35% and 93.1% for the FPN and the DMN, respectively, compared with the independent model. Moreover, if the threshold-induced HOIs were taken into account (the DG model, indicated by red diamonds in Fig. 2B,C), RAI could be further improved (30.9% and 32.4%) compared with the Ising model.

In Figure 2, D and E, we show the results for the same analyses but under the condition of noGSR. We found that the improvements in RAI achieved by the Ising model (over the independent model) were quite similar to the GSR condition (94.9% and 91.5% for the FPN and the DMN, respectively). However, in contrast to the GSR condition, now for the patterns that were mostly likely to occur, the Ising model gave much worse results, as indicted by large deviation of those data points from the equality line, suggesting an increased amount of HOI under such condition. Importantly, the DG model's prediction in the noGSR condition was as accurate as that in the GSR condition, resulting in the RAI increasing by >60% (67.6% and 64.9% for the FPN and the DMN, respectively) compared with the Ising model. These results indicated that the apparent HOIs in the noGSR condition can be satisfactorily explained by considering the thresholding operation with the DG model.

The examples analyzed above suggest that, although the Ising model gives good prediction for the GSR condition, a specific structure of HOIs due to thresholding are needed to explain the systems' behavior for the noGSR condition and the RAI as a measure to judge whether a model can give accurate prediction to the empirical data has important limitations. Therefore, we next performed more comprehensive analyses using three different measures and under various thresholding levels to further examine the importance of pairwise and HOIs for these two networks.

First, consistent with the above results, we found that the JS divergence, which is another way to measure the prediction accuracy, of the Ising model's prediction differed significantly under the GSR and noGSR conditions, whereas that of the DG model's prediction kept almost identical between the two conditions (Fig. 3A,B), suggesting that the GS does not introduce true HOIs. Second, to examine directly the information encoded by HOIs, we defined the higher-order information beyond a model as the entropy difference between the system and that model (see Materials and Methods for details). In Figure 3, C and D, we show that the higher-order information beyond the Ising model was much larger under noGSR condition. However, the higher-order information beyond the DG model was very similar between these two conditions, confirming that there was actually no information encoded by HOIs in the noGSR condition. Finally, we determined the strength of 3OIs directly in the systems and compared that with those predicted by the DG model. We found that, although much stronger 3OIs were observed empirically under the noGSR condition, this increase could be explained accurately by the DG model (Fig. 3E,F). Again, these results suggest that, in the noGSR condition, observed HOIs can largely be explained by the thresholding operation.

Figure 3.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 3.

Three quantitative measures of the HOIs in the data with and without GSR. The JS divergences in predicting pattern probabilities (A, B) and the higher order information beyond models (C, D) are plotted as a function of the thresholding level for the FPN (A, C) and the DMN (B,D). Results are shown for both the Ising (blue stars) and the DG (red diamonds) models. E, F, Mean strengths of 3OIs observed in the empirical data (black circles) and corresponding DG model (red diamonds) plotted as a function of the thresholding level for the FPN (E) and the DMN (F). Dotted lines indicate GSR; solid lines, noGSR. Note that both the JS divergence and higher-order information beyond models decrease with the threshold is due to fewer active patterns occurred with higher thresholds and the mean strengths of 3OIs were calculated based on actually occurred activity patterns. The binary data used here were extracted based on 10-mm-radius ROIs.

Comprehensive analyses of HOIs for large-scale functional networks

In addition, results similar to Figure 3 were obtained for other four networks, including the DAN, CON, SAL, and SMN. In Figure 4, we first show the diversity of the PCC of BOLD signals across all six networks with and without the GSR (Fig. 4A,B). To determine whether such a correlation structure can be well captured by studying the thresholded signals, we compared the correlation coefficients between all pairs of continuous BOLD signals with their counterparts calculated based on dichotomized signals. These two measurements strongly correlated with each other (r > 0.93, p < 0.01 for both GSR and noGSR data). In addition, the negative correlations in the continuous signals were mostly preserved after the binarization (Fig. 4C,D). These results suggest that the dichotomized signals indeed reflect the correlation structure of the original data. Next, we summarized the overall results concerning HOIs in all six networks using the two-fold cross-validation method. Specifically, the original dataset was separated into two equally sized sets. The Ising and the DG model fitted to the first set were used to predict the probabilities for the second set. In addition, we also used the observed pattern probabilities of the first set to predict directly that of the second one (Exp method). We found that, in both the GSR (Fig. 4E) and noGSR (Fig. 4F) conditions, the prediction of the DG model, measured by the JS divergence, was better than that of the Exp method for all six networks (p < 0.05, permutation test, see Materials and Methods for details). This is possible because, with the limited sampling (the first half of data), the DG model needs to estimate a much smaller number of parameters (only the first- and second-order statistics for a network with N nodes; ∼N2 parameters) compared with the naive model (Exp, ∼2N parameters). Because one half of the dataset equals ∼23 h recording, our results indicate that the DG model, with ∼N2 parameters, could characterize brain networks' behavior very accurately. Our results from all six networks suggest that the lack of HOIs may be a fundamental feature for all major functional networks of the resting brain.

Figure 4.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 4.

Comprehensive analyses on strength of HOIs in brain networks. A, B, Correlation matrix among 50 brain regions for the GSR and noGSR conditions. Among them, six major networks are marked by colored rectangles. C, D, Comparison between the PCCs measured based on continuous BOLD signals among all 50 ROIs in the empirical data (x-axis) and their binarized counterpart (y-axis) in both the GSR (C) and noGSR (D) conditions. E, F, Prediction accuracy of the Ising and the DG models for all six typical networks under the GSR (left) and noGSR (right) conditions by using the two-fold cross-validation method. The ROIs of six networks (color coded) are marked in the horizontal view of the brain (inset). Note that the differences in absolute values of divergence were due to different sizes of the networks. G, H, To test the effects of HOIs for internetwork interactions, the same analysis as in E and F but for five 10-node networks (Rand-1 to Rand-5) with nodes chosen randomly from all 50 ROIs (without replacement) under the GSR (G) and noGSR (H) conditions was conducted. In E–H, 100 two-fold cross-validation based on randomly separating all 100 subjects into two groups (50 subjects each) were performed. The data are presented as mean ± SD. Note that all pairwise comparisons are significant unless they are marked as N.S. I, J, Probability distribution of active size, defined as number of active regions within the same time bin, for the whole network consisted of 50 ROIs from all six typical networks, under the conditions of the GSR (I) and noGSR (J). Solid lines indicate distribution estimated directly from the empirical data; dashed lines, distribution predicted by the DG model. The BOLD signals used here were based on 10-mm-radius ROIs and the threshold used here was the mean plus 1 SD unless marked otherwise.

Because there were significant PCCs between the regions belonging to different networks (Fig. 4A,B), we also analyzed the effect of HOIs among regions belonging to different networks with the same method. Specifically, we studied the HOIs in five networks of 10 nodes, which were chosen randomly from all 50 ROIs (without replacement). We found that, in both the GSR and noGSR conditions, the prediction of the DG model was better than that of the Exp method (Fig. 4G,H; p < 0.05, permutation test), confirming that HOIs are also weak for interactions among different networks. Similarly, significantly better performance of the DG model can be demonstrated for 100 random 10-node networks, with nodes chosen randomly from all 50 ROIs (with replacement; Fig. 5A; p < 0.05, Kruskal–Wallis and Mann–Whitney U tests, see Materials and Methods for details). Furthermore, we found that the DG model could predict activities correctly for the entire network consisting of all 50 ROIs from these six networks (Fig. 4I,J), suggesting that weak HOIs hold true even for a global network of the resting brain.

Figure 5.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 5.

Robustness of DG model's performance in various conditions examined by the two-fold cross-validation method in the GSR and noGSR conditions tested in 100 10-node networks with nodes randomly selected (with replacement) from all 50 ROIs (A), across different ROI sizes (B), across different time bin widths with two thresholding methods (C, D), and with a different preprocessing method to remove motion artifacts (E) (cf. Fig. 4E,F). In A–E, data are presented as mean ± SD. All corresponding pairwise comparisons are significant except for one case marked as NS in A. The threshold used here was the mean plus 1 SD. In A and C–E, the BOLD signals used were based on 10-mm-radius ROIs. In B–D, the results for the FPN are shown as examples.

Because the current study was aimed primarily at examining the functional networks at the macroscopic level, in the results mentioned above, we focused on relatively large spatial scale (10-mm-radius ROIs) and a specific temporal scale (TR = 720 ms as the time bin width). First, to determine whether, in a finer scale of ROIs, the HOIs would be more pronounced, we performed the same analyses for both the PFN and DMN with 2- and 5-mm-radius ROIs. The results obtained at these two finer scales are very similar to the results shown in Figure 4, C and D, demonstrating the significantly better performance of the DG model compared with the Exp method (p < 0.05, permutation test; results of the FPN are shown in Fig. 5B). Second, for different temporal scales, similar conclusions could be drawn for data obtained with dramatically different time bin widths (p < 0.05, permutation test; results of the FPN are shown in Fig. 5C,D). These results indicate that the lack of strong HOIs is not dependent on specific spatial or temporal scale at which the BOLD signals were analyzed during the resting state.

Because the motion artifacts had been demonstrated to be important in analyzing brain networks' activities based on fMRI signals (Power et al., 2012; Van Dijk et al., 2012; Power et al., 2014), we scrutinized the data using a more stringent method for motion artifact rejection suggested in a recent study (Ciric et al., 2017). The results were highly consistent with the previous ones. We found that the performance of the DG model was significantly better than the Exp method in all 6 functional networks under both the GSR and noGSR conditions (Fig. 5E, p < 0.05, permutation test), indicating that our conclusions were not affected by motion artifacts in the original BOLD signals.

Effects of HOIs on subnetworks

It remains to be tested whether HOIs are strong at specific subnetworks; for example, among three nodes, within these typical functional networks. To address this issue, using the FPN and the DMN as examples, we examined all three-node subnetworks within each of them and compared the strengths of observed 3OIs with that predicted by the DG model. Because the networks extracted in the noGSR condition exhibit stronger HOIs, here, we focus the analysis on this condition. We found that ∼65% of 3OIs in the real data were not significantly different from the DG model's prediction, judged as its value within 95% confidence interval of 1000 length-matched samples generated by the DG model (see Materials and Methods for details). For the remaining ∼35% of 3OIs that could not be fully explained by the DG model (referred to as significant 3OIs hereafter), we found that they had very limited influence on corresponding three-node subnetworks' activities. To demonstrate this, in Fig. 6, A and B, we plot the observed versus predicted (by the DG model) pattern probabilities for all triplets chosen within the FPN and the DMN. For the subnetworks composed of the triplets with significant 3OIs, we examined how these 3OIs were important in shaping subnetwork's activities, by using the same two-fold cross-validation analysis. We did not find larger prediction error for subnetworks containing significant 3OIs. Furthermore, the DG model was still significantly superior compared with the Exp method for this group of subnetworks (Fig. 6A,B, insets; p < 0.05, Kruskal–Wallis and Mann–Whitney U tests), resembling the situation for subnetworks without significant 3OIs. These results demonstrated that, even for the subnetworks with significant HOIs, the DG model's prediction is accurate enough so that the error is smaller than statistical fluctuations in tens of hours of recording (reflected by the Exp method). Therefore, similar to the whole network, the HOIs that cannot be explained by the DG model play a very minor role, if any, in shaping activities at the subnetwork level. Very similar results were obtained when we investigated the 3OIs for subnetworks consisted of nodes belonging to different networks (Fig. 6C,D). Taken together, we did not find significant HOIs that are important in shaping activities of subnetworks at both the intranetwork and internetwork levels.

Figure 6.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 6.

Effects of HOIs on subnetworks under the noGSR condition. A, B, Observed 3OIs of 3-node subnetworks of the FPN (A) and the DMN (B) are plotted against the 3OIs of corresponding DG models. C, D, Same analysis as in A and B applied to subnetworks within the first two random networks analyzed in Figure 4H (C: Rand-1; D: Rand-2), showing the behavior of 3OIs in subnetworks with nodes belonging to different brain networks. The significant and nonsignificant 3OIs are marked with red crosses and blue circles, respectively. The dashed line represents the equality. Insets, Performance of three methods in predicting the activities of all 3-node networks (3NWs) that contains significant 3OIs was quantified by the two-fold cross-validation method. The bar diagrams represent mean ± SD of DJS. Note that all pairwise comparisons are significant unless they are marked as N.S. The binary data for the four networks used here were the same as those in Fig. 4.

Effects of spatial vicinity on HOIs

At the mesoscopic scale, that is, the networks composed of individual neurons, it has been reported that, the closer the neurons are to each other in space, the stronger the HOIs that will be found among them (Ohiorhenuan et al., 2010). To determine how spatial proximity can affect the HOIs among macroscopic functional networks, we analyzed the HOIs among networks at finer scale and in which individual nodes were much closer to each other. Specifically, we used 2 mm as the radius to define ROIs and chose 9 nodes within a single coarse-grained DMN or PFN node (10-mm-radius ROI). In total, 23 of such finer-scale networks were studied (see Materials and Methods for details). The results were highly consistent with what we found for spatially separated networks because the DG model performed significantly better than the Ising model across all threshold levels examined (Fig. 7A,B; p < 0.05, Kruskal–Wallis and Mann–Whitney U tests). Although there were apparently sizable HOIs exhibited in both the GSR and noGSR conditions, they could be satisfactorily explained by the DG model (Fig. 7C). To further determine whether the DG model's approximation of the network activities was accurate enough, we applied the same two-fold cross-validation test to these 23 finer scale networks. The results confirmed that the DG model was better than the Exp method in both the GSR and noGSR conditions (Fig. 7D, p < 0.05, Kruskal–Wallis and Mann–Whitney U tests). Therefore, the weak HOIs seem to be a general feature for functional brain networks extracted by BOLD signals at rest regardless of the spatial relationship among the constituent network nodes.

Figure 7.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 7.

HOIs in finer scale, spatial close networks. A–C, Three quantitative measures of the HOIs in networks with individual nodes in spatial neighborhood at a finer scale; 2 mm was used as the radius to define ROIs. Nine-node networks within individual coarse-grained ROIs (10-mm-radius ROI) were chosen to form a network. The results for all observable 23 such finer-scale, spatially close subnetworks are shown (mean ± SD). The analyses and meaning of markers are the same as in Fig. 3. D, Two-fold cross-validation method was used to examine the three methods' performance. The bar diagrams represent mean ± SD of DJS for individual methods across the 23 networks in the GSR and noGSR conditions. The threshold used here was the mean plus 1 SD. In A–D, all corresponding pairwise comparisons are significant except for one case marked as N.S. in D.

Behavior of HOIs in dynamic network models

Overview of the model

To better understand our empirical results, including the largely absent HOIs and the difference between the GSR and noGSR condition, we built a network model that covered the condition similar to the empirical data (see Materials and Methods for details) but extended to a much wider range of conditions to examine the behavior of HOIs thoroughly. In the model with predefined network structure (e.g., the structure of the FPN based on Gong et al., 2009 and Watanabe et al., 2013 as shown in Fig. 8A, top), the synaptic dynamics were generated by a mean-field neural network model (Wong and Wang, 2006; Deco et al., 2013) and then the Balloon–Windkessel hemodynamic model was used to convert the fast synaptic dynamics (Fig. 8A, middle) to slower BOLD signals (Fig. 8A, bottom) (Friston et al., 2003; Deco et al., 2013; see Materials and Methods for details).

Figure 8.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 8.

HOIs detected in a dynamic network model. A, Diagram showing the model with the same number of nodes and connectivity as the FPN (Gong et al., 2009; Watanabe et al., 2013) (top), synaptic activities averaged within individual nodes (Wong and Wang, 2006; Deco et al., 2013) (middle), and corresponding BOLD signals (Friston et al., 2003; Deco et al., 2013) (bottom). B, C, The mean strengths of 3OIs in simulated synaptic activities (B) as well as BOLD signals (C) are plotted as a function of global coupling, whereas no common input was injected into the system. D, E, Mean strengths of 3OIs in simulated synaptic activities (D) as well as BOLD signals (E) are plotted as a function of common input ratio, whereas global coupling was set to 0.1. F, G, Similar analyses for the simulated data as in D and E with GSR before calculating 3OIs. Different thresholds are color coded. Solid line indicates the results for the empirical data; dotted lines, results for the DG model.

Effects of global coupling and common inputs on HOIs

In this system, the connection strength between different regions (termed global coupling, assumed to be homogeneous) plays an important role in shaping the whole network's behavior. In Figure 8, B and C, we show the results of the HOIs detected in the model as a function of global coupling, with the same topology of interarea connections as shown in Figure 8A. We found that the mean strengths of 3OIs in simulated synaptic activities and BOLD signals were weak. Similar to the empirical results, such an absence of HOIs could be correctly approximated by the DG model. Importantly, the mean strength of 3OIs was not affected by the increase of global coupling (Fig. 8B,C) from 0 to 0.1, which was chosen to make the average pairwise correlation in the simulated signals cover the range of that in the empirical data (see Materials and Methods for details). It is noteworthy that, based on the comparison between the results of simulated synaptic activities and BOLD signals, we can also conclude that the absence of HOIs is not due to the lower temporal resolution of the BOLD signal. Instead, it is a property rooted in how neuronal populations among different brain areas interact with each other.

Based on the empirical data, we demonstrated that, under the noGSR condition, apparent HOIs are much stronger. The nature and origin of the GS in the brain are not completely understood yet, but the common inputs mediated through diffuse ascending projections from subcortical structures, including the locus ceruleus (Sara, 2009), have been suggested as possible sources. Therefore, to understand how global signal gives rise to apparently stronger HOIs, we introduced global signal in the model by injecting common inputs to the entire network. Specifically, we systemically changed the ratio of common input to all external inputs, mimicking the effects of global signal in the system. We found that, in both synaptic activities and BOLD signals, strong common inputs led to strong 3OIs (Fig. 8D,E). These results explain both weak HOIs in data with GSR and apparently strong HOIs in data without GSR. Moreover, consistent with our empirical results, the HOIs induced by common inputs could be well explained by the DG model. The mean strengths of 3OIs in empirical data can be largely predicted by the DG model for both synaptic activities and BOLD signals, indicating that they were not true HOIs intrinsic to the system. To further mimic the situation of GSR in the empirical data, we also preprocessed the simulated data with GSR similar to that was done in the empirical data. We found that at various levels of common inputs, both simulated synaptic activities and BOLD signals showed no significant HOIs that cannot be explained by the DG model, resembling the empirical results with the GSR condition (Fig. 8F,G).

Similar dynamic features between the empirical and simulated BOLD signals

To determine whether the parameter space explored above covered the network states having similar dynamic features with the real BOLD signals, we compared the strength of pairwise correlations between the empirical and simulated data. We found that, with the proper choice of global coupling and common input ratio, the model exhibited distribution of correlation strength highly similar to the empirical data (Fig. 9). Specifically, when the global coupling was chosen to be 0.1 (without common inputs), the distribution of PCC (Fig. 9A) and mean magnitude squared coherence (MSC, averaged across the frequency range of 0.01–0.1 Hz; Sun et al., 2004; Fig. 9C) observed in the simulation was similar to that of the empirical data with GSR. When the common input ratio was chosen to be 0.3, the distribution of PCC (Fig. 9B) and mean MSC (Fig. 9D) observed in the model was similar to that of the empirical data without GSR. These results suggested that the model indeed captured essential features of brain dynamics and therefore can be used to study the mechanisms underlying the lack of strong HOIs observed in the empirical data.

Figure 9.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 9.

Similarity in statistical features between the empirical and simulated data. A, B, Comparison between the distribution of PCCs measured for all ROIs of the FPN in the empirical data and in the simulation with the global coupling = 0.1 (see figure for details in parameter setting; CI, common input ratio), in both the GSR (A) and noGSR (B) conditions. C, D, Comparison between the distribution of mean MSC averaged across the frequency range of 0.01–0.1 Hz measured for all ROIs of the FPN in the empirical data band-pass filtered between 0.01 and 0.1 Hz and in the simulation with the global coupling = 0.1 in both the GSR (C) and noGSR (D) conditions. We note that other correlation measures, such as the Spearman correlation coefficient, gave essentially the same results.

Robustness of weak HOIs in various conditions

The lack of strong HOIs in the model was highly robust in both the GSR and noGSR conditions within a wide range of other key parameters such as ratio of inhibitory connections of the network (Fig. 10A), transmission delays among different regions (Fig. 10B), and time scale used for analysis (Fig. 10C). Specifically, we investigated how HOIs can be affected by the typological feature (e.g., the clustering coefficient) of the underlying network. We found that changing the clustering coefficient from 0.2 to 0.6 in the simulated networks (see Materials and Methods for details) did not introduce significant HOIs that cannot be explained by the DG model (Fig. 10D). Consistently, the DG model significantly outperformed the Exp method in the two-fold cross-validation regardless of the networks' clustering coefficient (p < 0.05, Mann–Whitney U tests; Fig. 10D, inset).

Figure 10.
  • Download figure
  • Open in new tab
  • Download powerpoint
Figure 10.

Robustness of weak HOIs in simulated BOLD signals in various conditions under both the GSR and noGSR conditions. A, Across different ratios of inhibitory connections under the assumption of zero transmission delay. B, Similar to A but with transmission delays randomly chosen from 0 to 50 ms for all connections. C, Across different time bin widths, with peak-value thresholding method. TRs = 240 ms. In A–C, the FPN was used as the underlying structure of the network. Global coupling was set to 0.1. Common input was set to 0. D, Across different clustering coefficients of the underlying network. For each value of clustering coefficient, 6 11-node networks were generated with the mean degree of 4. Data are presented as mean ± SD. Global coupling was set to 0.087. Common input was set to 0. Inset, Performance measured by DJS of the Exp method and the DG model in the two-fold cross-validation test plotted as a function of the clustering coefficients. The threshold used here was the mean plus 1 SD.

In summary, our simulation results demonstrated that, given a similar level of overall pairwise interaction strength as the empirical data, the network model exhibited very weak true HOIs in a wide range of parameter settings. Combined with the results regarding the apparent HOIs in the model due to strong common input, they captured well the characteristic of weak HOIs under the GSR condition and spurious HOIs under the noGSR condition in the empirical data.

Discussion

In the present study, we quantified the strength of HOIs in BOLD signals recorded during the resting state. The results demonstrated that HOIs are generally weak in brain networks and thus play little role in forming functional systems at the macroscopic scale. Our modeling results further suggest that the lack of HOIs may be an intrinsic property of the macroscopic brain dynamics during the resting state.

Methodological considerations

Is it possible that the lack of HOIs is due to limitations in our analysis, for example, the dichotomization of BOLD signals, movement artifacts in the signal, or the slow fluctuation of the fMRI signals?

To reduce the dimensionality of the data for HOI analysis, we dichotomized continuous BOLD signals. Similar approaches have been used in analyzing local field potentials (LFPs), magnetoencephalography, and fMRI signals (Petermann et al., 2009; Yu et al., 2011; Watanabe et al., 2013; Shriki et al., 2013). It is known that fMRI signals are closely related to LFPs (Logothetis et al., 2001; Goense and Logothetis 2008; Magri et al., 2012; Logothetis and Panzeri 2015), whereas previous work has demonstrated that neuronal interactions can be correctly revealed by studying binarized LFP signals (Petermann et al., 2009; Yu et al., 2011, 2013). Our results showing the monotonic relationship between correlations measured based on the continuous BOLD signals and their binarized counterparts provide more direct evidence that the correlation structure in the original signals was indeed well preserved after the thresholding. Importantly, we emphasize that the threholding operation introduced a specific structure of HOIs, which can be well captured by the DG model. Therefore, the good performance of the DG model actually indicates the lack of true HOIs in the original, continuous BOLD signals.

It is known that head movements during fMRI scanning can induce systematic, yet spurious correlations among different brain areas (Power et al., 2012; Van Dijk et al., 2012; Power et al., 2014), thereby artificially increasing the synchrony of the recorded signals. In the current study, we found that there was no intrinsic HOI in either the GSR or noGSR condition. Because GSR is an effective way to remove motion artifacts in the fMRI data (Power et al., 2014), these results provide evidence that our conclusion are not affected by possible movement artifacts. This point was further supported by the fact that consistent results were obtained by using a more stringent method of motion artifact rejection (Ciric et al., 2017).

Last, is it possible that HOIs are difficult to be detected given the time scale of fluctuations in fMRI signals? For example, the interactions mediated by slow signals may need a longer temporal scale to reveal. However, we increased the width of time bin used for analyzing both the empirical and simulated BOLD signals in a wide range, but the results of weak HOIs were obtained consistently. Another possibility is that the HOIs existing in the faster temporal scales (e.g., neuronal interactions) were obscured in the slow fluctuations of the BOLD signals. This is also unlikely because we have analyzed the simulated synaptic activities, which have faster temporal fluctuation yet exhibit the same weak HOIs as the corresponding BLOD signals. Together, these results strongly suggest that the conclusion of weak HOIs may not depend on specific time scale of the analyzed BOLD signals.

Pairwise interactions determine activities of macroscopic brain networks

For microscopic brain networks, that is, the network of individual neurons, it is known that, although HOIs are in general weak (Schneidman et al., 2006; Tang et al., 2008; Yu et al., 2008), in certain conditions such as being driven by specific types of inputs (Ganmor et al., 2011; Macke et al., 2011) or with constituent neurons in close proximity (Ohiorhenuan et al., 2010), HOIs can be strong and no longer negligible. Previously, few studies have examined the HOIs at the macroscopic level. Although many researches have applied various graph-theory-based analyses on functional networks extracted using BOLD signals, these results do not provide information regarding how different brain areas can interact at higher orders because the widely used FC analysis is fundamentally based on measuring pairwise interactions. Therefore, examining the structure of HOIs will provide information about the functional organization of the brain complementary to the existing knowledge obtained by FC analyses.

Recently, Watanabe et al. (2013) reported that, compared with the independent model, the pairwise maximal entropy model (the Ising model) could significantly improve the accuracy in approximating activities of two functional networks under the special condition of GSR, although we demonstrated here that the Ising model falls short in predicting pattern probabilities for the condition without GSR. Using the DG model, we demonstrated that the lack of HOIs is a feature generally true for macroscopic brain networks with different sizes, topologies, degrees of spatial proximity, and spatial scales. These results provide a sound basis for the conventional, pairwise-based FC analysis, explaining why the FC analysis can indeed extract the functional organizations of the brain (Honey et al., 2007; Honey et al., 2009). In addition, we note that the lack of HOI does not conflict with the complex interaction structure, reflected by clusters/modules consisting of multiple brain regions, that was reported recently (Bassett et al., 2011; Stolz 2014; Davison et al., 2015). Such clusters could be formed by HOIs or aggregated pairwise interactions. Our results strongly imply that the latter probably is the case. From this perspective, the current results reveal the level of complexity at which the brain networks are organized, suggesting that it may be unnecessary to go beyond the pairwise level to study the macroscopic interactions in the resting brain.

Moreover, we took special attention to compare the conditions with and without the GSR. Although it was reported that the GSR can improve the effectiveness in detecting functional correlations, thereby providing better correspondence between observed FC and anatomy (Fox et al., 2009), other studies have suggested that the GS itself may have physiological importance (Schölvinck et al., 2010) and removing it will cause unnecessary distortions to the signals (Saad et al., 2012) or impede the detection of disease-related changes (Yang et al., 2014). Here, we did not aim to address the nature of the GS nor the validity in removing it. Instead, we demonstrated that, with or without the GS, the BOLD signals at rest do not contain true HOIs. Therefore, conventional FC analysis can be appropriately applied in both conditions.

It is worth noting that the current results demonstrated the validity of FC analysis at the pairwise level, but we do not address the question of which pairwise measures should be used. Although the simple PCC has been widely used, it has been suggested that more sophisticated pairwise measures such as partial correlation (Salvador et al., 2005; Salvador et al., 2008), magnitude squared coherence (Sun et al., 2004), mutual information (Salvador et al., 2008; Salvador et al., 2010), and exchange interactions (Watanabe et al., 2013) could provide more accurate results by taking into account the nonlinear correlations or the aggregated effects of pairwise correlations of the entire network.

Dynamic mechanism underlying the weak HOIs

It has been demonstrated that dynamic models are instrumental in revealing the relationship between the anatomical connectivity and dynamics of brain networks (Zhou et al., 2006; Honey et al., 2007; Honey et al., 2009). Therefore, we examined the HOIs in a network model constructed to match the structural and dynamical characteristics of the real brain networks. Our simulation results suggest a possible explanation for the lack of HOIs at the macroscopic brain, that is, at this scale the dynamics are driven by background inputs and recurrent interactions that are both weak. It is a regime close to a vicinity of linear fluctuation, which can be reasonably approximated by a Gaussian linear process: the DG model. Consistently, we found that the amplitude distributions of the empirical BOLD signals are Gaussian like. This was also the case for our simulated BOLD signals provided that the global coupling was <0.1 and the common input was <0.3, in which condition we found no strong HOIs can be generated. Therefore, the dynamic of macroscopic brain networks that is close to the linear fluctuation regime likely limits the existence of strong HOIs. This view is also supported by the finding that a feedforward neuronal circuit model generally produces very weak HOIs at the mesoscopic level unless bimodal common inputs are imposed to the circuit (Barreiro et al., 2014), indicating deviation from the vicinity of linear fluctuation.

We note that the current model has made a number of assumptions, including the effectiveness of the mean-field approximation, spatial homogeneity in both the mean-field model and the Balloon–Windkessel hemodynamic model, as well as specific setting of many parameters. These approximations and similar parameter setting were commonly used in previous studies (Friston et al., 2003; Wong and Wang, 2006; Deco et al., 2013) and we have verified the effectiveness of our model by showing the similarity of dynamical features between the model and the empirical data. Nevertheless, it awaits future studies to investigate to what degree these assumptions and specific parameters may affect the behavior of HOIs in a network.

Finally, we emphasize that the lack of strong HOIs demonstrated here is only applicable to the resting state so far. During execution of tasks, input to the systems may be changed significantly (Ganmor et al., 2011). In addition, interactions may undergo specific fine tuning to perform coordinated information processing. More importantly, functional benefits of HOIs such as XOR operation (Schneidman et al., 2006) may play an essential role in such conditions. These issues need to be further explored by studying HOIs under task conditions and, if they are present, determining how they can emerge during the transition from the resting state to the task execution state.

Footnotes

  • This work was supported by the National Key Research and Development Program of China (Grant 2017YFA0105203), the Natural Science Foundation of China (Grants 81471368, 11505283, 91132301, and 91432302), the Strategic Priority Research Program of the Chinese Academy of Sciences (CAS Grants XDB02060003 and XDB02030300), the China Postdoctoral Science Foundation (Special Financial Grant 2015T80154), the Hundred-Talent Program of CAS (S.Y.), the CAS Key Research Program of Frontier Sciences (Grant QYZDJ-SSW-SMC019), and the Beijing Municipal Science and Technology Commission (Grant Z161100000216139).

  • The authors declare no competing financial interests.

  • Correspondence should be addressed to Shan Yu, Brainnetome Center, Institute of Automation, Chinese Academy of Sciences, 95 Zhongguancun East Road, Beijing 100190, China. shan.yu{at}nlpr.ia.ac.cn

References

  1. ↵
    1. Abeles M
    (1991) Corticonics: neural circuits of the cerebral cortex. Cambridge: Cambridge University.
  2. ↵
    1. Amari S,
    2. Nakahara H,
    3. Wu S,
    4. Sakai Y
    (2003) Synchronous firing and higher-order interactions in neuron pool. Neural Comput 15:127–142. doi:10.1162/089976603321043720 pmid:12590822
    OpenUrlCrossRefPubMed
  3. ↵
    1. Barreiro AK,
    2. Gjorgjieva J,
    3. Rieke F,
    4. Shea-Brown E
    (2014) When do microcircuits produce beyond-pairwise correlations? Front Comput Neurosci 8:10. pmid:24567715
    OpenUrlPubMed
  4. ↵
    1. Bassett DS,
    2. Wymbs NF,
    3. Porter MA,
    4. Mucha PJ,
    5. Carlson JM,
    6. Grafton ST
    (2011) Dynamic reconfiguration of human brain networks during learning. Proc Natl Acad Sci U S A 108:7641–7646. doi:10.1073/pnas.1018985108 pmid:21502525
    OpenUrlAbstract/FREE Full Text
  5. ↵
    1. Beckmann CF,
    2. DeLuca M,
    3. Devlin JT,
    4. Smith SM
    (2005) Investigations into resting-state connectivity using independent component analysis. Philos Trans R Soc Lond B Biol Sci 360:1001–1013. pmid:16087444
    OpenUrlAbstract/FREE Full Text
  6. ↵
    1. Biswal B,
    2. Yetkin FZ,
    3. Haughton VM,
    4. Hyde JS
    (1995) Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn Reson Med 34:537–541. doi:10.1002/mrm.1910340409 pmid:8524021
    OpenUrlCrossRefPubMed
  7. ↵
    1. Braitenberg V,
    2. Schüz A
    (2013) Anatomy of the cortex: statistics and geometry. New York: Springer Science and Business Media.
  8. ↵
    1. Brier MR,
    2. Thomas JB,
    3. Snyder AZ,
    4. Benzinger TL,
    5. Zhang D,
    6. Raichle ME,
    7. Holtzman DM,
    8. Morris JC,
    9. Ances BM
    (2012) Loss of intranetwork and internetwork resting state functional connections with Alzheimer's disease progression. J Neurosci 32:8890–8899. doi:10.1523/JNEUROSCI.5698-11.2012 pmid:22745490
    OpenUrlAbstract/FREE Full Text
  9. ↵
    1. Brunel N,
    2. Wang XJ
    (2001) Effects of neuromodulation in a cortical network model of object working memory dominated by recurrent inhibition. J Comput Neurosci 11:63–85. doi:10.1023/A:1011204814320 pmid:11524578
    OpenUrlCrossRefPubMed
  10. ↵
    1. Buckner RL,
    2. Andrews-Hanna JR,
    3. Schacter DL
    (2008) The brain's default network. Ann N Y Acad Sci 1124:1–38. doi:10.1196/annals.1440.011 pmid:18400922
    OpenUrlCrossRefPubMed
  11. ↵
    1. Bullmore E,
    2. Sporns O
    (2009) Complex brain networks: graph theoretical analysis of structural and functional systems. Nat Rev Neurosci 10:186–198. doi:10.1038/nrn2575 pmid:19190637
    OpenUrlCrossRefPubMed
  12. ↵
    1. Ciric R,
    2. Wolf DH,
    3. Power JD,
    4. Roalf DR,
    5. Baum GL,
    6. Ruparel K,
    7. Shinohara RT,
    8. Elliott MA,
    9. Eickhoff SB,
    10. Davatzikos C,
    11. Gur RC,
    12. Gur RE,
    13. Bassett DS,
    14. Satterthwaite TD
    (2017) Benchmarking of participant-level confound regression strategies for the control of motion artifact in studies of functional connectivity. Neuroimage 154:174–187. doi:10.1016/j.neuroimage.2017.03.020 pmid:28302591
    OpenUrlCrossRefPubMed
  13. ↵
    1. Daffner KR,
    2. Scinto LF,
    3. Weitzman AM,
    4. Faust R,
    5. Rentz DM,
    6. Budson AE,
    7. Holcomb PJ
    (2003) Frontal and parietal components of a cerebral network mediating voluntary attention to novel events. J Cogn Neurosci 15:294–313. doi:10.1162/089892903321208213 pmid:12683359
    OpenUrlCrossRefPubMed
  14. ↵
    1. Davison EN,
    2. Schlesinger KJ,
    3. Bassett DS,
    4. Lynall ME,
    5. Miller MB,
    6. Grafton ST,
    7. Carlson JM
    (2015) Brain network adaptability across task states. PLoS Comput Biol 11:e1004029. doi:10.1371/journal.pcbi.1004029 pmid:25569227
    OpenUrlCrossRefPubMed
  15. ↵
    1. Deco G,
    2. Ponce-Alvarez A,
    3. Mantini D,
    4. Romani GL,
    5. Hagmann P,
    6. Corbetta M
    (2013) Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. J Neurosci 33:11239–11252. doi:10.1523/JNEUROSCI.1091-13.2013 pmid:23825427
    OpenUrlAbstract/FREE Full Text
  16. ↵
    1. Deco G,
    2. Ponce-Alvarez A,
    3. Hagmann P,
    4. Romani GL,
    5. Mantini D,
    6. Corbetta M
    (2014) How local excitation–inhibition ratio impacts the whole brain dynamics. J Neurosci 34:7886–7898. doi:10.1523/JNEUROSCI.5068-13.2014 pmid:24899711
    OpenUrlAbstract/FREE Full Text
  17. ↵
    1. Dosenbach NU,
    2. Fair DA,
    3. Cohen AL,
    4. Schlaggar BL,
    5. Petersen SE
    (2008) A dual-networks architecture of top-down control. Trends Cogn Sci 12:99–105. doi:10.1016/j.tics.2008.01.001 pmid:18262825
    OpenUrlCrossRefPubMed
  18. ↵
    1. Feinberg DA,
    2. Moeller S,
    3. Smith SM,
    4. Auerbach E,
    5. Ramanna S,
    6. Gunther M,
    7. Glasser MF,
    8. Miller KL,
    9. Ugurbil K,
    10. Yacoub E
    (2010) Multiplexed echo planar imaging for sub-second whole brain FMRI and fast diffusion imaging. PLoS One 5:e15710. doi:10.1371/journal.pone.0015710 pmid:21187930
    OpenUrlCrossRefPubMed
  19. ↵
    1. Fox MD,
    2. Raichle ME
    (2007) Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat Rev Neurosci 8:700–711. doi:10.1038/nrn2201 pmid:17704812
    OpenUrlCrossRefPubMed
  20. ↵
    1. Fox MD,
    2. Snyder AZ,
    3. Vincent JL,
    4. Corbetta M,
    5. Van Essen DC,
    6. Raichle ME
    (2005) The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc Natl Acad Sci U S A 102:9673–9678. doi:10.1073/pnas.0504136102 pmid:15976020
    OpenUrlAbstract/FREE Full Text
  21. ↵
    1. Fox MD,
    2. Corbetta M,
    3. Snyder AZ,
    4. Vincent JL,
    5. Raichle ME
    (2006) Spontaneous neuronal activity distinguishes human dorsal and ventral attention systems. Proc Natl Acad Sci U S A 103:10046–10051. doi:10.1073/pnas.0604187103 pmid:16788060
    OpenUrlAbstract/FREE Full Text
  22. ↵
    1. Fox MD,
    2. Zhang D,
    3. Snyder AZ,
    4. Raichle ME
    (2009) The global signal and observed anticorrelated resting state brain networks. J Neurophysiol 101:3270–3283. doi:10.1152/jn.90777.2008 pmid:19339462
    OpenUrlAbstract/FREE Full Text
  23. ↵
    1. Friston KJ,
    2. Harrison L,
    3. Penny W
    (2003) Dynamic causal modelling. Neuroimage 19:1273–1302. doi:10.1016/S1053-8119(03)00202-7 pmid:12948688
    OpenUrlCrossRefPubMed
  24. ↵
    1. Ganmor E,
    2. Segev R,
    3. Schneidman E
    (2011) Sparse low-order interaction network underlies a highly correlated and learnable neural population code. Proc Natl Acad Sci U S A 108:9679–9684. doi:10.1073/pnas.1019641108 pmid:21602497
    OpenUrlAbstract/FREE Full Text
  25. ↵
    1. Glasser MF,
    2. Sotiropoulos SN,
    3. Wilson JA,
    4. Coalson TS,
    5. Fischl B,
    6. Andersson JL,
    7. Xu J,
    8. Jbabdi S,
    9. Webster M,
    10. Polimeni JR,
    11. Van Essen DC,
    12. Jenkinson M,
    13. Jenkinson M
    (2013) The minimal preprocessing pipelines for the Human Connectome Project. Neuroimage 80:105–124. doi:10.1016/j.neuroimage.2013.04.127 pmid:23668970
    OpenUrlCrossRefPubMed
  26. ↵
    1. Goense JB,
    2. Logothetis NK
    (2008) Neurophysiology of the BOLD fMRI signal in awake monkeys. Curr Biol 18:631–640. doi:10.1016/j.cub.2008.03.054 pmid:18439825
    OpenUrlCrossRefPubMed
  27. ↵
    1. Gong G,
    2. He Y,
    3. Concha L,
    4. Lebel C,
    5. Gross DW,
    6. Evans AC,
    7. Beaulieu C
    (2009) Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. Cereb Cortex 19:524–536. doi:10.1093/cercor/bhn102 pmid:18567609
    OpenUrlCrossRefPubMed
  28. ↵
    1. Greicius MD,
    2. Krasnow B,
    3. Reiss AL,
    4. Menon V
    (2003) Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc Natl Acad Sci U S A 100:253–258. doi:10.1073/pnas.0135058100 pmid:12506194
    OpenUrlAbstract/FREE Full Text
  29. ↵
    1. Greicius MD,
    2. Srivastava G,
    3. Reiss AL,
    4. Menon V
    (2004) Default-mode network activity distinguishes Alzheimer's disease from healthy aging: evidence from functional MRI. Proc Natl Acad Sci U S A 101:4637–4642. doi:10.1073/pnas.0308627101 pmid:15070770
    OpenUrlAbstract/FREE Full Text
  30. ↵
    1. Greicius MD,
    2. Flores BH,
    3. Menon V,
    4. Glover GH,
    5. Solvason HB,
    6. Kenna H,
    7. Reiss AL,
    8. Schatzberg AF
    (2007) Resting-state functional connectivity in major depression: abnormally increased contributions from subgenual cingulate cortex and thalamus. Biol Psychiatry 62:429–437. doi:10.1016/j.biopsych.2006.09.020 pmid:17210143
    OpenUrlCrossRefPubMed
  31. ↵
    1. Honey CJ,
    2. Kötter R,
    3. Breakspear M,
    4. Sporns O
    (2007) Network structure of Cereb Cortex shapes functional connectivity on multiple time scales. Proc Natl Acad Sci U S A 104:10240–10245. doi:10.1073/pnas.0701519104 pmid:17548818
    OpenUrlAbstract/FREE Full Text
  32. ↵
    1. Honey CJ,
    2. Sporns O,
    3. Cammoun L,
    4. Gigandet X,
    5. Thiran JP,
    6. Meuli R,
    7. Hagmann P
    (2009) Predicting human resting-state functional connectivity from structural connectivity. Proc Natl Acad Sci U S A 106:2035–2040. doi:10.1073/pnas.0811168106 pmid:19188601
    OpenUrlAbstract/FREE Full Text
  33. ↵
    1. Jenkinson M,
    2. Bannister P,
    3. Brady M,
    4. Smith S
    (2002) Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage 17:825–841. doi:10.1006/nimg.2002.1132 pmid:12377157
    OpenUrlCrossRefPubMed
  34. ↵
    1. Logothetis NK,
    2. Panzeri S
    (2015) Local field potential, relationship to BOLD signal. In: Encyclopedia of computational neuroscience (Jaeger D and Jung R, ed), pp 1560–1568. New York: Springer.
  35. ↵
    1. Logothetis NK,
    2. Pauls J,
    3. Augath M,
    4. Trinath T,
    5. Oeltermann A
    (2001) Neurophysiological investigation of the basis of the fMRI signal. Nature 412:150–157. doi:10.1038/35084005 pmid:11449264
    OpenUrlCrossRefPubMed
  36. ↵
    1. Macke JH,
    2. Berens P,
    3. Ecker AS,
    4. Tolias AS,
    5. Bethge M
    (2009) Generating spike trains with specified correlation coefficients. Neural Comput 21:397–423. doi:10.1162/neco.2008.02-08-713 pmid:19196233
    OpenUrlCrossRefPubMed
  37. ↵
    1. Macke JH,
    2. Opper M,
    3. Bethge M
    (2011) Common input explains higher-order correlations and entropy in a simple model of neural population activity. Phys Rev Lett 106:208102. doi:10.1103/PhysRevLett.106.208102 pmid:21668265
    OpenUrlCrossRefPubMed
  38. ↵
    1. Magri C,
    2. Schridde U,
    3. Murayama Y,
    4. Panzeri S,
    5. Logothetis NK
    (2012) The amplitude and timing of the BOLD signal reflects the relationship between local field potential power at different frequencies. J Neurosci 32:1395–1407. doi:10.1523/JNEUROSCI.3985-11.2012 pmid:22279224
    OpenUrlAbstract/FREE Full Text
  39. ↵
    1. Moeller S,
    2. Yacoub E,
    3. Olman CA,
    4. Auerbach E,
    5. Strupp J,
    6. Harel N,
    7. Uğurbil K
    (2010) Multiband multislice GE-EPI at 7 tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI. Magn Reson Med 63:1144–1153. doi:10.1002/mrm.22361 pmid:20432285
    OpenUrlCrossRefPubMed
  40. ↵
    1. Murphy K,
    2. Birn RM,
    3. Handwerker DA,
    4. Jones TB,
    5. Bandettini PA
    (2009) The impact of global signal regression on resting state correlations: are anti-correlated networks introduced? Neuroimage 44:893–905. doi:10.1016/j.neuroimage.2008.09.036 pmid:18976716
    OpenUrlCrossRefPubMed
  41. ↵
    1. Nakahara H,
    2. Amari S
    (2002) Information-geometric measure for neural spikes. Neural Comput 14:2269–2316. doi:10.1162/08997660260293238 pmid:12396564
    OpenUrlCrossRefPubMed
  42. ↵
    1. Nichols TE,
    2. Holmes AP
    (2002) Nonparametric permutation tests for functional neuroimaging: a primer with examples. Hum Brain Mapp 15:1–25. doi:10.1002/hbm.1058 pmid:11747097
    OpenUrlCrossRefPubMed
  43. ↵
    1. Ohiorhenuan IE,
    2. Victor JD
    (2011) Information-geometric measure of 3-neuron firing patterns characterizes scale dependence in cortical networks. J Comput Neurosci 30:125–141. doi:10.1007/s10827-010-0257-0 pmid:20635129
    OpenUrlCrossRefPubMed
  44. ↵
    1. Ohiorhenuan IE,
    2. Mechler F,
    3. Purpura KP,
    4. Schmid AM,
    5. Hu Q,
    6. Victor JD
    (2010) Sparse coding and high-order correlations in fine-scale cortical networks. Nature 466:617–621. doi:10.1038/nature09178 pmid:20601940
    OpenUrlCrossRefPubMed
  45. ↵
    1. Petermann T,
    2. Thiagarajan TC,
    3. Lebedev MA,
    4. Nicolelis MA,
    5. Chialvo DR,
    6. Plenz D
    (2009) Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc Natl Acad Sci U S A 106:15921–15926. doi:10.1073/pnas.0904089106 pmid:19717463
    OpenUrlAbstract/FREE Full Text
  46. ↵
    1. Power JD,
    2. Barnes KA,
    3. Snyder AZ,
    4. Schlaggar BL,
    5. Petersen SE
    (2012) Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. Neuroimage 59:2142–2154. doi:10.1016/j.neuroimage.2011.10.018 pmid:22019881
    OpenUrlCrossRefPubMed
  47. ↵
    1. Power JD,
    2. Mitra A,
    3. Laumann TO,
    4. Snyder AZ,
    5. Schlaggar BL,
    6. Petersen SE
    (2014) Methods to detect, characterize, and remove motion artifact in resting state fMRI. Neuroimage 84:320–341. doi:10.1016/j.neuroimage.2013.08.048 pmid:23994314
    OpenUrlCrossRefPubMed
  48. ↵
    1. Raichle ME,
    2. MacLeod AM,
    3. Snyder AZ,
    4. Powers WJ,
    5. Gusnard DA,
    6. Shulman GL
    (2001) A default mode of brain function. Proc Natl Acad Sci U S A 98:676–682. doi:10.1073/pnas.98.2.676 pmid:11209064
    OpenUrlAbstract/FREE Full Text
  49. ↵
    1. Roudi Y,
    2. Nirenberg S,
    3. Latham PE
    (2009) Pairwise maximum entropy models for studying large biological systems: when they can work and when they can't. PLOS Comput Biol 5:e1000380. doi:10.1371/journal.pcbi.1000380 pmid:19424487
    OpenUrlCrossRefPubMed
  50. ↵
    1. Saad ZS,
    2. Gotts SJ,
    3. Murphy K,
    4. Chen G,
    5. Jo HJ,
    6. Martin A,
    7. Cox RW
    (2012) Trouble at rest: how correlation patterns and group differences become distorted after global signal regression. Brain Connect 2:25–32. doi:10.1089/brain.2012.0080 pmid:22432927
    OpenUrlCrossRefPubMed
  51. ↵
    1. Salvador R,
    2. Suckling J,
    3. Coleman MR,
    4. Pickard JD,
    5. Menon D,
    6. Bullmore E
    (2005) Neurophysiological architecture of functional magnetic resonance images of human brain. Cereb Cortex 15:1332–1342. pmid:15635061
    OpenUrlCrossRefPubMed
  52. ↵
    1. Salvador R,
    2. Martínez A,
    3. Pomarol-Clotet E,
    4. Gomar J,
    5. Vila F,
    6. **Sarró S,
    7. Capdevila A,
    8. Bullmore E
    (2008) A simple view of the brain through a frequency-specific functional connectivity measure. Neuroimage 39:279–289. doi:10.1016/j.neuroimage.2007.08.018 pmid:17919927
    OpenUrlCrossRefPubMed
  53. ↵
    1. Salvador R,
    2. Anguera M,
    3. Gomar JJ,
    4. Bullmore ET,
    5. Pomarol-Clotet E
    (2010) Conditional mutual information maps as descriptors of net connectivity levels in the brain. Front Neuroinform 4:115. doi:10.3389/fninf.2010.00115 pmid:21151357
    OpenUrlCrossRefPubMed
  54. ↵
    1. Sara SJ
    (2009) The locus coeruleus and noradrenergic modulation of cognition. Nat Rev Neurosci 10:211–213. doi:10.1038/nrn2573 pmid:19190638
    OpenUrlCrossRefPubMed
  55. ↵
    1. Schneidman E,
    2. Still S,
    3. Berry MJ 2nd.,
    4. Bialek W
    (2003) Network information and connected correlations. Phys Rev Lett 91:238701. doi:10.1103/PhysRevLett.91.238701 pmid:14683220
    OpenUrlCrossRefPubMed
  56. ↵
    1. Schneidman E,
    2. Berry MJ 2nd.,
    3. Segev R,
    4. Bialek W
    (2006) Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440:1007–1012. doi:10.1038/nature04701 pmid:16625187
    OpenUrlCrossRefPubMed
  57. ↵
    1. Schölvinck ML,
    2. Maier A,
    3. Ye FQ,
    4. Duyn JH,
    5. Leopold DA
    (2010) Neural basis of global resting-state fMRI activity. Proc Natl Acad Sci U S A 107:10238–10243. doi:10.1073/pnas.0913110107 pmid:20439733
    OpenUrlAbstract/FREE Full Text
  58. ↵
    1. Seeley WW,
    2. Menon V,
    3. Schatzberg AF,
    4. Keller J,
    5. Glover GH,
    6. Kenna H,
    7. Reiss AL,
    8. Greicius MD
    (2007) Dissociable intrinsic connectivity networks for salience processing and executive control. J Neurosci 27:2349–2356. doi:10.1523/JNEUROSCI.5587-06.2007 pmid:17329432
    OpenUrlAbstract/FREE Full Text
  59. ↵
    1. Setsompop K,
    2. Gagoski BA,
    3. Polimeni JR,
    4. Witzel T,
    5. Wedeen VJ,
    6. Wald LL
    (2012) Blipped-controlled aliasing in parallel imaging for simultaneous multislice echo planar imaging with reduced g-factor penalty. Magn Reson Med 67:1210–1224. doi:10.1002/mrm.23097 pmid:21858868
    OpenUrlCrossRefPubMed
  60. ↵
    1. Shemesh Y,
    2. Sztainberg Y,
    3. Forkosh O,
    4. Shlapobersky T,
    5. Chen A,
    6. Schneidman E
    (2013) High-order social interactions in groups of mice. Elife 2:e00759. doi:10.7554/eLife.00759 pmid:24015357
    OpenUrlAbstract/FREE Full Text
  61. ↵
    1. Shriki O,
    2. Alstott J,
    3. Carver F,
    4. Holroyd T,
    5. Henson RN,
    6. Smith ML,
    7. Coppola R,
    8. Bullmore E,
    9. Plenz D
    (2013) Neuronal avalanches in the resting MEG of the human brain. J Neurosci 33:7079–7090. doi:10.1523/JNEUROSCI.4286-12.2013 pmid:23595765
    OpenUrlAbstract/FREE Full Text
  62. ↵
    1. Smith SM,
    2. Jenkinson M,
    3. Woolrich MW,
    4. Beckmann CF,
    5. Behrens TE,
    6. Johansen-Berg H,
    7. Bannister PR,
    8. De Luca M,
    9. Drobnjak I,
    10. Flitney DE,
    11. Niazy RK,
    12. Saunders J,
    13. Vickers J,
    14. Zhang Y,
    15. De Stefano N,
    16. Brady JM,
    17. Matthews PM
    (2004) Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage 23:S208–S219. doi:10.1016/j.neuroimage.2004.07.051 pmid:15501092
    OpenUrlCrossRefPubMed
  63. ↵
    1. Stolz B
    (2014) Computational topology in neuroscience. Master's thesis, University of Oxford, Oxford, United Kingdom.
  64. ↵
    1. Sun FT,
    2. Miller LM,
    3. D'Esposito M
    (2004) Measuring interregional functional connectivity using coherence and partial coherence analyses of fMRI data. Neuroimage 21:647–658. doi:10.1016/j.neuroimage.2003.09.056 pmid:14980567
    OpenUrlCrossRefPubMed
  65. ↵
    1. Tang A,
    2. Jackson D,
    3. Hobbs J,
    4. Chen W,
    5. Smith JL,
    6. Patel H,
    7. Prieto A,
    8. Petrusca D,
    9. Grivich MI,
    10. Sher A,
    11. Hottowy P,
    12. Dabrowski W,
    13. Litke AM,
    14. Beggs JM
    (2008) A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro. J Neurosci 28:505–518. doi:10.1523/JNEUROSCI.3359-07.2008 pmid:18184793
    OpenUrlAbstract/FREE Full Text
  66. ↵
    1. Van Dijk KR,
    2. Sabuncu MR,
    3. Buckner RL
    (2012) The influence of head motion on intrinsic functional connectivity MRI. Neuroimage 59:431–438. doi:10.1016/j.neuroimage.2011.07.044 pmid:21810475
    OpenUrlCrossRefPubMed
  67. ↵
    1. Van Essen DC, et al
    . (2012) The Human Connectome Project: a data acquisition perspective. Neuroimage 62:2222–2231. doi:10.1016/j.neuroimage.2012.02.018 pmid:22366334
    OpenUrlCrossRefPubMed
  68. ↵
    1. Wang XJ
    (2002) Probabilistic decision making by slow reverberation in cortical circuits. Neuron 36:955–968. doi:10.1016/S0896-6273(02)01092-9 pmid:12467598
    OpenUrlCrossRefPubMed
  69. ↵
    1. Watanabe T,
    2. Hirose S,
    3. Wada H,
    4. Imai Y,
    5. Machida T,
    6. Shirouzu I,
    7. Konishi S,
    8. Miyashita Y,
    9. Masuda N
    (2013) A pairwise maximum entropy model accurately describes resting-state human brain networks. Nat Commun 4:1370. doi:10.1038/ncomms2388 pmid:23340410
    OpenUrlCrossRefPubMed
  70. ↵
    1. Watts DJ,
    2. Strogatz SH
    (1998) Collective dynamics of ‘small-world’ networks. Nature 393:440–442. doi:10.1038/30918 pmid:9623998
    OpenUrlCrossRefPubMed
  71. ↵
    1. Wong KF,
    2. Wang XJ
    (2006) A recurrent network mechanism of time integration in perceptual decisions. J Neurosci 26:1314–1328. doi:10.1523/JNEUROSCI.3733-05.2006 pmid:16436619
    OpenUrlAbstract/FREE Full Text
  72. ↵
    1. Yang GJ,
    2. Murray JD,
    3. Repovs G,
    4. Cole MW,
    5. Savic A,
    6. Glasser MF,
    7. Pittenger C,
    8. Krystal JH,
    9. Wang XJ,
    10. Pearlson GD,
    11. Glahn DC,
    12. Anticevic A
    (2014) Altered global brain signal in schizophrenia. Proc Natl Acad Sci U S A 111:7438–7443. doi:10.1073/pnas.1405289111 pmid:24799682
    OpenUrlAbstract/FREE Full Text
  73. ↵
    1. Yu S,
    2. Huang D,
    3. Singer W,
    4. Nikolic D
    (2008) A small world of neuronal synchrony. Cereb Cortex 18:2891–2901. doi:10.1093/cercor/bhn047 pmid:18400792
    OpenUrlCrossRefPubMed
  74. ↵
    1. Yu S,
    2. Yang H,
    3. Nakahara H,
    4. Santos GS,
    5. Nikolić D,
    6. Plenz D
    (2011) Higher-order interactions characterized in cortical activity. J Neurosci 31:17514–17526. doi:10.1523/JNEUROSCI.3127-11.2011 pmid:22131413
    OpenUrlAbstract/FREE Full Text
  75. ↵
    1. Yu S,
    2. Yang H,
    3. Shriki O,
    4. Plenz D
    (2013) Universal organization of resting brain activity at the thermodynamic critical point. Front Syst Neurosci 7:42. doi:10.3389/fnsys.2013.00042 pmid:23986660
    OpenUrlCrossRefPubMed
  76. ↵
    1. Zanto TP,
    2. Gazzaley A
    (2013) Fronto-parietal network: flexible hub of cognitive control. Trends Cogn Sci 17:602–603. doi:10.1016/j.tics.2013.10.001 pmid:24129332
    OpenUrlCrossRefPubMed
  77. ↵
    1. Zhan Y,
    2. Ma J,
    3. Alexander-Bloch A,
    4. Xu K,
    5. Cui Y,
    6. Feng Q,
    7. Jiang T,
    8. Liu Y
    (2016) Longitudinal study of impaired intra- and internetwork brain connect in subjects at high risk for Alzheimer's disease. J Alzheimers Dis 52:913–927. doi:10.3233/JAD-160008 pmid:27060962
    OpenUrlCrossRefPubMed
  78. ↵
    1. Zhou C,
    2. Zemanová L,
    3. Zamora G,
    4. Hilgetag CC,
    5. Kurths J
    (2006) Hierarchical organization unveiled by functional connectivity in complex brain networks. Phys Rev Lett 97:238103. doi:10.1103/PhysRevLett.97.238103 pmid:17280251
    OpenUrlCrossRefPubMed
  79. ↵
    1. Zhou J,
    2. Greicius MD,
    3. Gennatas ED,
    4. Growdon ME,
    5. Jang J,
    6. Rabinovici GD,
    7. Kramer JH,
    8. Weiner M,
    9. Miller BL,
    10. Seeley WW
    (2010) Divergent network connectivity changes in behavioural variant frontotemporal dementia and Alzheimer's disease. Brain 133:1352–1367. doi:10.1093/brain/awq075 pmid:20410145
    OpenUrlCrossRefPubMed
  80. ↵
    1. Zielinski BA,
    2. Gennatas ED,
    3. Zhou J,
    4. Seeley WW
    (2010) Network-level structural covariance in the developing brain. Proc Natl Acad Sci U S A 107:18191–18196. doi:10.1073/pnas.1003109107 pmid:20921389
    OpenUrlAbstract/FREE Full Text
View Abstract
Back to top

In this issue

The Journal of Neuroscience: 37 (43)
Journal of Neuroscience
Vol. 37, Issue 43
25 Oct 2017
  • Table of Contents
  • Table of Contents (PDF)
  • About the Cover
  • Index by author
  • Advertising (PDF)
  • Ed Board (PDF)
Email

Thank you for sharing this Journal of Neuroscience article.

NOTE: We request your email address only to inform the recipient that it was you who recommended this article, and that it is not junk mail. We do not retain these email addresses.

Enter multiple addresses on separate lines or separate them with commas.
Weak Higher-Order Interactions in Macroscopic Functional Networks of the Resting Brain
(Your Name) has forwarded a page to you from Journal of Neuroscience
(Your Name) thought you would be interested in this article in Journal of Neuroscience.
CAPTCHA
This question is for testing whether or not you are a human visitor and to prevent automated spam submissions.
Print
View Full Page PDF
Citation Tools
Weak Higher-Order Interactions in Macroscopic Functional Networks of the Resting Brain
Xuhui Huang, Kaibin Xu, Congying Chu, Tianzi Jiang, Shan Yu
Journal of Neuroscience 25 October 2017, 37 (43) 10481-10497; DOI: 10.1523/JNEUROSCI.0451-17.2017

Citation Manager Formats

  • BibTeX
  • Bookends
  • EasyBib
  • EndNote (tagged)
  • EndNote 8 (xml)
  • Medlars
  • Mendeley
  • Papers
  • RefWorks Tagged
  • Ref Manager
  • RIS
  • Zotero
Respond to this article
Request Permissions
Share
Weak Higher-Order Interactions in Macroscopic Functional Networks of the Resting Brain
Xuhui Huang, Kaibin Xu, Congying Chu, Tianzi Jiang, Shan Yu
Journal of Neuroscience 25 October 2017, 37 (43) 10481-10497; DOI: 10.1523/JNEUROSCI.0451-17.2017
del.icio.us logo Digg logo Reddit logo Twitter logo CiteULike logo Facebook logo Google logo Mendeley logo
  • Tweet Widget
  • Facebook Like
  • Google Plus One

Jump to section

  • Article
    • Abstract
    • Introduction
    • Materials and Methods
    • Results
    • Discussion
    • Footnotes
    • References
  • Figures & Data
  • Info & Metrics
  • eLetters
  • PDF

Keywords

  • default mode network
  • frontoparietal network
  • functional connectivity
  • pairwise correlation
  • resting-state fMRI

Responses to this article

Respond to this article

Jump to comment:

No eLetters have been published for this article.

Related Articles

Cited By...

More in this TOC Section

Research Articles

  • Prelimbic Cortex Activity during a Distress Tolerance Task Predicts Cocaine-Seeking Behavior in Male, But Not Female Rats
  • 2-AG-Mediated Control of GABAergic Signaling Is Impaired in a Model of Epilepsy
  • Identification of the Acid-Sensitive Site Critical for Chloral Hydrate (CH) Activation of the Proton-Activated Chloride Channel
Show more Research Articles

Systems/Circuits

  • 2-AG-Mediated Control of GABAergic Signaling Is Impaired in a Model of Epilepsy
  • Expression of a Form of Cerebellar Motor Memory Requires Learned Alterations to the Activity of Inhibitory Molecular Layer Interneurons
  • Quantitative Fluorescence Analysis Reveals Dendrite-Specific Thalamocortical Plasticity in L5 Pyramidal Neurons during Learning
Show more Systems/Circuits
  • Home
  • Alerts
  • Visit Society for Neuroscience on Facebook
  • Follow Society for Neuroscience on Twitter
  • Follow Society for Neuroscience on LinkedIn
  • Visit Society for Neuroscience on Youtube
  • Follow our RSS feeds

Content

  • Early Release
  • Current Issue
  • Issue Archive
  • Collections

Information

  • For Authors
  • For Advertisers
  • For the Media
  • For Subscribers

About

  • About the Journal
  • Editorial Board
  • Privacy Policy
  • Contact
(JNeurosci logo)
(SfN logo)

Copyright © 2023 by the Society for Neuroscience.
JNeurosci Online ISSN: 1529-2401

The ideas and opinions expressed in JNeurosci do not necessarily reflect those of SfN or the JNeurosci Editorial Board. Publication of an advertisement or other product mention in JNeurosci should not be construed as an endorsement of the manufacturer’s claims. SfN does not assume any responsibility for any injury and/or damage to persons or property arising from or related to any use of any material contained in JNeurosci.