Exponential Random Graph Modeling for Complex Brain Networks

Abstract

Exponential random graph models (ERGMs), also known as p* models, have been utilized extensively in the social science literature to study complex networks and how their global structure depends on underlying structural components. However, the literature on their use in biological networks (especially brain networks) has remained sparse. Descriptive models based on a specific feature of the graph (clustering coefficient, degree distribution, etc.) have dominated connectivity research in neuroscience. Corresponding generative models have been developed to reproduce one of these features. However, the complexity inherent in whole-brain network data necessitates the development and use of tools that allow the systematic exploration of several features simultaneously and how they interact to form the global network architecture. ERGMs provide a statistically principled approach to the assessment of how a set of interacting local brain network features gives rise to the global structure. We illustrate the utility of ERGMs for modeling, analyzing, and simulating complex whole-brain networks with network data from normal subjects. We also provide a foundation for the selection of important local features through the implementation and assessment of three selection approaches: a traditional p-value based backward selection approach, an information criterion approach (AIC), and a graphical goodness of fit (GOF) approach. The graphical GOF approach serves as the best method given the scientific interest in being able to capture and reproduce the structure of fitted brain networks.

About the Authors:

Sean L. Simpson

* E-mail: [email protected]

Affiliation: Department of Biostatistical Sciences, Wake Forest University School of Medicine, Winston-Salem, North Carolina, United States of America

Satoru Hayasaka

Affiliations Department of Biostatistical Sciences, Wake Forest University School of Medicine, Winston-Salem, North Carolina, United States of America, Department of Radiology, Wake Forest University School of Medicine, Winston-Salem, North Carolina, United States of America

Paul J. Laurienti

Affiliation: Department of Radiology, Wake Forest University School of Medicine, Winston-Salem, North Carolina, United States of America

Introduction

Brain networks

Whole-brain connectivity analyses are gaining prominence in the neuroscientific literature due to the need to understand how various regions of the brain interact with one another. The inherent complexity in the way these regions interact necessitates studying the brain as a whole rather than just its individual parts. The application of network and graph theory to the brain has facilitated these whole-brain analyses and helped to uncover new insights into the structure and function of the nervous system. Structural and functional connectivity studies have revealed that the brain exhibits the small-world properties [1]–[4]. These properties are characterized by tight local clustering and efficient long distance connections as described in the seminal work of [5]. Network models based on a given small-world property or other local property (e.g., node degree ()) have mostly been utilized as a means to describe various brain networks. However, in order to gain deeper insights into the complex neurobiological interactions and changes that occur in many neurological conditions and disorders, analysis methods that enable systematically assessing several properties simultaneously are needed given the statistical dependencies among these properties [6], [7]. The exponential random graph models discussed in this paper provide one such analysis approach.

Exponential random graph models

Exponential random graph models (ERGMs), also known as p* models [8]–[11], have been utilized extensively in the social science literature to analyze complex network data as discussed in [12], [13] and others. However, the literature on their use in biological networks (especially brain networks) has remained sparse. Descriptive models based on a specific feature of the network such as characteristic path length () and clustering coefficient () have dominated connectivity research in neuroscience [14]. The few inferential studies have employed relatively rudimentary testing techniques such as the ANOVA used in [15] to examine group differences based on one of these features. ERGMs provide a statistically principled approach to the systematic exploration of several features simultaneously and how they interact to form the global network architecture. They allow parsimoniously modeling the probability mass function (pmf) for a given class of graphs based on a set of explanatory metrics (local features). The pmf can then be used to determine the probability that any given graph is drawn from the same distribution as the observed graph. These models enable achieving an efficient representation of complex network data structures and allow examining the way in which a network's global structure and function depend on its local structure. That is, they provide a means of assessing how and to what extent combinations of local (brain) structures produce global network properties.

In ERGMs networks are analogous to a multivariate response variable in regression analysis, with the explanatory metrics quantifying local features of the network such as how clustered connections are (short distance communication) or how well the network transmits information globally (long distance communication). Fitted parameter values from the model can then be utilized to understand particular emergent behaviors of the network (how local features give rise to the global structure). These values can also be used to simulate random realizations of networks that retain constitutive characteristics of the original network.

A more intuitive way to view ERGMs in the brain network context are as models that quantify the relative significance of various graph/network measures (, , , etc.), or their analogues, in explaining the overall network structure, thus enabling generative conjectures about global architecture. These models provide several benefits for brain network researchers. They allow asking specific questions about processes that may give rise to the network architecture via the inclusion of explanatory metrics of choice. ERGMs inherently account for any confounding bias, like the (, )-dependence of network measures (where is the number of nodes and the average degree) detailed in [6], when the potential confounding variables are included in the model. The stochastic nature of the model allows understanding and quantifying the uncertainty (an intrinsic feature of complex biological processes) associated with our observed brain network(s) [12]. Simulations based on ERGM fits to brain networks (sets of selected network measures and their parameter estimates) can provide insight into biological variability via the distribution of possible brain networks produced. However, currently, the computational intensiveness of fitting ERGMs may preclude their use with very large networks (e.g., voxel-based networks with tens of thousands of nodes) and certain combinations of network measures.

