Changes in connectivity profiles define functionally distinct regions in human medial frontal cortex

Methods

Data Acquisition. Diffusion-weighted, blood oxygenation level-dependent (BOLD) fMRI, and T1-weighted images were acquired in nine healthy subjects (ages 24-35 years, five male) on a 1.5 T Siemens Sonata MR scanner with maximum gradient strength of 40 mT·m^-1. All subjects gave informed written consent in accordance with ethical approval from the Oxford Research Ethics Committee.

Diffusion-weighted data were acquired by using echo planar imaging (72 × 2 mm thick axial slices, matrix size 128 × 104, field of view 256 × 208 mm², giving a voxel size of 2 × 2 × 2 mm). The diffusion weighting was isotropically distributed along 60 directions by using a b value of 1,000 s·mm^-2. For each set of diffusion-weighted data, 5 vol with no diffusion weighting were acquired at points throughout the acquisition. Three sets of diffusion-weighted data were acquired for subsequent averaging to improve signal-to-noise ratio. The total scan time for the diffusion-weighted imaging (DWI) protocol was 45 min.

BOLD fMRI data were acquired by using echo planar imaging (20 × 5 mm thick axial slices positioned from the top of the brain, matrix size 128 × 128, field of view 256 × 256 mm², giving a voxel size of 2 × 2 × 5 mm, repetition time = 2.5 s, 341 vol, echo time = 45 ms, and flip angle = 90°).

Subjects practiced the fMRI tasks before entering the scanner. Blocks (30-s duration) of rest (i) alternated with blocks of finger tapping (ii) or serial subtraction (counting backward in threes) (iii) in an ABACACAB cycle repeated 3.5 times. The current task was indicated by the word “rest,” “move,” or “count” displayed on a projection screen at the foot of the scanner bed viewed via a mirror. During move blocks, subjects pressed buttons with the fingers of their right hand in a repeating 1234321 sequence at a frequency of ≈4 Hz. During count blocks, subjects counted covertly backwards in threes from a three-digit reference number displayed on the screen for 2 s before the start of the counting block. To ensure task compliance, at the end of each counting block a red screen instructed subjects to report the number they had reached by pressing buttons with their index figure to indicate tens and middle finger to indicate units (e.g., 47 = index finger button four times and middle finger button seven times). The total scan time for the fMRI protocol was ≈15 min.

A T1-weighted anatomical image was acquired by using a flash sequence (repetition time = 12 ms, echo time = 5.65 ms, and flip angle = 19°, with elliptical sampling of kspace, giving a voxel size of 1 × 1 × 1 mm in 5 min and 5 sec).

DWI Analysis. Diffusion data were corrected for eddy currents and head motion by using affine registration to a reference volume (26). Data from the three acquisitions were averaged to improve the signal-to-noise ratio. Probability distributions on fiber direction were calculated at each voxel by using described methods (27). Probabilistic tractography was then performed from voxels within specified seed masks (12, 27).

A medial frontal cortex mask was defined on the group average T1-weighted image by using the program fslview, which can be accessed at www.fmrib.ox.ac.uk/fsl (Fig. 1). The mask included gray matter on the medial wall and extended from y = -22 to y = 30 [Montreal Neurological Institute (MNI) coordinates] and from a short distance above the cingulate sulcus (as visible on the average brain) to the dorsal surface of the brain. Two single-slice masks were used for initial parcellation: an axial slice (MNI Z = 58) and a sagittal slice (MNI X = -2). These masks were transformed into the space of each subject's diffusion data by using the program flirt (26). For each subject, probabilistic tractography was run from all voxels in this seed mask (12, 27). During tractography, knowledge of location was maintained in three spaces: Seed space (MNI), the space of the original DWI data, and a low-resolution (5 × 5 × 5 mm³) space. All diffusion information used in tractography was sampled on the original diffusion imaging grid. Probabilities of connection from each seed voxel (at 2 × 2 × 2 mm³ resolution) to every other voxel in the low-resolution brain (at 5 × 5 × 5 mm³) were binarized and stored in a matrix, A, of dimensions (number of seed voxels × number of voxels in the rest of the brain). The connectivity profiles were binned at this lower resolution for reasons of data storage. The cross correlation matrix of A was computed, resulting in a symmetric matrix, B, of dimensions (number of seeds × number of seeds) in which the (i,j)th element value is the correlation between the connectivity profile of seed i and the connectivity profile of seed j. The nodes in B were permuted by using a spectral reordering algorithm (Supporting Text, which is published as supporting information on the PNAS web site, ref. 28, and D.J.H., unpublished work) that finds the reordering that minimizes the sum of element values multiplied by the squared distance of that element from the diagonal, hence forcing large values toward the diagonal. If the data contains clusters (representing seed voxels with similar connectivity), then these clusters will be apparent in the reordered matrix and break points between clusters will represent locations where connectivity patterns change. Clusters were identified by eye as groups of elements that were strongly correlated with each other and weakly correlated with the rest of the matrix. Elements that did not clearly belong to a single cluster were left unclassified.

