Abstract
The human cerebral cortex is a complex network of functionally specialized regions interconnected by axonal fibers, but the organizational principles underlying cortical connectivity remain unknown. Here, we report evidence that one such principle for functional cortical networks involves finding a balance between maximizing communication efficiency and minimizing connection cost, referred to as optimization of network costefficiency. We measured spontaneous fluctuations of the blood oxygenation leveldependent signal using functional magnetic resonance imaging in healthy monozygotic (16 pairs) and dizygotic (13 pairs) twins and characterized costefficient properties of brain network functional connectivity between 1041 distinct cortical regions. At the global network level, 60% of the interindividual variance in costefficiency of cortical functional networks was attributable to additive genetic effects. Regionally, significant genetic effects were observed throughout the cortex in a largely bilateral pattern, including bilateral posterior cingulate and medial prefrontal cortices, dorsolateral prefrontal and superior parietal cortices, and lateral temporal and inferomedial occipital regions. Genetic effects were stronger for costefficiency than for other metrics considered, and were more clearly significant in functional networks operating in the 0.09–0.18 Hz frequency interval than at higher or lower frequencies. These findings are consistent with the hypothesis that brain networks evolved to satisfy competitive selection criteria of maximizing efficiency and minimizing cost, and that optimization of network costefficiency represents an important principle for the brain's functional organization.
Introduction
Functional magnetic resonance imaging (fMRI) studies have shown that regional fluctuations of spontaneous brain activity, measured in the absence of an explicit task (the socalled resting state), are highly organized, being correlated across spatially distributed networks in a manner that recapitulates the topography of taskevoked functional coactivation patterns (Fox and Raichle, 2007; Smith et al., 2009). These networks are robust (Damoiseaux et al., 2006), are linked to synchronized oscillations in neuronal activity (He et al., 2008), and represent a major source of variance in taskevoked activation and behavior (Fox et al., 2007), suggesting that they reflect a functionally important aspect of brain organization (Fox and Raichle, 2007).
The nontrivial configuration of spontaneous brain activity suggests that certain principles of network organization may confer adaptive benefits. One proposed principle involves satisfying competitive selection criteria of maximizing communication efficiency and minimizing connection costs (Bullmore and Sporns, 2009). In complex networks such as the brain, communication efficiency is enhanced when the number of connections, or path length, linking any two nodes is low, resulting in faster, more accurate information transmission (Latora and Marchiori, 2001). In anatomical brain networks, a few longrange projections can dramatically increase efficiency (Buzsáki et al., 2004), but this comes at the cost of increased wiring length. In principle, network efficiency could be maximized by adding more connections (Achard and Bullmore, 2007), but the metabolic costs associated with forming and maintaining each brain connection (Laughlin and Sejnowski, 2003) limit the total number sustainable. This balance, between efficiency maximization and connection cost minimization, may be construed as one concerned with optimization of network costefficiency.
Cost and efficiency constraints on anatomical brain networks have been studied widely (Cherniak et al., 2004; Chen et al., 2006; Kaiser and Hilgetag, 2006; Bassett et al., 2010), but there are reasons to believe that costefficiency optimization is also an important principle for functional brain organization. Restingstate functional networks are correlated with underlying anatomical connectivity (Vincent et al., 2007; Honey et al., 2009) and show a smallworld organization (Achard et al., 2006) characterized by high efficiency and low cost (Latora and Marchiori, 2001, 2003). Moreover, costefficiency may index functional network integrity, being affected by aging and disease and positively correlated with cognitive performance (Achard and Bullmore, 2007; Bassett et al., 2009). Consequently, it has been proposed that both functional and structural brain network organization may represent the outcome of evolutionary selection for neural connections conferring high efficiency at low cost (Kaiser and Hilgetag, 2006; Bullmore and Sporns, 2009; Bassett et al., 2010). This hypothesis predicts that individual differences in functional network costefficiency should be under genetic control and therefore demonstrably heritable.
We tested this prediction using restingstate fMRI (rsfMRI) in healthy monozygotic (MZ) and dizygotic (DZ) twins. Twins provide an ideal natural experiment that enables determination of genetic and environmental contributions to a given phenotype. MZ twins have identical genomes, whereas DZ twins share, on average, 50% of their genes. Thus, a phenotype is deemed to be under genetic control, or heritable, if MZ twin covariation exceeds DZ twin covariation.
Materials and Methods
Participants.
Thirtyfour (18 male) monozygotic (MZ) and twentysix (12 male) samesex dizygotic (DZ) psychiatrically and neurologically healthy twins were recruited for scanning. Twin pairs were registrants with the Australian Twin Registry (ATR), and were recruited through mailouts and advertisements. All participants were screened using structured clinical interviews (National Institute of Mental Health, 1992). To confirm zygosity, blood samples extracted from the twins were assessed using an AmpFLSTR Profile Plus PCR Amplification Kit (Applied Biosystems). Amplicon allele size ranges were between 100 and 360 bp for each of the following loci: AMEL (X, Y alleles), D13S317 (8 alleles), D18S51 (23 alleles), D21S11 (24 alleles), D3S1358 (8 alleles), D5S818 (10 alleles), D7S820 (10 alleles), D8S1179 (12 alleles), FGA (28 alleles), and vWA (14 alleles). There were no differences between MZ and DZ twins in terms of age (MZ mean = 37.81 years, SD = 13.62 years; DZ mean = 43.15, SD = 9.91 years; t_{(27)} = −1.16, p = 0.26) or intelligence, as estimated using the Wechsler Adult Reading Test (The Psychological Corporation, 2001) (MZ mean = 107.30, SD = 7.80; DZ mean = 106.88, SD = 6.07; F_{(1,27)} = 0.04, p = 0.84). Data from one male MZ twin pair were excluded due to one member moving excessively in the scanner, resulting in 32 MZ twins (16 pairs) and 26 DZ twins (13 pairs) being included in the final analysis. All participants gave written, informed consent before participating, in accordance with local ethics committee guidelines.
Image acquisition and preprocessing.
Participants were scanned using a 1.5 T Siemens Avanto located at St Vincent's Hospital, Melbourne, Australia. Eight hundred and forty echoplanar imaging volumes providing T2*weighted blood oxygenation leveldependent (BOLD) contrast were acquired with the following parameters: number of slices, 21; voxel dimensions, 3.4 × 3.4 × 4 mm (1 mm interslice gap); flip angle, 78°; echo time, 26 ms; time to repetition, 1420 ms. Participants were instructed to lie quietly in the scanner with their eyes closed without falling asleep. All participants confirmed that they did not fall asleep in subsequent debriefing. T1weighted anatomical scans were also acquired: voxel dimensions = 0.94 × 0.94 × 1.5; flip angle = 15°; echo time = 3.93 ms; time to repetition = 1930 ms.
For each individual, functional volumes were realigned using a rigidbody transformation to correct for geometric displacements associated with head movements and rotations. Temporal motion correction was then performed by regressing the current and lagged first and secondorder displacements against the time series of the realigned images, as implemented in CamBa (http://wwwbmu.psychiatry.cam.ac.uk/software). The residuals of this regression were used for further analysis. Finally, the realigned, temporally corrected images were spatially normalized to the International Consortium for Brain Mapping echoplanar imaging template supplied with SPM5 (http://www.fil.ion.ucl.ac.uk/spm/software/spm5) using a 12parameter affine transformation (Jenkinson and Smith, 2001) (http://www.fmrib.ox.ac.uk/fsl).
Cortical parcellation.
We used a graph analytic approach (Bullmore and Sporns, 2009; Rubinov and Sporns, 2010) to measure the topological properties of each individual's functional brain network. In this framework, the brain is modeled as a graph comprising N nodes connected by M edges. The nodes can be defined at the voxel level (van den Heuvel et al., 2008; Hayasaka and Laurienti, 2010), region level (Achard et al., 2006), or intermediate scales (Hagmann et al., 2007; Fornito et al., 2010; Zalesky et al., 2010). For the current analysis, we parcellated each participant's cortex into 1041 discrete regions using an approach similar to that used by previous authors (Hagmann et al., 2007; Honey et al., 2009). We have previously found that this approximate resolution of analysis provides a good tradeoff between computational speed, spatial resolution, and signaltonoise ratio (Fornito et al., 2010).
To generate the parcellation, each participant's T1weighted anatomical image was used to generate surface (triangular) mesh models of the gray/white and gray/CSF boundaries of the left and right hemisphere cortical ribbons using validated algorithms (Dale et al., 1999; Fischl et al., 1999) implemented in freely available software (http://surfer.nmr.mgh.harvard.edu/fswiki/Home). The white matter surface of each hemisphere was then parcellated using an algorithm that incrementally merged adjacent pairs of triangles forming the white matter surface mesh. The input to the algorithm was the desired number of regions, denoted with N, and the white matter surface triangulation. Let F denote the number of triangles comprising the triangulation. The goal of the algorithm was to achieve a parcellation that maximized the uniformity and compactness of each regions' surface area.
The algorithm comprised two distinct phases. In the first phase, the algorithm was serially performed at a total of F − N merge operations. To complete each merge operation, a pair of adjacent triangles was merged together to form a region. In this sense, a region was simply a set of adjacent triangles. Merge operations were also performed between adjacent pairs of regions, in which case two distinct regions became one. It can be easily shown that performing exactly F − N merge operations ensured that exactly N regions remained at the end of the first phase of the algorithm. In particular, if each triangle is thought of as a distinct region, whenever a pair of regions was merged (regardless of which pair), the total number of regions was reduced by one, and thus after F − N merge operations, N regions must have remained.
While the uniformity in surface area between regions was reasonable, this uniformity was further improved during the second phase of the algorithm, which subdivided any region with a surface area at least twice as large that of the smallest region. This rule guaranteed a strict bound on the uniformity of the final parcellation; namely, the smallest region was guaranteed to be greater than or equal to half the size of the largest region. A result of the second phase of the algorithm was that the final number of regions comprising the parcellation was always slightly greater than the desired number of regions, N. In particular, if M regions had a surface area of size at least twice as large as the smallest region, then each of these regions was subdivided into at least two separate regions, and hence the total number of regions comprising the final parcellation was at least N + M. Note that a subdivided region was itself subdivided if its surface area was still at least twice as large as the surface area of the smallest region.
In the present study, the algorithm was executed with N = 500 regions per hemisphere. At the completion of phase 2, the left hemisphere comprised 596 regions, while the right comprised 603. To sample regional activity from the functional volumes, the surfacebased parcellations were converted into volumetric representations and projected out to fill the gray matter of the entire cortical ribbon, defined via atlasbased tissue segmentation of the T1weighted image (Fischl et al., 2002). This volumebased parcellation was dilated by one voxel to ensure adequate coverage of the cortical mantle. Each participant's functional volumes were then registered to their anatomical image via an affine transformation (Jenkinson and Smith, 2001). To minimize noise introduced by variable EPI coverage across individuals, we excluded ROIs containing <70% nonzero intensity voxels, resulting in a final template comprising 1041 cortical regions across both hemispheres (Fig. 1A).
Temporal filtering.
The custom template was used to extract the mean time series of each of the 1041 defined cortical regions from each participant's motioncorrected, spatially normalized functional volumes. These mean regional time series were then decomposed into four distinct frequency intervals via the maximal overlap discrete wavelet transform (Percival and Walden, 2000), which is well suited to nonstationary, 1/ftype signals such as fMRI measurements (Bullmore et al., 2004). The four frequency intervals (wavelet scales) were as follows: scale 1, 0.18–0.35 Hz; scale 2, 0.09–0.18 Hz; scale 3, 0.04–0.09 Hz; and scale 4, 0.02–0.04 Hz (Fig. 1B). We then used the CompCor technique to correct the waveletfiltered regional time series for fluctuations in nonneuronal physiological processes (e.g., fluctuations in cardiac and respiratory rates) (Behzadi et al., 2007). Briefly, this approach involved the following steps for each participant: (1) identification of voxels showing the highest temporal standard deviation (top 2%) in the functional data; (2) principal component analysis of these voxels' time courses; (3) identification, via permutation testing, of the components accounting for a significant proportion of variance in the data; (4) extraction of reference time courses for each significant principal component; and (5) correction, via linear regression, of each regional time series for covariance with the principal component time courses identified in step 4. The results are comparable to those obtained using modelbased correction procedures requiring physiological monitoring (Behzadi et al., 2007).
Network characterization.
To estimate the costefficiency of each individual's cortical functional network, analysis of regional time series progressed according to three broad stages: (1) estimation of interregional functional connectivity; (2) generation of a graphbased representation of network topology; and (3) topological analysis of the network. The details of each of these steps are provided in the following. Detailed accounts of basic principles underlying graph analysis and its application to neuroimaging data have been published previously (Albert and Barabási, 2002; Newman, 2003; Bullmore and Sporns, 2009; Rubinov and Sporns, 2010).
Functional connectivity estimation.
To generate cortical functional networks for each participant, correlations between the waveletfiltered, noisecorrected time series of each possible pair of the 1041 cortical regions sampled were computed separately for each frequency interval (wavelet scale), resulting in four {1041 × 1041} functional connectivity matrices for each individual (Fig. 1C). A considerable range of correlation values, extending from <0.01 to >0.80, was observed at all frequency intervals (supplemental Fig. S1, available at www.jneurosci.org as supplemental material). Average functional connectivity was greater at lower frequencies, consistent with prior observations (Cordes et al., 2001; Achard et al., 2006), but there was also a substantial increase in the estimation error (computed according to Whitcher et al., 2000) associated with lower frequency correlations: from scales 1 to 4, there was an average ∼3fold increase in mean estimation error, compared to an ∼1.1fold increase in average connectivity (supplemental Fig. S1, available at www.jneurosci.org as supplemental material). The increase in estimation error likely reflects the greater statistical difficulties associated with precisely estimating lowfrequency properties from relatively short time series (Achard et al., 2008; see also Discussion).
Graphbased representations of brain network topology.
The topology of a complex network can be efficiently represented as a graph of N nodes connected by M edges. In the case of brain networks, the nodes represent different brain regions, and the edges represent some measure of interaction between them. In the current analysis, the brain regions represented the 1041 cortical areas from which mean regional fMRI time series were sampled, and the edges represented the frequencyspecific wavelet correlations between these time series, computed as described in the previous section.
When networks are generated from continuous association measures, such as wavelet correlations, it is customary to apply a threshold to remove spurious edges and emphasize key topological properties of the network. As this threshold can be arbitrary, we examined network measures across the full range of possible thresholds representing connection densities 5% < k < 40%; i.e., from 5% to 40% of all possible connections (Fig. 1D). The 5% lower bound was chosen to avoid excess network fragmentation at lower costs; the 40% upper bound was chosen because topological properties of brain networks tend toward randomness at higher costs due to the inclusion of potentially spurious functional associations (Bassett et al., 2008).
At each connection density, the graphs were binarized such that suprathreshold correlations were set to 1 and subthreshold correlations to 0, resulting in a series of unweighted, undirected adjacency matrices representing network connectivity at that threshold [for a similar approach, see Achard et al. (2006) and Achard and Bullmore (2007)]. Graph models were generated from these adjacency matrices by drawing an edge between two nodes if and only if they were linked by a suprathreshold connection (Fig. 1E).
Topological analysis.
Once a graphbased representation has been obtained, numerous measures describing different topological properties of the network can be computed (Bullmore and Sporns, 2009; Rubinov and Sporns, 2010). In this study, we were principally interested in characterizing network costefficiency. Four measures were integral to our definition of network costefficiency: global communication efficiency (E_{global}), regional communication efficiency (E_{regional}), connection distance (D_{global}), and connection density (k). The connection density, k, simply represented the proportion of edges present in the thresholded graph relative to the maximum possible 541,320 connections, i.e., the percentage of suprathreshold edges. The global efficiency of the graph G was calculated following the definition given by Latora and Marchiori (2003): where N is the number of nodes in G and L_{ij} is the shortest path length between nodes i and j. Similarly, the efficiency of each node, i, was calculated as follows: In general, E_{global} increases as more connections are added to the network (i.e., as k increases), until it reaches its maximum value when the network is completely connected.
Prior work estimating connection costs in functional networks has simply taken the number of connections present in the network as an index of cost (Achard and Bullmore, 2007; Bassett et al., 2009). However, in physically embedded networks such as the brain, the cost of connecting different nodes scales in proportion to the spatial distance between them (Latora and Marchiori, 2003), and introducing spatial information into the calculation of network costs can alter the results obtained when compared with simple binary measures (Kaiser and Hilgetag, 2006). In this regard, the net connection cost of a network may be viewed as arising from the costs associated with making connections, measured by the connection density of the network, k, and the summed physical distance of these connections, which we refer to as the connection distance, D. In the absence of a complete map of interregional human brain connectivity that would allow derivation of accurate anatomical estimates of D, we followed the approach of Kaiser and Hilgetag (2006) and calculated D as the summed Euclidean distance between regions connected by a suprathreshold connection. Thus, global connection cost was calculated as follows:
where d_{ij} =
Previous studies taking the difference between efficiency and binary connection cost as an index of network costefficiency have found that this difference shows a clear peak at a specific connection density, beyond which the costs associated with adding more connections to the network outweigh the benefits of increased efficiency (Achard and Bullmore, 2007; Bassett et al., 2009). A similar relationship was observed for E_{global} and D_{global} in our data. Figure 2, A and B, shows how these two measures scale as a function of connection density, k, at each wavelet scale. Both measures increase as more connections are added to the network, but the increase in E_{global} is much steeper at sparser densities (k < 30%). Due to this discrepancy, taking the difference between E_{global} and D_{global} as a simple index of costefficiency (Achard and Bullmore, 2007; Bassett et al., 2009) reveals a concave curve with a clear peak, max(E_{global} − D_{global}) (see Fig. 2C). The height of this peak, max(E_{global} − D_{global}), and the connection density at which it occurred, k_{max}, varied across individuals (see Fig. 2C, inset). Thus, for each individual, the network defined by k_{max} represented his or her most costefficient network configuration: the pivotal connection density at which the cost associated with adding more connections to the network outweighed the benefits of increased efficiency. We therefore defined individual differences in global network costefficiency as follows: Regional costefficiency was similarly calculated as follows: where E_{regional} and D_{regional} were calculated from the network defined by each individual's k_{max}.
We normalize by k_{max} in Equations 5 and 6 to compensate for differences in the connection densities associated with each individual's maximally costefficient network, and ensure that when all else was equal, networks with fewer connections were considered to be more costefficient. Thus, by this definition, network connection costs arise from both the distance and number of connections in the network. By exploiting the behavior of the E_{global} − D_{global} as a function of k, the measure provides a datadriven, individually tailored approach to threshold selection because each participant's network is thresholded based on his or her unique maximally costefficient network configuration. This approach is therefore more sensitive to individual differences in network organization, which is likely to be an important consideration in twin studies of brain function (Koten et al., 2009). Related measures based on this individually tailored approach have shown stronger correlations with cognitive performance than other network properties (Bassett et al., 2009), supporting the behavioral relevance of our costefficiency metric.
We note that the interpretation of connection costs in functional networks can be somewhat ambiguous. In general, the operational costs of a network can be defined in terms of its metabolic requirements, and for a brain network these requirements will scale with the amount of wiring required to connect the network (Laughlin and Sejnowski, 2003). In anatomical networks, connection distance can be clearly related to the wiring cost of the network, since wiring cost scales with the total fiber length used to connect the network. In functional networks, the interpretation is less clear because functional connectivity between two regions may emerge in the absence of a direct anatomical link, presumably through polysynaptic pathways (Vincent et al., 2007; Honey et al., 2009). Nonetheless, it is reasonable to assume that greater wiring and metabolic resources are required to propagate information between two regions that are farther apart than between two spatially adjacent areas (Laughlin and Sejnowski, 2003), whether they are connected directly or through multiple pathways. Thus, connection distance can still be viewed as a proxy estimate of the connection cost of a functional network. Future work integrating fMRI with diffusionweighted imaging (Skudlarski et al., 2008; Honey et al., 2009; Zalesky and Fornito, 2009) will facilitate more precise characterization of functional connection costs.
In addition to network costefficiency, we examined genetic influences on several other commonly studied topological properties of brain networks, assessed across the range of connection densities 5% < k < 40%, in 5% increments. These included the component measures making up our costefficiency index—E_{global} and D_{global}, and E_{global} − D_{global}—in addition to several other network properties commonly studied in the literature: local efficiency (E_{local}), normalized path length (λ), normalized clustering coefficient (γ), and network smallworldness (σ). E_{local} is a measure of local information processing or fault tolerance, and is computed as the efficiency of the subgraph defined by an index node's neighbors after removal of that node (Latora and Marchiori, 2003). To estimate λ, an index of the global integration of the network, the mean minimum path length between nodes in the observed networks was computed and normalized by the corresponding value computed using a comparable random network. To estimate γ, the clustering coefficient, a measure of the local cliquishness of network connectivity (i.e., the mean connectedness between each node's neighbors), was calculated in the observed networks and also normalized by its corresponding value in a matched random network. For the normalization, we used the average measure computed from 20 random graphs generated using a rewiring algorithm that ensures matched network size, connection density, and degree distribution (Maslov and Sneppen, 2002). Smallworldness was assessed using the scalar summary γ/λ=σ (Humphries et al., 2006). In smallworld networks, λ ≈ 1 and γ > 1, so a nonrandom, smallworld topology is apparent when σ > 1. Formal definitions for each of these measures can be found elsewhere (Rubinov and Sporns, 2010).
Structural equation modeling.
Similarity between twin pairs (i.e., twin covariance) was assessed using the intraclass correlation coefficient, computed separately for MZ (r_{MZ}) and DZ (r_{DZ}) twins. As a general rule, genetic influences on a given trait are apparent if r_{MZ} is substantially greater than r_{DZ}. However, we used structural equation modeling to more formally quantify the proportion of variance in each studied phenotype attributable to additive genetic (A), shared environmental (C), and nonshared environmental (E) effects (Neale et al., 2003) (http://www.vcu.edu/mx/). In this context, A, C, and E represent latent variables, and their variances (a^{2}, c^{2}, and e^{2}, respectively) reflect the proportion of phenotypic variation attributable to their effects. Model path coefficients were estimated using maximumlikelihood fitting. We first estimated parameters for the full ACE model, and then examined the fit of nested submodels comprising AE, CE, and E variance components. Model fit was indexed by −2 times the loglikelihood, which follows a χ^{2} distribution. However, to avoid inappropriate distributional assumptions on model fits and parametric effects in this relatively modest sample, we used permutation testing for statistical inference (see below). In this context, p > 0.05 for a given reduced model indicates a fit that is not significantly different from the full ACE model, and should thus be accepted under the rule of parsimony. The significance of A was tested with the p value for the difference in fit between CE and ACE models; the significance of C was tested by the p value for the difference in fit between the AE and ACE models (Christian et al., 1995). In both cases, p < 0.05 indicated a significant effect of that parameter. This modeling procedure assumes that twin covariances are positive, and that r_{MZ} > r_{DZ}, consistent with a genetic influence on the phenotype under consideration. These assumptions were not met for all the measures studied in our analyses, and so models were not estimated in these cases.
Inferential statistics.
At each wavelet scale, or frequency interval, for each phenotypic measure, twin pairings and zygosity (MZ or DZ) labels were permuted 5000 times, with the only constraint being that the original pairings were not broken. Let T_{i}_{1} and T_{i}_{2} denote the first and second members of twin pair i. If under random reassignment T_{i}_{1} became paired with T_{j}_{1}, then T_{i}_{2} would be paired with T_{j}_{2}. Permuting twin pairs in this way ensured exchangeability of observations. At each iteration, the full ACE model, in addition to nested submodels was fitted to the permuted data. The relevant model fit statistics were then extracted to create an empirical null distribution for each parameter. The p value of each effect was given by the percentile of the observed statistic on its corresponding null distribution. Ninetyfive percent confidence intervals (CIs) for each parameter were generated using a bootstrapping procedure whereby MZ and DZ pairs were randomly resampled with replacement to produce a surrogate dataset from which genetic parameters were reestimated. This procedure was iterated 1000 times to yield a sampling distribution for each parameter, from which CIs were calculated according to established procedures (Efron, 1987).
Regional genetic effects were thresholded using a permutationbased, clusterwise procedure. For each of the 1041 regions, we permuted the data 1000 times and reestimated the genetic models, leading to a total of 1.041 × 10^{6} permuted datasets. At each permutation, clusters of regions showing genetic effects were formed by merging adjacent nodes surviving a primary clusterforming threshold of p < 0.05. To obtain p values for each cluster that were corrected for multiple comparisons, we followed the procedure advocated by Nichols and Holmes (2002). Specifically, at each permutation, the size of each cluster was computed, and the maximum cluster size obtained at each permutation was retained to form a null distribution against which the observed cluster sizes were evaluated. Evaluating the observed data against the null distribution of the maximal statistic ensured that the resulting p values were corrected for multiple comparisons (Nichols and Holmes, 2002). This method is analogous to the clusterwise correction procedures commonly implemented in traditional voxelbased fMRI studies (Poline et al., 1997; Bullmore et al., 1999). The advantage of clusterbased approaches is that they are more sensitive to spatially extended genetic signals than thresholding based on independent analysis of regional effects.
Results
Costefficiency of human cortical functional networks
Figure 2A–C illustrates how E_{global}, D_{global}, and E_{global} − D_{global} scaled as a function of connection density, k, at each wavelet scale. The positive E_{global} − D_{global} values observed at each frequency interval indicate that the functional networks were characterized by high efficiency and relatively low connection cost; i.e., they were costefficient. Notably, average CE_{global} was higher for wavelet scales 1 and 2 (0.18–0.35 Hz and 0.09–0.18 Hz, respectively) than for lowerfrequency scales (Fig. 2D).
Genetic influences on global network costefficiency
To consider genetic effects on global network costefficiency, we first examined the MZ and DZ twin correlations for this phenotype at each wavelet scale. For scale 1 (0.18–0.35 Hz), twin correlations were virtually identical (r_{MZ} = 0.44 and r_{DZ} = 0.45), suggesting minimal genetic influence. In contrast, twin correlations at scale 2 (0.09–0.18 Hz) pointed to a pronounced genetic effect, with the MZ correlation being substantially larger than the DZ correlation (r_{MZ} = 0.61 and r_{DZ} = 0.23). Structural equation modeling (Neale et al., 2003) was used to partition the phenotypic variance in network costefficiency into additive genetic (A), shared environmental (C), and unique environmental/error (E) components, represented by the parameters a^{2}, c^{2}, and e^{2}, respectively (see Materials and Methods). Additive genetic effects accounted for 60% of the phenotypic variance in network costefficiency at wavelet scale 2 (a^{2} = 0.60, 95% CI: 0.17–0.83), with unique environmental effects accounting for the remainder (e^{2} = 0.40, 95% CI: 0.19–0.72) and shared environmental influences being negligible (c^{2}=0, 95% CI: 0–0.20). Accordingly, a model incorporating only additive genetic and unique environmental components provided the best fit to the data (−2LL_{AE} = 97.18, χ_{(54)}^{2} = 0, p = 1.0 vs −2LL_{CE} = 99.370, χ_{(54)}^{2} = 2.19, p = 0.02).
As with scale 1, twin correlations at scales 3 (0.04–0.09 Hz) and 4 (0.02–0.04 Hz) indicated a lack of genetic effects, being very low for both MZ (scale 3: r_{MZ} = 0.10; scale 4: r_{MZ} = −0.19) and DZ (scale 3: r_{DZ} = −0.09; scale 4: r_{DZ} = −0.38) twins. The low and negative twin correlations observed at scales 3 and 4 likely reflect the aforementioned measurement error associated with estimating functional connectivity at these lower frequencies. As a result, genetic parameters were not estimable at these scales (see Discussion).
Genetic influences on regional network costefficiency
We further examined genetic influences on the regional costefficiency, CE_{regional}, of each cortical node at wavelet scale 2 (0.09–0.18 Hz). Estimates of r_{DZ} were relatively low across most cortical regions (median = 0.20, interquartile range = 0.18), whereas r_{MZ} values were consistently higher (median = 0.54, interquartile range = 0.13), particularly in regions of medial and lateral prefrontal cortex, where r_{MZ} ≥ 0.6 (Fig. 3).
Genetic models were estimable for 90% of nodes (see Materials and Methods) and indicated that genetic influences were not distributed homogeneously throughout the cortex: regional heritability (a^{2}) estimates ranged from <0.10 to 0.81, with a median of 0.51 (Fig. 3). Regions showing statistically significant genetic effects were found throughout the cortex, in a largely bilateral pattern (Fig. 4). These bilateral clusters were found in posterior cingulate, dorsomedial prefrontal cortex (extending onto the superior frontal gyrus), the posterior middle frontal, postcentral, and lingual gyri, the superior parietal lobule, and the ventroposterior bank of the superior temporal sulcus (extending into adjacent temporal and occipital cortices). Additional leftlateralized clusters were found in ventromedial prefrontal, anterior cingulate, primary visual, and parahippocampal cortices; rightlateralized clusters were found in the precuneus, paracentral lobule, and parietal operculum. No statistically significant regional effects of shared environment (c^{2}) were observed.
Genetic influences on other network properties
To compare genetic effects on cortical costefficiency with other network properties, we also examined twin correlations and genetic parameter estimates for each of the component measures comprising our CE_{global} metric (i.e., E_{global}, D_{global}, and E_{global} − D_{global}) (supplemental Fig. S2, available at www.jneurosci.org as supplemental material), in addition to other commonly studied topological measures: E_{local}, λ, γ, and σ (supplemental Fig. S4, available at www.jneurosci.org as supplemental material). In general, higherfrequency networks (i.e., wavelet scales 1 and 2) showed lower path length and higher clustering, and therefore higher values of the smallworldness parameter, σ (supplemental Fig. S3, available at www.jneurosci.org as supplemental material). In particular, networks at wavelet scale 2 (0.09–0.18 Hz) were associated with the highest smallworldness values for k < 20%. Networks operating at lower frequencies (wavelet scales 3 and 4) were generally associated with higher local efficiency (supplemental Fig. S3, available at www.jneurosci.org as supplemental material).
Consistent with the costefficiency findings, twin correlations were generally low for scales 3 and 4 across all measures, and genetic models were not estimable in most cases. Twin correlations and genetic effects were higher for certain measures in scale 1 and 2 networks, but the results were often dependent on the connection density of the networks. For example, the proportion of variance in E_{local} attributable to additive genetic effects at scale 2 showed >4fold variation for different values of k (supplemental Fig. S4, available at www.jneurosci.org as supplemental material). In some cases, heritability estimates were comparable to those obtained for costefficiency, though examination of nodal properties suggested the estimates were unreliable. For example, for wavelet scale 2, a^{2} = 0.63 for E_{global} at k = 10% (supplemental Fig. S2, available at www.jneurosci.org as supplemental material), but nodal genetic parameters were estimable for only 10% of regions due to low and/or negative twin correlations (see Materials and Methods). This suggests that the global effect was an unreliable and unrepresentative index of nodal characteristics. Together, these findings indicate that costefficiency of cortical functional connectivity was the most reliably heritable measure of all network organizational properties studied, at both global and regional levels.
Discussion
The principles underlying cortical network organization have long been debated, largely due to difficulties associated with adequately modeling the brain's largescale connectivity architecture. Our results indicate that one important principle for brain functional networks involves the optimization of connection costefficiency. Cost efficiency of cortical functional connectivity was under strong genetic control, with 60% of the variance in global network costefficiency attributable to genetic factors. In some regions, >70% of the variance was accounted for by genetic influences. These heritability estimates are comparable to those reported for wholebrain and regional gray matter volume, which range between 0.66 and 0.97 (Thompson et al., 2001; Peper et al., 2007); cognitive abilities, which range between 0.30 and 0.80 (McClearn et al., 1997; Bouchard, 1998; Boomsma et al., 2002); and personality and psychopathology measures, which range between 0.50 and 0.60 (Boomsma et al., 2002). Thus, costefficient organization of cortical functional networks shows heritable influences commensurate with those reported for some of the most widely studied neurobehavioral phenotypes to date.
The heritability of brain network costefficiency
The benefits of enhancing network communication efficiency are intuitive: lower mean path length between nodes results in faster information transmission and lower probability of signal corruption. Accordingly, increased brain network efficiency has been associated with better cognitive performance (Li et al., 2009; van den Heuvel et al., 2009; Zalesky et al., 2011). Adding more connections will always increase efficiency, but the resulting benefits must be balanced with constraints on the brain's energy expenditure (Laughlin and Sejnowski, 2003). In anatomical brain networks, this constraint has been operationalized as a pressure to minimize wiring costs (Chklovskii et al., 2002; Cherniak et al., 2004; Kaiser and Hilgetag, 2006).
Our findings suggest that costefficiency optimization is also an important principle of functional brain organization, consistent with evidence that functional network costefficiency predicts cognitive performance and is affected by aging and disease (Achard and Bullmore, 2007; Bassett et al., 2009). There are substantial technical challenges associated with directly measuring metabolic costs in functional networks. The cost measure we used, operationalized in terms of brain network topology, represents a surrogate marker based on the reasonable assumption that the metabolic cost of each connection scales in proportion to its distance (Laughlin and Sejnowski, 2003). However, it remains an indirect cost estimate, and combined analysis of rsfMRI and diffusionimaging data (e.g., Honey et al., 2009) should refine this measure in future. Such work could also determine whether genetic effects on functional network costefficiency are independent of, or an epiphenomenon of, influences on anatomical network organization, given the correlation between the two (Honey et al., 2009; Zalesky and Fornito, 2009), and the demonstrated genetic effects on white matter architecture (Chiang et al., 2009).
Significant effects were observed throughout the cortex in a largely bilateral pattern. Many of these regions, including left ventromedial prefrontal cortex, right precuneus, and bilateral posterior cingulate, dorsomedial prefrontal and inferior parietal cortices, form part of the socalled defaultmode network. This ensemble of regions shows high metabolic activity and functional connectivity when individuals are not engaged in an explicit task (Raichle et al., 2001; Buckner et al., 2008), and its constituent regions act as major connectivity hubs in functional brain networks (Buckner et al., 2009). Additional regions showing significant genetic influences were found in lateral prefrontal and parietal regions, areas known to activate during attentionally demanding tasks (Smith et al., 2009). Many of these regions are implicated as part of a structural core of anatomical brain network connectivity (Hagmann et al., 2008). Thus, one hypothesis is that these network nodes comprise a genetically determined backbone of structural and functional connectivity that facilitates rapid, efficient interregional communication for low connection cost. The widespread, bilateral distribution of these nodes is consistent with their involvement as critical relay points that integrate information from diverse cortical regions.
Glahn et al. (2010) recently reported significant heritability of functional connectivity within the default mode network, and that connectivity within the network was influenced by shared genetic factors. Notably, heritability estimates for their connectivity phenotype were <0.43. These are considerably lower than those observed in the current study, and those reported in a recent electroencephalographic study of brain network topology, which were generally >0.60 (Smit et al., 2008). This discrepancy may reflect stronger genetic influences on higherorder organizational principles of functional brain networks, such as costefficiency optimization. Alternatively, it may reflect sampling and analytic differences, as Glahn et al. (2010) used an extended pedigree design.
Some topological properties related to our costefficiency measure also showed high heritability. At some connection densities, up to 63% of the variance in global efficiency, E_{global}, was attributable to additive genetic effects. This is to be expected, given that our costefficiency metric, CE_{global}, and the measures used to derive it are related by definition. However, heritability estimates in these component measures were less stable for two reasons. First, they were highly contingent on the connection density at which they were studied. An attractive property of our CE_{global} measure is that it is computed using individually tailored, datadriven thresholds, which likely increased our sensitivity to detect genetic effects (Koten et al., 2009). Second, nodal genetic parameters for these measures could only be estimated in a minority of cases, as most regions were associated with r_{MZ} or r_{DZ} < 0, or r_{DZ} > r_{MZ}. Such findings may reflect noisy estimation of twin correlations in a relatively small sample, but the fact that this trend was greatly reversed when examining network costefficiency (e.g., at scale 2, genetic parameters were estimable for 90% of nodes for CE_{global}, compared with only 10% nodes for E_{global} at k = 10%) indicates that our CE_{global} measure better captured the genetic variance in network topology. Larger samples may be required to appropriately characterize genetic effects on other measures.
Genetic effects across different frequency intervals
Twin correlations were low and genetic models unreliable for most network properties studied at wavelet scales 3 and 4, corresponding to the commonly studied <0.1 Hz range. One explanation for this finding concerns the error associated with estimating time series correlations, which increases at lower frequencies (Whitcher et al., 2000; Achard et al., 2008) (supplemental Fig. S1, available at www.jneurosci.org as supplemental material). This may be a particularly important consideration when estimating very lowfrequency (e.g., <0.1 Hz) properties from limited time series, as acquired during a typical fMRI experiment. Heritability is the proportion of variability in a trait attributable to genetic factors. If measurement error for the trait is large, the proportion of variation attributable to genetic effects will be lower than in a comparable trait measured with less error. This will reduce power for detecting significant genetic influences. In addition, lowerfrequency networks were less costefficient (Fig. 2D) and less “smallworld,” showing both reduced clustering and higher path length (supplemental Fig. S3, available at www.jneurosci.org as supplemental material). These findings are consistent with the lower frequency networks being constructed from noisier functional connectivity estimates. Thus, our findings suggest either that topological properties of lowfrequency networks are not significantly heritable or that more reliable estimation of these properties using longer time series may be required to identify any putative genetic influences.
Genetic effects were strongest and most robustly estimated at wavelet scale 2, corresponding to the 0.09–0.18 Hz interval. Most rsfMRI experiments filter out signals in this frequency range because signal power is typically highest below 0.1 Hz, and noise from physiological sources can contaminate higherfrequency signal fluctuations (Cordes et al., 2001). While our scale 1 measurements may have been affected by such noise, the scale 2 results are unlikely to be contaminated as the 0.09–0.18 Hz band is still lower than typical cardiac (∼1 Hz) and respiratory (∼0.3 Hz) cycles (Lowe et al., 1998).
To control for residual noise, particularly from aliased components of undersampled physiological cycles, our data were corrected using a databased method shown to be as effective as modelbased procedures that incorporate physiological recordings (Behzadi et al., 2007). These considerations, combined with the strong genetic effects observed, suggest that the genetic effects we observed on functional network costefficiency do indeed exert their influence on neuronal dynamics.
The functional significance of 0.09–0.18 Hz BOLD signal oscillations has seldom been investigated, but recent work has found that functional connectivity in this frequency range can distinguish between healthy and patient populations (Garrity et al., 2007; Malinen et al., 2010) and varies in accordance with changing task contexts (Salvador et al., 2008). Collectively, these findings suggest that genetically and behaviorally relevant variance in restingstate BOLD signal fluctuations is manifest at frequencies >0.10 Hz, consistent with evidence that these fluctuations may possess a broader frequency spectrum than previously thought (for a discussion, see Cole et al., 2010).
Conclusions
The strong genetic effects observed here support the hypothesis that functional brain networks are organized according to competitive selection criteria of maximizing communication efficiency and minimizing connection costs. An important avenue of future research will involve understanding how such effects relate to other putative organizational principles for the brain (Simon, 1962), and to genetic influences on anatomical connectivity.
Footnotes

This work was supported by a National Health and Medical Research Council (NHMRC) Project Grant and the Ian Potter Foundation. A.F. was supported by an NHMRC CJ Martin Fellowship (ID: 454797). A.Z. is supported by the Australian Research Council, Discovery Project 0986320. M.Y. and S.J.W. were supported by an NHMRC Clinical Career Development Award. D.S.B. was supported by the National Institutes of Health/Cambridge Graduate Training program. We thank Drs. Sophie Achard and Sarah Medland for helpful advice.
 Correspondence should be addressed to Dr. Alex Fornito, Melbourne Neuropsychiatry Centre, Levels 2 and 3, 161 Barry Street, Carlton South, VIC 3053, Australia. fornitoa{at}unimelb.edu.au