Here we illustrate the utility of ERGMs for modeling, analyzing, and simulating complex whole-brain network. We also provide a foundation for the development of a “best assessment” ERGM for analyzing complex brain networks. Appropriate statistical comparisons between networks (or groups of networks) via ERGMs necessitates establishing one model (set of explanatory metrics/local features) in order to extract comparable parameter estimates due to the dependence of these features on each other. Toward this end, we assess three potential methods of feature selection for ERGMs in the brain network context. These approaches include a traditional p-value based backward selection approach, an information criterion approach (AIC), and a graphical goodness of fit (GOF) approach. Although the latter two techniques have been discussed in the context of ERGMS [16], [17], no detailed comparisons have been performed to determine whether the approaches generally produce the same “best” model/set of features.

Materials and Methods

Ethics statement

This study included 10 volunteers representing a subset of a previous study [18]. The study protocol, including all analyses performed here, was approved by the Wake Forest University School of Medicine Institutional Review Board. All subjects gave written informed consent in accordance with the Declaration of Helsinki.

Data and network construction

Our data include whole-brain functional connectivity networks for 10 normal subjects aged 20–35 (5 female, average age 27.7 years old [4.7 SD]). Each network is comprised of 90 nodes corresponding to the 90 brain regions (90 ROIs-Regions of Interest) defined by the Automated Anatomical Labeling atlas (AAL; [19]). The whole-brain networks were constructed based on fMRI images using graph theory methods. For each subject, 120 images were acquired during 5 minutes of resting using a gradient echo echoplanar imaging (EPI) protocol with TR/TE = 2500/40 ms on a 1.5 T GE twin-speed LX scanner with a birdcage head coil (GE Medical Systems, Milwaukee, WI). The acquired images were motion corrected, spatially normalized to the MNI (Montreal Neurological Institute) space and re-sliced to 4 mm voxel size using an in-house processing script based on SPM99 package (Wellcome Trust Centre for Neuroimaging, London, UK). The resulting images were not smoothed in order to avoid artificially introducing local spatial correlation [20].

The first step in performing the network construction was to generate a whole brain connectivity matrix, or adjacency matrix . This is a binary matrix where is the number of nodes representing 90 ROIs. The matrix notes the presence or absence of a connection between any two nodes ( and ). The determination of a connection between and was done by calculating a partial correlation coefficient adjusted for motion and physiological noises (see [21] for further details).

An unweighted, undirected network was then generated for each subject by applying a threshold to the correlation matrix to yield an adjacency matrix . In order to compare data across people, it is necessary to generate comparable networks. The network was defined so that the relationship between the number of nodes and the average node degree is the same across different subjects. In particular, the network was defined so that  = ()() is the same across subjects, with . This relationship is based on the path length of a random network with nodes and average degree [4], [5], and can be re-written as . Our analysis includes ERGM fits to these thresholded whole-brain functional connectivity networks for each subject (an example of which is shown in Figure 1).

[Figure omitted. See PDF.]

Figure 1. Network of subject 10 in brain space.

https://doi.org/10.1371/journal.pone.0020039.g001

Model definition

Exponential random graph models have the following form [22]:(1)Here is an ( nodes) random symmetric adjacency matrix representing a brain network from a particular class of networks, with if an edge exists between nodes and and otherwise. Nodes represent locations in the brain (e.g., ROIs) and edges represent functional or structural connections between them. We statistically model the probability mass function (pmf) of this class of networks as a function of the prespecified network features defined by the -dimensional vector g. This vector of explanatory metrics consists of covariates that are functions of the network and can contain any graph statistic (e.g., number of paths of length two) or node statistic (e.g., brain location of the node). The parameter vector , associated with g, quantifies the relative significance of the network features in explaining the structure of the network after accounting for the contribution of all other network features in the model and must be estimated. More specifically, indicates the change in the log odds of an edge existing for each unit increase in the corresponding explanatory metric. If the value corresponding to a given metric is large and positive, then that metric plays a considerable role in explaining the network architecture and is more prevalent than in the null model (random network with the probability of an edge existing ()). Conversely, if the value is large and negative, then that metric still plays a considerable role in explaining the network architecture but is less prevalent than in the null model. Consequently, inferences can be made about whether certain local features/substructures are observed in the network more than would be expected by chance enabling hypothesis development regarding the biological processes that produce these structural properties. The normalizing constant ensures that the probabilities sum to one. This approach allows representing the global network structure by locally specified explanatory metrics, thus providing a means to examine the nature of networks that are likely to emerge from these effects.

The goal in defining g is to identify local metrics that concisely summarize the global (whole-brain) network structure. Table 1 defines a subset of mathematically compatible explanatory network metrics (for further details see [16], [23], [24]). Several analogs to these metrics for directed graphs have been detailed by [25]. The GWD, GWESP, and GWDSP statistics discussed in [17] help address degeneracy issues illuminated in [22] and [26]. These issues concern the shape of the estimated pmf (e.g., a pmf in which only a few graphs have nonzero probability) and can lead to lack of model convergence and unreliable results. As noted by [16], [27], the most appropriate explanatory metrics vary by network type. Thus, an exploration of which network metrics best characterize brain networks has great appeal. Once the most appropriate statistics have been established, parameter profiles can be utilized to classify and compare whole-brain networks. These parameter profile comparisons require the use of a uniform set of explanatory metrics for all networks (due to metric interdependencies) and balanced networks (same number of nodes for all networks) due to the dependence of the metrics on network size.