Population maps of resulting clusters in brain space were derived by binarizing clusters corresponding to putative SMA (most posterior cluster) and pre-SMA (most anterior cluster) for each subject and averaging these binarized clusters across subjects so that voxel values in the population maps indicated the proportion of subjects in whom a cluster was present at that point.

Diffusion image analysis was carried out on an Intel 2.4-GHz processor. Approximate processing times were as follows: Fitting local diffusion parameters at each brain voxel, 24 h per subject; probabilistic tractography, 10 s per voxel; and reordering a cross-correlation matrix of 400 × 400, 90 min (for larger matrices, time scales at approximately N³).

fMRI Analysis. fMRI data were analyzed by using programs contained in the fmrib software library, which can be accessed at www.fmrib.ox.ac.uk/fsl. The following preprocessing steps were applied: motion correction by using the program mcflirt, removal of nonbrain structures by using the program BET (29), spatial smoothing by using a Gaussian kernel of full width at half maximum of 3 mm, mean-based intensity normalization of all volumes by the same factor, temporal high-pass filtering by using Gaussian-weighted least-squares fitting with a filter of σ = 57.5 s. Time series statistical analysis was carried out by using the program film with local autocorrelation correction (30). Z (Gaussianized T) statistic images were thresholded by using Gaussian Random Field theory-based maximum height thresholding with a corrected significance threshold of P = 0.01. Registration to standard space was carried out by using the program flirt (26).

Assessing Structure-Function Correspondence. For each subject, voxels that were suprathreshold for either functional task and were contained within the anatomically defined medial frontal volume of interest defined a seed mask for further connectivity analyses. Cross-correlation connectivity matrices for these voxels were derived, reordered, and mapped back onto the brain as before to define putative SMA and pre-SMA for each subject. Centers of gravity of fMRI activation clusters and of connectivity-defined clusters were found and the distances between connectivity- and functionally-defined clusters were compared by using paired Student t tests. To characterize the connections of SMA and pre-SMA, the connectivity distributions for all voxels within each region were found for each subject. These profiles were then averaged across subjects and were mapped onto the average T1-weighted brain in standard space.

Results

For each subject, DWI data were used to perform probabilistic tractography (12, 27) from voxels within medial frontal “seed” masks. Probabilities of connection from each seed voxel to every other brain voxel were binarized and stored in a matrix, A, whose cross correlation matrix, B, was found. Elements in matrix B therefore express the correlation in connectivity profile between medial frontal seed points. The nodes in matrix B were permuted by using a spectral reordering algorithm (Supporting Text, ref. 28, and D.J.H., unpublished work) that forces large values toward the diagonal. If the data contain clusters (representing seed voxels with similar connectivity), then these clusters will be apparent in the reordered matrix and break points between clusters will represent locations where connectivity patterns change. Note that if such structure is not present in the original data then the reordered matrix will not have a clustered organization.

Connectivity-Based Division of Medial Frontal Cortex. We first defined single-slice orthogonal seed masks on the medial frontal cortex in the axial (MNI Z = 58) or sagittal (MNI X = -2) plane on the group average T1-weighted image (Fig. 1). These masks were registered to each subject's DWI data for generation of connectivity matrices. Reordered connectivity cross-correlation matrices contained clearly identifiable clusters in all subjects (Fig. 2, and Figs. 6 and 7, which are published as supporting information on the PNAS web site). Note that such structure will only be apparent in the reordered matrices if there is clustered organization in the data. The reordered matrices were divided into two or three clusters. When these clusters were mapped back onto the brain they corresponded to discrete regions situated along the anterior-posterior axis of the medial frontal cortex (Figs. 2, 6, and 7). The border between the most anterior and most posterior cluster was located close to the vertical line extending from the anterior commissure (Y = 0), suggesting that the regions correspond to SMA and pre-SMA. To test this hypothesis directly, we compared subregions defined on the basis of connectivity to functional activation sites during tasks designed to involve SMA or pre-SMA selectively.

