Genetic Influences on Cost-Efficient Organization of Human Cortical Functional Networks

Participants.

Thirty-four (18 male) monozygotic (MZ) and twenty-six (12 male) same-sex 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 mail-outs 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₍₂₇₎ = −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 echo-planar imaging volumes providing T2*-weighted blood oxygenation level-dependent (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. T1-weighted 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 rigid-body 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 second-order displacements against the time series of the realigned images, as implemented in CamBa (http://www-bmu.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 echo-planar imaging template supplied with SPM5 (http://www.fil.ion.ucl.ac.uk/spm/software/spm5) using a 12-parameter 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 trade-off between computational speed, spatial resolution, and signal-to-noise ratio (Fornito et al., 2010).

To generate the parcellation, each participant's T1-weighted 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 surface-based parcellations were converted into volumetric representations and projected out to fill the gray matter of the entire cortical ribbon, defined via atlas-based tissue segmentation of the T1-weighted image (Fischl et al., 2002). This volume-based 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% non-zero intensity voxels, resulting in a final template comprising 1041 cortical regions across both hemispheres (Fig. 1A).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Overview of main processing steps involved in generating graph models of cortical functional connectivity. The cortex of each individual's functional volumes was parcellated using a custom template comprising 1041 regions (A). The mean time series of each region was extracted and decomposed into four distinct frequency bands using a wavelet transform (B). These filtered time series were corrected for noise, and pairwise correlations between regional wavelet coefficient series were computed to generate a {1041 × 1041} functional connectivity matrix at each wavelet scale (frequency interval) for each participant (C). These functional connectivity matrices were then thresholded at different connection densities (examples of 5%, 20%, and 40% are presented) and binarized to generate adjacency matrices (D) from which graph models of network topology were constructed (E).

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 motion-corrected, 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/f-type 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 wavelet-filtered regional time series for fluctuations in non-neuronal 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 model-based correction procedures requiring physiological monitoring (Behzadi et al., 2007).

Network characterization.

To estimate the cost-efficiency 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 graph-based 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 wavelet-filtered, noise-corrected 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 ∼3-fold increase in mean estimation error, compared to an ∼1.1-fold 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 low-frequency properties from relatively short time series (Achard et al., 2008; see also Discussion).

Graph-based 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 frequency-specific 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 graph-based 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 cost-efficiency. Four measures were integral to our definition of network cost-efficiency: 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 = (xi−xj)2+(yi−yj)2+(zi−zj)2 and x, y, and z define the coordinates of each region's centroid in stereotactic space. Similarly, the regional connection cost was calculated as follows: Each of the connection distance measures above was normalized by its corresponding value in a completely connected graph (i.e., when k = 100%) to scale it between 0 and 1.

Previous studies taking the difference between efficiency and binary connection cost as an index of network cost-efficiency 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 cost-efficiency (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 cost-efficient 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 cost-efficiency as follows: Regional cost-efficiency 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 cost-efficient network, and ensure that when all else was equal, networks with fewer connections were considered to be more cost-efficient. 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 data-driven, individually tailored approach to threshold selection because each participant's network is thresholded based on his or her unique maximally cost-efficient 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 cost-efficiency 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 diffusion-weighted 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 cost-efficiency, 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 cost-efficiency 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 small-worldness (σ). 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). Small-worldness was assessed using the scalar summary γ/λ=σ (Humphries et al., 2006). In small-world networks, λ ≈ 1 and γ > 1, so a nonrandom, small-world 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², c², and e², respectively) reflect the proportion of phenotypic variation attributable to their effects. Model path coefficients were estimated using maximum-likelihood 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 log-likelihood, which follows a χ² 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_i1 and T_i2 denote the first and second members of twin pair i. If under random reassignment T_i1 became paired with T_j1, then T_i2 would be paired with T_j2. 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. Ninety-five 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 permutation-based, 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⁶ permuted datasets. At each permutation, clusters of regions showing genetic effects were formed by merging adjacent nodes surviving a primary cluster-forming 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 voxel-based fMRI studies (Poline et al., 1997; Bullmore et al., 1999). The advantage of cluster-based approaches is that they are more sensitive to spatially extended genetic signals than thresholding based on independent analysis of regional effects.