[Figure omitted. See PDF.]

Table 1. Subset of explanatory network metrics.

https://doi.org/10.1371/journal.pone.0020039.t001

It is important to note that ERGMs can be thought of as a way of parameterizing models for networks, and are not a “kind” of network model in the way “model” is traditionally used in the brain network literature. Most other network models, in theory, should have an equivalent ERGM expression (though that specific expression may not be convenient, parsimonious, etc.). For instance, an ERGM with just the Edges metric (Table 1) in the formulation is equivalent to the Erdos-Renyi model. Thus, ERGMs allow parameterizations that subsume most (if not all) other network models.

Fitting of the ERGM in equation 1 is normally done with either Markov chain Monte Carlo maximum likelihood estimation (MCMC MLE) or maximum pseudo-likelihood estimation (MPLE) ([28] contains details). Model fits with MPLE are much simpler computationally than MCMC MLE fits and afford higher convergence rates with large networks. However, properties of the MPLE estimators are not well understood, and the estimates tend to be less accurate than those of MCMC MLE. Here we employ MCMC MLE to fit the model in equation 1 given that there were no convergence issues. See [29] for further details about this estimation approach which can be implemented in the statnet package [23] for the R statistical computing environment.

Model selection

In order to establish the most appropriate set of explanatory metrics for each subject's brain network and provide a foundation for the development of a “best assessment” ERGM for analyzing complex brain networks, we implemented and assessed three model/metric selection methods. They include a traditional p-value based backward selection approach [30], an information criterion approach (AIC, [31]), and a graphical goodness of fit (GOF) approach [17]. The latter two techniques are used most often for metric selection in ERGMs [16], [17]; and, to our knowledge, no detailed comparisons have been performed to determine whether the approaches generally produce the same “best” model. The p-value approach is based on removing metrics that are not statistically significant. Whereas, the AIC approach selects the set of metrics that produce the estimated distribution most likely to have resulted in the observed data with a penalty for additional metrics to ensure parsimony. Alternatively, the graphical GOF method allows subjectively selecting the set of explanatory metrics that produces the model most able to capture and reproduce certain topological properties of the observed network (see Appendix S1 for more details). For each approach ERGMs were fitted to the 90-node unweighted, undirected brain networks of the 10 subjects discussed previously. The potential explanatory metrics for each of the 10 networks are listed by category in Table 2. The categories were chosen based on properties of brain networks that are regarded as important in the literature [14]. These metrics are analogous to typical brain network metrics (e.g., clustering coefficient ()) but have been developed to be statistically compatible with ERGMs. Figure 2 illustrates the calculation of the less widely used of these statistics, namely GWESP, GWNSP, and GWDSP, on a six-node example network. The distribution of the unweighted analogues of these metrics (ESP, NSP, and DSP) is given for simplicity. The weighted versions simply sum the values of the distribution giving less weight to those with more shared partners. For this example we note that the network has 1 set of connected nodes with 1 shared partners (ESP), 5 sets with 1 shared partner (ESP), 1 set with 2 shared partners (ESP), and 0 sets with 3 or 4 shared partners (ESP and ESP). Further details on the metrics are provided in Table 1 and [27]. The parameters associated with GWESP, GWDSP, GWNSP, and GWD were all assumed to be fixed and known (for reasons outlined in [17]) and set to based on preliminary analyses as this value generally led to better fitting models according to all selection methods. The three aforementioned model selection approaches are outlined in Appendix S1.

[Figure omitted. See PDF.]

Figure 2. Six-node example network.

The edgewise, nonedgewise, and dyadwise shared partner distributions are (ESP, , ESP)(1, 5, 1, 0, 0), (NSP, , NSP)(1, 4, 3, 0, 0), and (DSP, , DSP)(2, 9, 4, 0, 0) respectively.

https://doi.org/10.1371/journal.pone.0020039.g002

[Figure omitted. See PDF.]

Table 2. Explanatory network metrics by category.

https://doi.org/10.1371/journal.pone.0020039.t002

Results

We implemented the model selection procedures delineated in the previous section and Appendix S1 for each of the 10 subjects using the statnet package [23] for the R statistical computing environment. The resulting models (for each approach) and their corresponding parameter estimates are displayed in Table 3. These estimates quantify the relative significance of the given metric in explaining the overall network structure; and, more specifically, they specify how much the log odds of an edge existing increases for each unit increase in the corresponding metric. For example, the final graphical GOF model for subject 10 shows that GWESP is the most important metric (other than the number of edges) in describing the structure of the subject's network given the larger absolute value of the parameter estimate. Additionally, the positiveness of the estimate associated with GWESP indicates that an edge that closes a triangle is more likely to exist than it would by chance (i.e., the network has more clustering than a random network where the probability of an edge is ) for the family of networks represented by subject 10's fitted model. As evidenced by the results in Table 3, the three model selection methods can lead to very different “best” models. The disparate final model GOF plots that can result from the three different model selection approaches are exhibited in Figures 3 and 4 (for subjects 2 and 8). Again, our aim here is not to judge the three selection methods, but to highlight the fact that they can lead to disparate final models/sets of features. These model selection approaches have been used seemingly arbitrarily in the literature; and, to our knowledge, no detailed comparisons have been performed to determine whether the approaches generally produce the same “best” model/set of features. For our purposes we recommend the graphical GOF approach as the standard and will use it in future analyses given that our main scientific interest lies in being able to capture and reproduce the structure of the fitted brain networks. With the exception of subject 8, the graphical GOF approach produces reasonably good fits for all subjects. The remaining best graphical selection model GOF plots are shown in Figures 5–12.