Medial Wall fMRI Activations. We acquired BOLD fMRI data while subjects performed blocks of finger tapping or serial subtraction (counting backwards in threes) alternating with rest. These two functional tasks were selected because previous studies have shown that finger tapping activates the SMA (24), whereas serial subtraction activates the pre-SMA (31, 32) in the superior medial frontal cortex. Both tasks activated the superior medial frontal cortex in all subjects (Fig. 3a, and Fig. 8, which is published as supporting information on the PNAS web site). In some subjects, there was overlap between activated clusters for the two tasks. In all subjects, medial wall activations during finger tapping were more posterior and superior than those during serial subtraction, as anticipated if finger tapping involves the SMA and serial subtraction involves the pre-SMA.

Testing Structure-Function Correspondence. Medial superior frontal voxels activated in either task were entered into a connectivity analysis for each individual subject. For all subjects, the resulting reordered cross-correlation matrices contained clusters of similar connectivity that were defined by an investigator blind to the fMRI results (Fig. 3b). In all nine subjects, two connectivity clusters were identified, and in two of nine subjects, an additional, smaller cluster was found between the other two (see Fig. 8). When the clusters were mapped back onto the brain, they appeared as distinct regions along medial frontal cortex and the anterior and posterior connectivity clusters corresponded closely to the activated SMA and pre-SMA volumes during fMRI (Figs. 3 a and c and 8).

The centers of gravity for superior medial frontal counting-related activations (pre-SMA) colocalized with the centers of the most anterior connectivity-defined clusters, whereas the centers of movement-related activations (SMA) colocalized with the most posterior connectivity-defined clusters (Fig. 4). For all subjects, the center of the most posterior connectivity cluster was closer to functionally defined SMA (median distance = 2.23 mm, and range = 0.42-5.30) compared with pre-SMA (median distance = 8.02 mm, and range = 4.53-13.29) (P = 0.002), whereas for all subjects the center of the anterior connectivity cluster was closer to functionally defined pre-SMA (median distance = 3.07 mm, and range = 1.5-5.42) than SMA (median distance = 9.13 mm, and range = 3.47-13.17); (P = 0.002).

Connections from SMA and pre-SMA Regions. The finding that regions corresponding to SMA and pre-SMA form distinct clusters in the reordered cross-correlation matrices reflects their different connectivity profiles. To characterize the connectivity profile of each area, we mapped the connectivity distributions from all voxels within putative SMA or pre-SMA for each subject onto the average T1-weighted brain template. These distributions were then averaged across all subjects (Fig. 5). As predicted from literature from nonhuman primates (17, 18, 33-37), connections from SMA were found to the corticospinal tract, precentral gyrus, and ventrolateral thalamus (Fig. 5 a and c), whereas connections from pre-SMA were found to the superior frontal gyrus, medial parietal cortex, inferior frontal cortex, and anterior thalamus (Fig. 5 b and c). Consistent with a recent DWI study in humans (21), pathways from SMA connected to more posterior sites in the putamen than those from pre-SMA. More unexpectedly, connections were seen from SMA to orbitofrontal cortex (data not shown) and from pre-SMA to the external capsule/insula (Fig. 5c). Although there is evidence for a weak orbitofrontal connection from SMA in macaque (38), consistent with its presence here in humans, it also is possible that the orbitofrontal and insula connections originated from superior parts of the cingulate sulcus (35, 39) included in the medial frontal seed mask. The borders between SMA/pre-SMA and the cingulate motor areas are difficult to define and the anatomy of this region is highly variable between subjects (9, 40). The inferior border of our medial frontal mask was located a short distance above the cingulate sulcus on the group average anatomical image, corresponding to Z = 50 at its most caudal end, Z = 46 at the level of the vertical line from the anterior commissure, and Z = 38 at its most rostral end. The connection from SMA to orbitofrontal cortex was more commonly seen in subjects in whom putative SMA extended below the level of Z = 48. Similarly, the connection from pre-SMA to external capsule/insula was most commonly seen in subjects in whom pre-SMA extended below the level of Z = 40.

Discussion

By using a generalizable method for discriminating between gray matter regions based on differences in connectivity, we identified a sharp change in connectivity profile along the superior medial frontal cortex. The specific connectivity of the more posterior region to motor and premotor cortex and the corticospinal tract and of the more anterior region to the inferior frontal gyrus, medial parietal, and superior frontal cortex suggested anatomical homology to SMA and pre-SMA in the macaque brain (17, 18, 33-36). There was an excellent correspondence between these regions defined by connectivity and those identified as SMA or pre-SMA by functional criteria. These results therefore directly establish a close relationship between connectional anatomy and function in the superior medial frontal cortex and add further support to the general principle that variations in connectivity reflect functional specialization (10).