[Figure omitted. See PDF.]

Figure 3. Goodness-of-fit plots for the final models for subject 2.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks. (a) P-value model. (b) AIC model. (c) Graphical model.

https://doi.org/10.1371/journal.pone.0020039.g003

[Figure omitted. See PDF.]

Figure 4. Goodness-of-fit plots for the final models for subject 8.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks. (a) P-value model. (b) AIC model. (c) Graphical model.

https://doi.org/10.1371/journal.pone.0020039.g004

[Figure omitted. See PDF.]

Figure 5. Goodness-of-fit plots for the final graphical selection model for subject 3.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks.

https://doi.org/10.1371/journal.pone.0020039.g005

[Figure omitted. See PDF.]

Figure 6. Goodness-of-fit plots for the final graphical selection model for subject 5.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks.

https://doi.org/10.1371/journal.pone.0020039.g006

[Figure omitted. See PDF.]

Figure 7. Goodness-of-fit plots for the final graphical selection model for subject 9.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks.

https://doi.org/10.1371/journal.pone.0020039.g007

[Figure omitted. See PDF.]

Figure 8. Goodness-of-fit plots for the final graphical selection model for subject 10.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks. This exemplifies a good fitting model.

https://doi.org/10.1371/journal.pone.0020039.g008

[Figure omitted. See PDF.]

Figure 9. Goodness-of-fit plots for the final graphical selection model for subject 12.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks.

https://doi.org/10.1371/journal.pone.0020039.g009

[Figure omitted. See PDF.]

Figure 10. Goodness-of-fit plots for the final graphical selection model for subject 13.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks.

https://doi.org/10.1371/journal.pone.0020039.g010

[Figure omitted. See PDF.]

Figure 11. Goodness-of-fit plots for the final graphical selection model for subject 16.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks.

https://doi.org/10.1371/journal.pone.0020039.g011

[Figure omitted. See PDF.]

Figure 12. Goodness-of-fit plots for the final graphical selection model for subject 21.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks.

https://doi.org/10.1371/journal.pone.0020039.g012

[Figure omitted. See PDF.]

Table 3. Final model estimates by model selection approach for each subject.

https://doi.org/10.1371/journal.pone.0020039.t003

Despite the obvious importance of Edges (as evidenced by the absolute values of its parameter estimates in Table 3) in the models, the overlap between the simulated and observed networks in the GOF plots is not merely an effect of pure connectivity, but also an effect of network organization. As mentioned in the Materials and Methods section, an ERGM with just an Edges metric is equivalent to the Erdos-Renyi random graph. Thus, due to the small worldness of brain networks, models of this type will not capture the tight local clustering/regional specificity (among other properties) present in these networks [5]. Figure 13 illustrates this point by exhibiting the disparate GOF plots for an Edges only model and the final graphical selection model for subject 10. Clearly the Edges only model is unable to capture the regional specificity (edge-wise shared partners distribution) and global processing (minimum geodesic distance distribution) properties of brain networks which are well embodied by the final graphical selection model.

[Figure omitted. See PDF.]

Figure 13. Goodness-of-fit plots for the Edges only and final graphical selection models for subject 10.

The vertical axis is the logit of relative frequency, the solid lines represent the statistics of the observed network, and the boxplots represent the distributions of the 100 simulated networks. (a) Edges only model. (b) Final graphical model.

https://doi.org/10.1371/journal.pone.0020039.g013

The bolded explanatory metrics in Table 3 are those contained in at least half of the “best” subject network models based on the graphical GOF approach. Examining the uniformity of the selected explanatory metrics across subjects in this way is needed for the development of a “best assessment” ERGM for the reasons detailed in the Introduction and Materials and Methods sections. Examination of these metrics leads to an overall ERGM for whole-brain networks that requires a Connectedness metric (Edges), a Local Efficiency metric (GWESP), and a Global Efficiency metric (GWNSP). That is,(2)These three metrics having the most influential impact on overall functional brain network organization in these subjects seems consistent with our biological understanding of the brain. The number of functional connections present (Edges) is clearly instrumental in information transfer while also playing a role in brain network organization [6]. Clustering (GWESP) is another critical feature of brain network architecture that allows the efficient local processing of information. The consistently negative values associated with GWNSP indicate that if two brain areas are not functionally connected, they are less likely to have shared connections with other regions than they would by chance. That is, two regions are less likely than by chance to have a 2-path as the shortest path between them. Speculatively, this may result from the brain having direct connections when necessary, but allowing for slightly longer global connections (3-paths, etc.) to maintain efficiency otherwise. Additionally, the synergistic combination of these metrics engenders networks that well capture the geodesic (global efficiency), shared partner (local efficiency), degree, and triad census (motifs) distributions of brain networks as evidenced by several of the GOF plots in Figures 3–13.

Group-based network comparisons can potentially be performed by comparing the mean of the estimated , , and values among groups via hypothesis testing or classification techniques. It is important to note that if one were to just compare the mean of the estimated (Edges) values among groups, for instance, potential confounding from the GWESP and GWNSP would be inherently accounted for given that the estimates account for all other metrics in the model. In the hypothesis testing framework one can exploit the fact that the 's are approximate MLEs and thus asymptotically have a Gaussian distribution. Approximate T-tests and/or F-tests can then be employed. Investigating the individual differences in final models among subjects is also important. Although parameter values cannot be directly compared when different models are fitted, the disparate fits themselves may elucidate biologically interesting differences among groups or individual subjects.

Here we implement our best assessment ERGM from equation 2 to illustrate its utility for comparing groups of networks. The subjects were split into a younger (aged 20–26) and slightly older (aged 29–35) group (5 subjects each) in order to assess if there were any discernible differences between their brain networks. Other studies have shown that older adults tend to have less clustering and slightly more connections than their younger counterparts [32], [33]. However, direct comparisons have not been done on groups of subjects this close in age to establish whether these changes tend to commence immediately or take effect at older ages. Moreover, these studies did not consider the potential confounding effects of other network metrics when assessing these differences. As evidenced by the results of our analysis exhibited in Table 4, the two groups differ significantly in (the GWNSP parameter) with the younger group having a more negative value. That is, if two nodes are not functionally connected, they are more likely to have shared connections with other nodes in the brain networks of the older group. Biologically, this could be the result of the older brain maintaining two-path connections between brain areas that have lost their direct connections; however, this interpretation is purely speculative at this point. Interestingly, there is not a statistically significant difference between the groups for the Edges or GWESP parameter, with the trend being for the older subjects' networks to have more connections and clustering. These findings run counter to those in the literature and may stem from the fact that our analysis accounts for some of the confounding that arises from network metric dependencies [6], [7]. These disparate findings could also just be a result of the closeness in age of the two groups or random variability given our small sample size. As noted by [7], larger and methodologically more comparable future investigations are needed to resolve many of the contradictory findings in functional connectivity studies.

[Figure omitted. See PDF.]

Table 4. Results of ERGM parameter estimate comparisons between younger and older subjects.

https://doi.org/10.1371/journal.pone.0020039.t004

In addition to model representation and comparison, ERGMs also provide a statistically sound method for simulating complex brain networks as is done for the GOF plots. To illustrate their utility in this context we simulated 100 networks based on the fitted ERGM of subject 10. We then calculated several descriptive metrics commonly used in the neuroimaging literature for the observed and simulated networks to assess the utility of the simulated networks within the neuroscientific context. Table 5 displays the results of these computations for Clustering coefficient (), Characteristic path length (), Local Efficiency (), Global Efficiency (), and Mean Nodal Degree () (see [4], [34] for details on these metrics). As evidenced by the results in this table, the simulated networks are very similar to the observed network. Hence ERGMs render an approach to simulating scientifically meaningful brain networks.

[Figure omitted. See PDF.]

Table 5. Network metrics of observed and simulated networks from subject 10.

https://doi.org/10.1371/journal.pone.0020039.t005

Discussion

Our analyses in the previous section illustrate the utility of ERGMs for modeling, analyzing, and simulating complex whole-brain networks. We have also provided a foundation for the development of a best assessment ERGM for the classification and comparison of brain networks via the evaluation of three model/feature selection approaches. The graphical GOF approach serves as the best method given the scientific interest in being able to capture and reproduce the structure of the fitted networks. The greatest appeal of modeling brain networks with ERGMs lies in their ability to efficiently represent this complex network data and allow examining the way in which a network's global structure and function depend on local structural components.

There are a myriad of ways in which ERGMs can potentially be useful for brain network researchers. As previously discussed and demonstrated, groups of networks can be statistically compared and classified (by disease status, age, task, etc.) based on several network features simultaneously. The models also provide a way of exploring which local features of brain networks are most important in explaining their global architecture. As noted by many authors [35]–[38], an analysis approach that can capture the network characteristics from a group of subjects' brain networks is needed. ERGMs provide a potential solution since one could average the parameter profiles, , of a group and then simulate “representative” networks based on this averaged profile. Preliminary work has shown this approach to be quite effective. These representative networks can serve as null networks against which other networks and network models can be compared, as visualization tools, and as a means for characterizing properties of network metrics in a group (e.g., community structure). ERGMs, in general, will also serve to both accommodate the ever increasing complexity of whole-brain analyses and inform future statistical models for whole-brain research.