The strategy we present here provides an approach to in vivo investigation of human brain organization. Previously, the only structural boundaries easily measurable in vivo were features such as gyri and sulci that do not correspond well to microstructural borders (4, 6). Conventionally defined microstructural borders based on measures such as cyto-, myelo-, and receptor architecture can only be assessed post mortem. Although recent high-resolution MRI studies have shown some promise for detection of prominent myeloarchitectonic boundaries in specific regions of visual cortex (41, 42), it is not clear whether the approach will be more generally useful in regions with less distinct differences in myeloarchitecture. However, changes in connectivity characterize differences between many different brain regions (10, 11, 17) and could therefore provide a basis for defining microstructural boundaries much more generally.

Diffusion imaging data previously has been used to parcellate subcortical gray matter on the basis of remote connectivity (12) or local diffusion properties (43). For example, classification of thalamic voxels according to the cortical target with which they showed the highest probability of connection resulted in clusters that we proposed correspond to thalamic nuclei or nuclear groups (12, 44). The approach presented here represents a fundamentally different strategy for using connectivity information to parcellate gray matter. First, it does not rely on prior knowledge of what constitute meaningful divisions in connectivity targets. The regions in the medial frontal cortex were detected by using changes in connectivity profiles alone without the need for prior knowledge of what those profiles are. Having detected a change in connectivity, we were able in this case to trace the connections of the two defined regions to their major targets to characterize the different spatial distributions of connectivity profiles. This procedure allowed the areas to be related to homologous regions in the macaque brain. However, accurate and complete tracing to all final targets is not essential simply for discrimination of connectionally distinct regions. This finding is particularly important when tracing distributions from regions of cerebral cortex where low starting anisotropy and the presence of fiber crossing and complexity can make longer distance tractography difficult. However, the accuracy of the underlying tractography process will be reflected in the sensitivity of the approach proposed here. The inclusion of more sophisticated tracking techniques that attempt to model and trace through areas of complex fiber architecture would enable the detection of connectional differences that are invisible to the tracking methods used to drive the results in this paper. In addition, localized image artifacts could lead to spurious differences in connectivity profiles, and therefore, improvements in data quality and artifact reduction will further increase confidence and sensitivity.

The ability to assess borders of connectivity-defined cortical regions noninvasively allows for the individual variation and basis of brain structure to be addressed in another way. It already is known that there is substantial variability in both structural (4, 6, 40) and functional (45, 46) brain anatomy between individuals, but a more complete description demands analysis of substantially greater numbers of individuals than can be studied conveniently by using classical histological methods post mortem. However, there are limits to the information that can be provided by diffusion tractography not only due to limitations in current technology such as the difficulty of tracking in the presence of complex fiber architecture (47, 48) but also due to fundamental limitations of diffusion imaging data. For example, it is not possible to differentiate between anterograde and retrograde connections or to determine whether a connection is direct or indirect. Therefore, classical approaches to identifying connections and defining histology remain crucially important for gaining a full understanding of cortical anatomy.

Clusters in the reordered connectivity matrices were identified by eye and were manually defined in the present study. In the case of connectivity data from the superior medial frontal cortex, matrices could be clearly divided in this way into two or three clusters after reordering (see Figs. 6-8). However, it would be desirable to objectively determine both the number of clusters and the location of break points in reordered matrices. This determination will become increasingly important as the approach is applied to larger cortical volumes. There are a number of established techniques for data clustering, some of which have been successfully applied for example to brain data, for example, to objectively define cytoarchitectonic borders in human cortex (49). Future work will be needed to explore the suitability of such techniques for objectively identifying clusters in connectivity data.

Definition of connectional and functional boundaries in the same brain allows direct testing of structure-function relationships. Previous studies (50, 51) have extrapolated from population maps of cortical areas based on cytoarchitecture in one group of subjects to functionally defined regions of cortex in another group to infer relations between structure and function. A general limitation of this strategy is the uncertain interpretation of group correlations arising from interindividual variations in both the structural and functional borders. Determining whether variations in functional anatomy reflect variations in structural anatomy within an individual has previously been limited by the difficulty of measuring both structural and functional anatomy in the same subject. The potentially generalizable methods described here allow the connectional anatomy of the neocortex to be related both to gross anatomical features and to functional activation patterns within individuals and over relatively large populations.

Supplementary Material