A computational limitation of note for brain network researchers is that MCMC MLE fits of ERGMs can be computationally intensive and may fail to converge with more spatially resolved networks than the 90 ROI ones used here. This fitting algorithm has been shown to handle networks of several thousand nodes [17]; however, its effectiveness is more dependent on the number and topological structure of the edges than the node count [23]. Future work will examine the scalability of ERGMs fitted with MCMC MLE in the context of brain networks. As convergence issues arise with more finely parcellated networks, MPLE fits may serve as an appropriate alternative [16].

Another potential issue of note is that the original data's variability may affect the resulting ERGM fits. A given subject may exhibit variability of the connections in their brain networks at different times of day due to experimental or physiological reasons. Assessment of the robustness of ERGM fits to this within-subject variability is important and will be the focus of future investigations.

In addition to the utility of ERGMs in the research context, the potential implications of their use in the clinical context are profound as they can aid in elucidating system level functional features/neurological processes (represented by the explanatory network metrics) that play a role in various cognitive disorders. For instance, several authors have shown that schizophrenics have less local efficiency in their brain networks [7], [35], [36] (which would correspond to a smaller parameter estimate for GWESP in the ERGM framework) than control subjects. ERGMs enable empirically examining how this difference in efficiency affects global brain structure and comparing these emergent whole-brain brain networks between schizophrenics and controls. For example, one could simulate networks based on model fits to schizophrenics and controls to see how this difference affects the variability of the resulting networks. This comparison may give us insight into the neurological mechanisms that lead to schizophrenia (e.g., lack of local neuronal communication leads to less stability in global structure for schizophrenics).

Aside from the aforementioned clinical and biological work that can be done with the models, there are also many possible directions for future methodological research involving the analysis of complex brain networks with ERGMs. Approximating the small-sample distribution of may prove useful for hypothesis testing frameworks in which appealing to asymptotic normality may not be appropriate. Developing methods for quantifying GOF plots to remove subjectivity and allow for analytical comparisons of the graph will be valuable. The approach should allow some flexibility in determining how to weight the four comparison statistics with respect to their relative importance to the scientific context. Developing novel explanatory network metrics rooted in both the biology of the brain and the mathematics of ERGMs will engender better best assessment models for network comparison. A corresponding hybrid model selection approach where models are penalized for using many covariates and the GOF plots are assessed will prove useful in maintaining parsimony as the number of relevant explanatory metrics increases. The extension of ERGMs to directed and/or weighted brain networks will prove beneficial as construction of these network types gains feasibility.

Supporting Information

[Figure omitted. See PDF.]

Appendix S1.

Model Selection

https://doi.org/10.1371/journal.pone.0020039.s001

(PDF)

Acknowledgments

We thank the editor and referees for their comments that considerably improved the paper. An earlier version of this manuscript can be found at arxiv.org (Simpson et al., arXiv:1007.3230v1 [stat.AP]).

Author Contributions

Conceived and designed the experiments: SLS SH PJL. Performed the experiments: SLS SH PJL. Analyzed the data: SLS SH PJL. Wrote the paper: SLS SH PJL.

Citation: Simpson SL, Hayasaka S, Laurienti PJ (2011) Exponential Random Graph Modeling for Complex Brain Networks. PLoS ONE6(5): e20039. https://doi.org/10.1371/journal.pone.0020039

References

1. Iturria-Medina Y, Sotero RC, Canales-Rodriguez EJ, Aleman-Gomez Y, Melie-Garcia L (2008) Studying the human brain anatomical network via diffusion-weighted MRI and graph theory. Neuro Image 40: 1064–1076.Y. Iturria-MedinaRC SoteroEJ Canales-RodriguezY. Aleman-GomezL. Melie-Garcia2008Studying the human brain anatomical network via diffusion-weighted MRI and graph theory.Neuro Image4010641076

2. Gong G, He Y, Concha L, Lebel C, Gross DW, et al. (2009) Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. Cerebral Cortex 19: 524–536.G. GongY. HeL. ConchaC. LebelDW Gross2009Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography.Cerebral Cortex19524536

3. Bassett DS, Bullmore E (2006) Small-world brain networks. Neuroscientist 12: 512–523.DS BassettE. Bullmore2006Small-world brain networks.Neuroscientist12512523

4. Stam CJ, Reijneveld JC (2007) Graph theoretical analysis of complex networks in the brain. Nonlinear Biomedical Physics 1: 3.CJ StamJC Reijneveld2007Graph theoretical analysis of complex networks in the brain.Nonlinear Biomedical Physics13

5. Watts DJ, Strogatz SH (1998) Collective dynamics of small-world networks. Nature 393: 440–442.DJ WattsSH Strogatz1998Collective dynamics of small-world networks.Nature393440442

6. van Wijk BCM, Stam CJ, Daffertshofer A (2010) Comparing brain networks of different size and connectivity density using graph theory. PLoS ONE 5(10): e13701.BCM van WijkCJ StamA. Daffertshofer2010Comparing brain networks of different size and connectivity density using graph theory.PLoS ONE510e13701

7. Lynall ME, Bassett DS, Kerwin R, McKenna PJ, Kitzbichler M, et al. (2010) Functional Connectivity and brain networks in schizophrenia. The Journal of Neuroscience 30(28): 9477–9487.ME LynallDS BassettR. KerwinPJ McKennaM. Kitzbichler2010Functional Connectivity and brain networks in schizophrenia.The Journal of Neuroscience302894779487

8. Frank O, Strauss D (1986) Markov graphs. Journal of the American Statistical Association 81: 832–842.O. FrankD. Strauss1986Markov graphs.Journal of the American Statistical Association81832842

9. Pattison PE, Wasserman S (1999) Logit models and logistic regressions for social networks: II. Multivariate relations. British Journal of Mathematical and Statistical Psychology 52: 169–194.PE PattisonS. Wasserman1999Logit models and logistic regressions for social networks: II. Multivariate relations.British Journal of Mathematical and Statistical Psychology52169194

10. Robins GL, Pattison PE, Wasserman S (1999) Logit models and logistic regression for social networks: III. Valued relations. Psychometrika 64: 371–394.GL RobinsPE PattisonS. Wasserman1999Logit models and logistic regression for social networks: III. Valued relations.Psychometrika64371394

11. Wasserman S, Pattison PE (1996) Logit models and logistic regressions for social networks: I. An introduction to Markov graphs and p*. Psychometrika 61: 401–425.S. WassermanPE Pattison1996Logit models and logistic regressions for social networks: I. An introduction to Markov graphs and p*.Psychometrika61401425

12. Robins GL, Pattison PE, Kalish Y, Lusher D (2007) An introduction to exponential random graph (p*) models for social networks. Social Networks 29: 173–191.GL RobinsPE PattisonY. KalishD. Lusher2007An introduction to exponential random graph (p*) models for social networks.Social Networks29173191

13. Robins GL, Snijders T, Wang P, Handcock M, Pattison PE (2007) Recent developments in exponential random graph (p*) models for social networks. Social Networks 29: 192–215.GL RobinsT. SnijdersP. WangM. HandcockPE Pattison2007Recent developments in exponential random graph (p*) models for social networks.Social Networks29192215

14. Bullmore E, Sporns O (2009) Complex brain networks: graph theoretical analysis of structural and functional systems. Nature Reviews Neuroscience 10: 186–198.E. BullmoreO. Sporns2009Complex brain networks: graph theoretical analysis of structural and functional systems.Nature Reviews Neuroscience10186198

15. Meunier D, Achard S, Morcom A, Bullmore E (2009) Age-related changes in modular organization of human brain functional networks. Neuro Image 44: 715–723.D. MeunierS. AchardA. MorcomE. Bullmore2009Age-related changes in modular organization of human brain functional networks.Neuro Image44715723

16. Saul ZM, Filkov V (2007) Exploring biological network structure using exponential random graph models. Bioinformatics 23: 2604–2611.ZM SaulV. Filkov2007Exploring biological network structure using exponential random graph models.Bioinformatics2326042611

17. Hunter DR, Goodreau SM, Handcock MS (2008) Goodness of fit of social network models. Journal of the American Statistical Association 103: 248–258.DR HunterSM GoodreauMS Handcock2008Goodness of fit of social network models.Journal of the American Statistical Association103248258

18. Peiffer AM, Hugenschmidt CE, Maldjian JA, Casanova R, Srikanth R, et al. (2009) Aging and the interaction of sensory cortical function and structure. Human Brain Mapping 30: 228–240.AM PeifferCE HugenschmidtJA MaldjianR. CasanovaR. Srikanth2009Aging and the interaction of sensory cortical function and structure.Human Brain Mapping30228240

19. Tzourio-Mazoyer N, Landeau B, Papathanassiou D, Crivello F, Etard O, et al. (2002) Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuro Image 15: 273–289.N. Tzourio-MazoyerB. LandeauD. PapathanassiouF. CrivelloO. Etard2002Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain.Neuro Image15273289

20. van den Heuvel MP, Stam CJ, Boersma M, Hulshoff Pol HE (2008) Small-world and scale free organization of voxel-based resting-state functional connectivity in the human brain. Neuro Image 43: 528–539.MP van den HeuvelCJ StamM. BoersmaHE Hulshoff Pol2008Small-world and scale free organization of voxel-based resting-state functional connectivity in the human brain.Neuro Image43528539

21. Hayasaka S, Laurienti PJ (2010) Comparison of characteristics between region- and voxel based network analysis in resting-state fMRI. Neuro Image 50: 499–508.S. HayasakaPJ Laurienti2010Comparison of characteristics between region- and voxel based network analysis in resting-state fMRI.Neuro Image50499508

22. Handcock MS (2002) Statistical models for social networks: Inference and degeneracy. In: Breiger R, Carley K, Pattison PE, editors. Dynamic Social Network Modelling and Analysis: Workshop Summary and Papers. Washington, DC: National Academy Press. pp. 229–240.MS Handcock2002Statistical models for social networks: Inference and degeneracy.R. BreigerK. CarleyPE PattisonDynamic Social Network Modelling and Analysis: Workshop Summary and Papers.Washington, DCNational Academy Press229240

23. Handcock MS, Hunter DR, Butts CT, Goodreau SM, Morris M (2008) statnet: Software tools for the statistical modeling of network data. MS HandcockDR HunterCT ButtsSM GoodreauM. Morris2008statnet: Software tools for the statistical modeling of network data.URL http://statnetproject.org. URL http://statnetproject.org.

24. Snijders TAB, Pattison PE, Robins GL, Handcock MS (2006) New specifications for exponential random graph models. Sociological Methodology 36: 99–154.TAB SnijdersPE PattisonGL RobinsMS Handcock2006New specifications for exponential random graph models.Sociological Methodology3699154

25. Robins G, Pattison P, Wang P (2009) Closure, connectivity and degree distributions: exponential random graph (p*) models for directed social networks. Social Networks 31: 105–117.G. RobinsP. PattisonP. Wang2009Closure, connectivity and degree distributions: exponential random graph (p*) models for directed social networks.Social Networks31105117

26. Rinaldo A, Fienberg SE, Zhou Y (2009) On the geometry of discrete exponential families with application to exponential random graph models. Electronic Journal of Statistics 3: 446–484.A. RinaldoSE FienbergY. Zhou2009On the geometry of discrete exponential families with application to exponential random graph models.Electronic Journal of Statistics3446484

27. Morris M, Handcock MS, Hunter DR (2008) Specification of exponential-family random raph models: terms and computational aspects. Journal of Statistical Software 24: M. MorrisMS HandcockDR Hunter2008Specification of exponential-family random raph models: terms and computational aspects.Journal of Statistical Software24

28. van Duijn MAJ, Gile KJ, Handcock MS (2009) A framework for the comparison of maximum pseudo-likelihood and maximum likelihood estimation of exponential family random graph models. Social Networks 31: 52–62.MAJ van DuijnKJ GileMS Handcock2009A framework for the comparison of maximum pseudo-likelihood and maximum likelihood estimation of exponential family random graph models.Social Networks315262

29. Hunter DR, Handcock MS, Butts CT, Goodreau SM, Morris M (2008) ergm: A package to fit, simulate and diagnose exponential-family models for networks. Journal of Statistical Software 24: DR HunterMS HandcockCT ButtsSM GoodreauM. Morris2008ergm: A package to fit, simulate and diagnose exponential-family models for networks.Journal of Statistical Software24

30. Muller KE, Fetterman BA (2002) Regression and ANOVA: an integrated approach using SAS software. Cary, , NC: SAS Institute Inc. KE MullerBA Fetterman2002Regression and ANOVA: an integrated approach using SAS softwareCary, , NCSAS Institute Inc

31. Akaike H (1974) A new look at the statistical model identification. IEEE Transaction on Automatic Control AC-19: 716–723.H. Akaike1974A new look at the statistical model identification.IEEE Transaction on Automatic ControlAC-19716723

32. Gaal ZA, Boha R, Stam CJ, Molnar M (2009) Age-dependent features of EEG-reactivity Spectral, complexity, and network characteristics. Neuroscience Letters 479: 79–84.ZA GaalR. BohaCJ StamM. Molnar2009Age-dependent features of EEG-reactivity Spectral, complexity, and network characteristics.Neuroscience Letters4797984

33. Gong G, Rosa-Neto P, Carbonell F, Chen ZJ, He Y, Evans AC (2009) Age- and gender related differences in the cortical anatomical network. Journal of Neuroscience 29: 15684–15693.G. GongP. Rosa-NetoF. CarbonellZJ ChenY. HeAC Evans2009Age- and gender related differences in the cortical anatomical network.Journal of Neuroscience291568415693

34. Rubinov M, Sporns O (2010) Complex network measures of brain connectivity: uses and interpretations. Neuroimage 52(3): 1059–1069.M. RubinovO. Sporns2010Complex network measures of brain connectivity: uses and interpretations.Neuroimage52310591069

35. Rubinov M, Knock SA, Stam CJ, Micheloyannis S, Harris AWF, et al. (2007) Small-world properties of nonlinear brain activity in schizophrenia. Human Brain Mapping 30: 403–416.M. RubinovSA KnockCJ StamS. MicheloyannisAWF Harris2007Small-world properties of nonlinear brain activity in schizophrenia.Human Brain Mapping30403416

36. Alexander-Bloch AF, Gogtay N, Meunier D, Birn R, Clasen L, et al. (2010) Disrupted modularity and local connectivity of brain functional networks in childhood-onset schizophrenia. Frontiers in Systems Neuroscience 4: Article 147.AF Alexander-BlochN. GogtayD. MeunierR. BirnL. Clasen2010Disrupted modularity and local connectivity of brain functional networks in childhood-onset schizophrenia.Frontiers in Systems Neuroscience4Article 147

37. Meunier D, Lambiotte R, Fornito A, Ersche KD, Bullmore ET (2009) Hierarchical modularity in human brain functional networks. Frontiers in Neuroinformatics 3: 37.D. MeunierR. LambiotteA. FornitoKD ErscheET Bullmore2009Hierarchical modularity in human brain functional networks.Frontiers in Neuroinformatics337

38. Joyce KE, Laurienti PJ, Burdette JH, Hayasaka S (2010) A new measure of centrality for brain networks. PLoS ONE 5(8): e12200.KE JoycePJ LaurientiJH BurdetteS. Hayasaka2010A new measure of centrality for brain networks.PLoS ONE58e12200