Abstract

We combined quantitative relaxation rate (R1= 1/T1) mapping—to measure local myelination—with fMRI-based retinotopy. Gray–white and pial surfaces were reconstructed and used to sample R1 at different cortical depths. Like myelination, R1 decreased from deeper to superficial layers. R1 decreased passing from V1 and MT, to immediately surrounding areas, then to the angular gyrus. High R1 was correlated across the cortex with convex local curvature so the data was first “de-curved”. By overlaying R1 and retinotopic maps, we found that many visual area borders were associated with significant R1 increases including V1, V3A, MT, V6, V6A, V8/VO1, FST, and VIP. Surprisingly, retinotopic MT occupied only the posterior portion of an oval-shaped lateral occipital R1 maximum. R1 maps were reproducible within individuals and comparable between subjects without intensity normalization, enabling multi-center studies of development, aging, and disease progression, and structure/function mapping in other modalities.

Introduction

Cortical areas are best defined and recognized on the basis of multiple converging techniques (Allman and Kaas 1971). Five measures developed during invasive experiments on animals include 1) receptotopic organization (e.g. retinotopy), 2) architectonic features, 3) connection patterns, 4) neurophysiological properties, and 5) effects of localized lesions. These measures have each been adapted to human brains. Studies on the effects of brain damage in humans have a long pedigree, recently supplemented by transcranial magnetic stimulation studies in normal subjects. Cortical-surface-based functional magnetic resonance imaging retinotopy (Engel et al. 1994; Sereno et al. 1995; DeYoe et al. 1996) is now established and often combined with functional studies, which constitute the overwhelming majority of in vivo human neuroimaging experiments. Connection patterns are beginning to be addressed in vivo with diffusion-imaging-based fiber-tracking methods. But despite the historical precedence of postmortem human architectonic studies, in vivo human architectonics has perhaps been the least well-adapted of the 5 measures.

Postmortem cortical parcellation remains an difficult problem a century after the publication of Brodmann's classic cytoarchitectonic maps of the cortex in humans and other species. Related early work using cytoarchitectonic and myeloarchitectonic methods soon resulted in a plethora of human cortical area maps (reviews: Zilles and Amunts 2010; Geyer et al. 2011). But human cortical parcellation then went into decline, especially after the blistering critique of Lashley and Clark (1946). Bailey and von Bonin (1951), for example, divided the human neocortex into only 6 areas, even blurring the sharp border between V1 and V2 originally recognized by Meynert (1867).

Despite subsequent advances in mapping cortical areas in non-human primates (Merzenich and Kaas 1980), specifically human cortical parcellation languished (though see Braak (1980) on “pigmento”-architectonics). The growth of human neuroimaging in the 1980s had the effect of elevating Brodmann's map—after a rough translation from Brodmann's 2D summary images to a 3D atlas—to a world standard. Despite its acknowledged shortcomings (e.g. the middle temporal area (MT/V5) did not appear), the moderate resolution of early volume-averaged neuroimaging data did not initially demand better. Two decades later, as modern postmortem human cortical architectonics augmented by immunohistochemistry began to appear, the situation has improved considerably (e.g. Malikovic et al. 2007; Caspers et al. 2012).

However, a persistent stumbling block to combining structural and functional data has been to integrate postmortem parcellations based on high-resolution microscopic measures with retinotopic, connectional, and functional data from in vivo brains. Variability between postmortem and in vivo brains—whether at the single subject or group level—introduces irreducible uncontrolled systematic variation.

Most previous attempts at in vivo architectonics have only been able to consider few localized regions of the cortex (Clark et al. 1992; Walters et al. 2003; Eickhoff et al. 2005; Sigalovsky et al. 2006; Geyer et al. 2011). One reason is signal-to-noise constraints due to in vivo scan time limits. But a more serious problem is the difficulty of reliably detecting small signal differences (1–4%) between different parts of non-quantitative T1-weighted (T1w) images that have large artifactual variations in brightness and contrast due to uncorrected inhomogeneities in the local radiofrequency transmit and receive fields. Widely applied post hoc histogram-based normalization methods (e.g. Dale et al. 1999) and T1w/T2weighted (T2w) ratio methods (Glasser and Van Essen 2011) remain sensitive to the effects of uncorrected B1 “transmit” inhomogeneities (flip-angle variation) on brightness and contrast, which affect both within- and between-individual comparisons. Ratio methods must also contend with vessel artifacts and local spatial distortion that differ between the 2 scan types, and unfavorable error propagation inherent in ratio estimation. Finally, differences in local cortical curvature (convex vs. concave) systematically affect the laminar and myeloarchitecture of the cortex, even within a cortical area (Fig. 1; Smart and McSherry 1986; Annese et al. 2004).

Figure 1.

Relaxation rate (R1) as function of cortical depth, area, and curvature. (A) Overall cortical average R1 (equivalent T1 on right y-axis) decreases monotonically from gray–white boundary (depth fraction 0.0) to a slight plateau at middle depths (0.3–0.6, beginning where second derivative in lower graph crosses zero), and then resumes its decrease in superficial layers (0.7–1.0; error bars: ±1 standard error of the mean over subjects). (B) Cross-ROI differences in average cortical R1 shown as line for 8 depths (from 0.1 near white matter to 0.9 near pia); y-axes, error bars as in (A); ROIs: Angular (angular gyrus), angular-fs (FreeSurfer angular gyrus label), MT low and MT high [non-overlapping low- and high-probability MT labels (Malikovic et al. 2007; Fischl et al. 2008), V1-fs (FreeSurfer V1)]. All matched-paired t-tests on hypothesized differences significant (P < 0.05) except where marked “m” (P < 0.1), “=” (no significant difference), or “-” (difference opposite prediction). (C) Vertex-wise correlation (adjusted R2) of R1 and curvature as function of depth (error bars as before over subjects). Scatter plot inset at right is from a single subject at depth 0.5. For comparison, left inset shows myelin-stained section of human cortex with reduced myelination in concave areas (compare 2 white arrows, from Annese et al. 2004).

Figure 1.

Relaxation rate (R1) as function of cortical depth, area, and curvature. (A) Overall cortical average R1 (equivalent T1 on right y-axis) decreases monotonically from gray–white boundary (depth fraction 0.0) to a slight plateau at middle depths (0.3–0.6, beginning where second derivative in lower graph crosses zero), and then resumes its decrease in superficial layers (0.7–1.0; error bars: ±1 standard error of the mean over subjects). (B) Cross-ROI differences in average cortical R1 shown as line for 8 depths (from 0.1 near white matter to 0.9 near pia); y-axes, error bars as in (A); ROIs: Angular (angular gyrus), angular-fs (FreeSurfer angular gyrus label), MT low and MT high [non-overlapping low- and high-probability MT labels (Malikovic et al. 2007; Fischl et al. 2008), V1-fs (FreeSurfer V1)]. All matched-paired t-tests on hypothesized differences significant (P < 0.05) except where marked “m” (P < 0.1), “=” (no significant difference), or “-” (difference opposite prediction). (C) Vertex-wise correlation (adjusted R2) of R1 and curvature as function of depth (error bars as before over subjects). Scatter plot inset at right is from a single subject at depth 0.5. For comparison, left inset shows myelin-stained section of human cortex with reduced myelination in concave areas (compare 2 white arrows, from Annese et al. 2004).

Here, we present a detailed analysis that combines functional and structural data. To overcome the inherent limitations of non-quantitative structural imaging listed above, quantitative T1 maps were acquired for all subjects, providing a measure of local myeloarchitecture (Schmierer et al. 2004; Draganski et al. 2011). The quantification allows for the direct comparison of datasets acquired on different subjects at different times and different MRI sites. The T1 maps were then combined with retinotopic mapping across multiple sessions in the same group of subjects. High-resolution cortical-surface reconstructions of the gray–white matter and pial surfaces were made from the structural scans and used to analyze, visualize, and fuse the T1 and retinotopy data. Since the longitudinal relaxation rate, R1(=1/T1), positively correlates with the level of myelination and the brightness of T1w scans as commonly used for morphometry, we henceforth refer exclusively to it.

Materials and Methods

Six subjects (ages 22–55, 3 female) with normal or corrected-to-normal vision participated in all parts of the study. Experimental protocols were approved by local ethics committees, and participants gave informed and signed written consent.

Structural Imaging

Structural images were acquired on a whole-body Tim Trio system (nominal field strength 3T, actual 2.89T, Siemens Healthcare, Erlangen, Germany), with body transmit and 32-channel receive head coil at the Wellcome Trust Centre for Neuroimaging. Proton density-weighted (PDw) and T1w images were acquired using an in-house multi-echo 3D FLASH pulse sequence (Weiskopf et al. 2011; voxel size: 0.8 × 0.8 × 0.81 mm3, field of view (FOV) = 256 × 216 × 194 mm, matrix = 320 × 270 × 240, repetition time (TR) = 23.7 ms, excitation flip angle: 6° (PDw) or 28° (T1w)). Acquisition was speeded up by 2× GRAPPA parallel imaging in the phase encoding and 6/8 Partial Fourier in the partition direction. To improve image quality, 4 gradient echoes were acquired (echo delay times (TE) = 2.2, 4.75, 7.3, 9.85 ms) after each excitation pulse. Each session consisted of four 10 min 31s acquisitions (2 PDw and 2 T1w) and 2 shorter scans to estimate field inhomogeneities (see below). Quantitative R1 maps were estimated from the PDw and T1w images according to the formalism developed by Helms et al. (2008) including a correction for imperfect RF spoiling (Preibisch and Deichmann 2009). Recent applications of this method have demonstrated its robustness (Helms et al. 2008, 2009; Draganski et al. 2011). To correct for the effect of radio frequency (RF) transmit inhomogeneities on R1 maps, maps of the transmit field B1+ were acquired using a 3D echo-planar imaging (EPI) spin-echo (SE)/stimulated echo (STE) method (Lutti et al. 2010, 2012; FOV = 256 × 192 × 192 mm3, matrix = 64 × 48 × 48, TESE/TESTE = 39.38/72.62 ms, TR = 500 ms, acquisition time 3 min 48 s), which was corrected for off-resonance effects using a B0 fieldmap. For further details, see Supplementary Information.

Functional Imaging and Analysis

Standard T2*-weighted EPI scans (3.2-mm isotropic resolution, 128 or 256 volumes each) and a T1w alignment scan with the same orientation and slice block center were acquired on a 1.5 whole-body Tim Avanto system (Siemens Healthcare) at the Birkbeck/UCL Centre for NeuroImaging, with body transmit and 32-channel receive head coil (see Supplementary Information for details). Five of 6 subjects completed at least 8 retinotopy scans (4 scans in 1 subject) and 4 ipsilateral mapping scans (2560 volumes per subject) across 3–5 additional sessions. See Supplementary Information for details of surface-based retinotopy analysis.

In-house OpenGL/Xlib software drove a video front-projection direct-view system that stimulated the visual field to an eccentricity of at least 57° of visual angle in all directions from central fixation. Balanced (clockwise/counterclockwise, in/out) phase-encoded polar angle (wedge) and eccentricity (ring) stimuli contained medium-luminance-contrast colored checkerboards, optical flow fields, and simultaneous monitor-for-upside-face-among-right-side-up and monitor-for-number-among-letters tasks to maintain continuous peripheral attention (eight 64 second cycles). Ipsilateral field mapping stimuli used low contrast moving versus stationary gratings avoiding the center-of-gaze and the vertical meridian (see Supplementary Information for details).

Cortical-Surface Reconstruction and Sampling of R1 Values Within Cortical Ribbon

Cortical surfaces were reconstructed with FreeSurfer (v5.0.0) from the aligned (AFNI 3dAllineate, hand-inspected) average of the 2 high-resolution T1w scans. We initially used quantitative R1 scans but experienced localized segmentation failures because some boundaries between the pial surface, the CSF, and the skull have different contrast than what is currently assumed by FreeSurfer algorithms (see Discussion).

R1 data sets were sampled along the normal to each gray–white matter surface vertex in steps of 10% of cortical thickness (thickness estimated in FreeSurfer; Fischl and Dale 2000) and then smoothed tangentially at each depth with a 4-mm full width half maximum (FWHM) 2D kernel. The FreeSurfer estimate of local curvature (smoothed with a 1-mm FWHM 2D kernel) was used as a linear predictor of R1 values at each depth. Vertex-wise residuals from this regression were used as “de-curved” estimates of R1 values whose units are directly comparable to raw de-meaned R1 values. See Supplementary Information for details of R1 region of interest (ROI) choice.

Results

R1 as a Function of Cortical Depth

Postmortem studies of cortical myelination demonstrated a consistent pattern of myelination versus depth across most cortical areas. Myelination is highest immediately above the white matter, there are 2 bands of high myelination in deeper layers (inner/deeper and outer/shallower stria of Baillarger), and finally, there is a stepwise reduction in myelination in superficial layers. To investigate whether this overall pattern was apparent in our R1 data, we sampled the R1 values for each surface vertex at different fractional cortical depths, and then averaged these profiles across all vertices.

Average R1 values (Fig. 1A) show an extremely consistent decrease from deep to superficial layers (small to large depth fractions, z). A moderate plateau in deeper cortical layers begins around depth fraction 0.3, where the second spatial derivative with respect to z (∂2R1/∂z2), shown immediately below, crosses zero (an inflection point). The pattern of R1 with depth closely resembles the overall profile of myelination seen in postmortem histology. It is considerably smoother, not distinguishing the 2 stria of Baillarger, which is a reflection of much coarser sampling units of in vivo human MRI (0.8 mm3 here) compared with histology (∼0.01 mm3), as well as the averaging of R1 across the entire cortex of both hemispheres.

R1 Variation over Probabilistically Defined Cortical ROIs

For an initial verification of regional differences in R1 as a measure of myelination, we defined 5 ROIs with large expected myelination differences. The first was the FreeSurfer probabilistic V1 label (V1-fs). The second and third were MT high and MT low (non-overlapping high and low probability of cross-subject overlap regions). Finally, we included 2 angular gyrus labels—a smaller one, angular, defined on our subjects (non-visually responsive region just superior to MT) and a larger one, angular-fs, from the FreeSurfer parcellation.

Based on the histological literature, we predicted 1) that R1 in V1 and the MT ROIs should exceed R1 in the angular gyrus ROIs at all cortical depths, and that V1 should have the highest R1, 2) R1 in high-probability MT should exceed R1 in surrounding lower probability MT in deeper cortical layers (depth fractions 0.1–0.6). The systematic R1 curves across areas at different depths illustrated in Figure 1B confirmed these predictions. All matched-paired t-tests were significant at P < 0.05 except where noted, with “m” indicating marginal (0.05 < P < 0.10) differences, “=” indicating no significant difference, and “-” indicating difference opposite to predicted direction. ROI-based analyses showed that at the level of macroscopic landmarks, our quantitative estimate of R1 parallels known regional differences in myelination.

Relationship of R1 to Local Cortical Curvature

The human cortex has deep sulci, but also a complex pattern of concavity and convexity, including many convex regions buried within major sulci. Postmortem studies in humans show that myeloarchitecture varies significantly with local cortical convexity (Annese et al. 2004); more convex regions are thicker and more myelinated, especially in middle and upper layers. We therefore examined the relationship between R1 and local curvature as a function of cortical depth. Figure 1C shows that there is a moderately strong correlation between local curvature and R1 at middle cortical depth fractions (up to a maximum average adjusted R2= 0.14) that falls off as one approaches the white matter and the pial surface. The bottom right inset scatterplot includes a least-squares fit line for R1 against vertex-wise curvature; vertices in convex regions of the cortex (often but not exclusively on gyri) have higher R1 than vertices in more concave regions. Since R1 varies more at the gray–white matter border and at the pial border than it does within the cortex itself, small errors in tissue segmentation during surface reconstruction can have a large effect on laminar R1 estimates near those boundaries. The fact that R1/curvature correlation is strongest at middle sampling depths suggests that segmentation errors are not responsible for our result.

Because local curvature is also related to local cortical thickness (average vertex-wise adjusted correlation between curvature and thickness: Adjusted R2 = 0.081), local cortical thickness might explain some of the variation in R1. However, the strong correlation between convexity and R1 survives the addition of vertex-wise thickness as a covariate. To minimize the effects of local curvature in our subsequent analyses of R1 variation across the cortical sheet, we therefore initially regressed R1 against curvature, and then used the “de-curved” residual values as our estimates of regional differences in cortical myelination.

Test–Retest Consistency of R1 Across the Cortex

The top 3 rows of Figure 2 show medial and lateral views of the cortical-surface-based cross-subject average (Fischl et al. 1999) of R1 from 6 subjects (2 R1 scans each), displayed on 1 of the 6 subjects' surfaces. The top row (gray–white matter surface), second row (pial surface), and third through sixth rows (inflated surface) show same magnification views of one hemisphere. The heat scale overlay in all rows but row 3 shows “de-curved” R1 values sampled from voxels lying at a cortical depth fraction of 0.5, while row 3 is not “de-curved” (but still de-meaned). The effects of “de-curving” can be seen by comparing rows 3 and 4 (e.g. the thin line posterior to MT disappears). Though regional variations in R1 are small (<0.04s−1), they are remarkably consistent over repeated measurements (bottom 2 rows). R1 data were masked on the midline and in the anterior insula—in the first case because single hemisphere surfaces cut into the corpus callosum and thalamus, and in second case because of the difficulty of accurately reconstructing the gray–white matter surface superficial to the claustrum. Small regions of high cortical R1 in posterior cingulate and anterior insula may thus have been omitted here.

Figure 2.

Cross-subject surface average relaxation rate (R1) test/retest. Top left inset redrawn from Flechsig (1920) shows early myelinating areas in gray. Rows 1–6 show spherical morph cross-subject average ΔR1 (difference above mean) projected back onto the gray–white, pial, and inflated surfaces of 1 subject. R1 values in each subject were sampled at a point halfway between the individual's gray–white matter and pial surfaces. Row 3 shows the data before removing variance due to local cortical curvature; that variance is removed in all other rows. Rows 5 and 6 show excellent scan/rescan reliability of average. Lower middle inset shows same magnification lateral occipital cutout with individual subjects' R1 maxima (yellow contours), superimposed on MT postmortem probability map in purple (Malikovic et al. 2007; Fischl et al. 2008), showing that the individual and average R1 maxima differ subtly from the postmortem map. See text for midline and anterior insula mask.

Figure 2.

Cross-subject surface average relaxation rate (R1) test/retest. Top left inset redrawn from Flechsig (1920) shows early myelinating areas in gray. Rows 1–6 show spherical morph cross-subject average ΔR1 (difference above mean) projected back onto the gray–white, pial, and inflated surfaces of 1 subject. R1 values in each subject were sampled at a point halfway between the individual's gray–white matter and pial surfaces. Row 3 shows the data before removing variance due to local cortical curvature; that variance is removed in all other rows. Rows 5 and 6 show excellent scan/rescan reliability of average. Lower middle inset shows same magnification lateral occipital cutout with individual subjects' R1 maxima (yellow contours), superimposed on MT postmortem probability map in purple (Malikovic et al. 2007; Fischl et al. 2008), showing that the individual and average R1 maxima differ subtly from the postmortem map. See text for midline and anterior insula mask.

These data show a strong resemblance to Flechsig's (Flechsig 1920) survey of perinatal infant myelogenesis, re-rendered in the top left inset of Figure 2. Flechsig found that early myelinating cortical regions (dark shading) tended to be the most densely stained for myelin in the adult cortex. These regions are clearly visible in our R1 data and include primary motor (M-I) and somatosensory (S-I) areas, primary and secondary auditory areas, as well as a number of visual areas described in more detail in the next section.

Relation of R1 to Retinotopic Maps

Visual areas were mapped in the same set of 6 subjects using a wide-field direct-view front-projection system. In 2–3 additional scan sessions for each subject, we mapped polar angle, eccentricity, and ipsilateral low-contrast visual field representations (Huk et al. 2002). We then performed cortical-surface-based cross-subject averages of complex-valued (amplitude and phase) data (Sereno and Huang 2006; Hagler et al. 2007) and calculated visual field sign from the polar angle and eccentricity averages to identify borders in early areas (Sereno et al. 1995). A detailed pattern of correspondences between R1 and retinotopic borders emerged when the 2 data sets were fused on the same surface.

The retinotopic and R1 averages from both hemispheres in Figures 3 and 4 are illustrated in identical lateral, posterior, medial, and inferior poses (inflated surface) to aid comparison and are described together. Iso-R1 borders (ΔR1= 0.020s−1) around R1 maxima were traced and superimposed on the polar angle data using thin white dashes to aid comparisons (see Supplementary Information for a movie that smoothly varies view and color scale contrast, more clearly illustrating the systematic fine level differences in the data shown in Fig. 4). Retinotopic borders (vertical meridian → circles, horizontal meridian → thick black dashes) are illustrated on both figures as is the posterior boundary of ipsilateral responses to low contrast moving gratings (thick yellow dashes in top row).

Figure 3.

Cross-subject surface average retinotopic maps. Spherical morph average polar angle maps are projected back to the inflated right and left hemispheres of one subject, and shown in lateral (top), posterior (middle-top), medial (middle-bottom), and inferior (bottom) views. The posterior boundary of ipsilateral visual responses is marked by a thick yellow dashed line. Thin dotted lines indicate the boundaries of regions with high quantitative R1 values traced from Figure 4. Vertical and horizontal meridians traced from field sign calculations (not shown) are shown as lines of small circles and thick black dashes. Generalized visual and multisensory area names (V1, V2, V3, VP/V3v, V6, V6A, V8/VO1, V3A, MT/V5, FST, LIP (multiple), VIP (multiple), PrCu [pre-cuneus visual area]) were drawn judiciously from the existing conflicting literature.

Figure 3.

Cross-subject surface average retinotopic maps. Spherical morph average polar angle maps are projected back to the inflated right and left hemispheres of one subject, and shown in lateral (top), posterior (middle-top), medial (middle-bottom), and inferior (bottom) views. The posterior boundary of ipsilateral visual responses is marked by a thick yellow dashed line. Thin dotted lines indicate the boundaries of regions with high quantitative R1 values traced from Figure 4. Vertical and horizontal meridians traced from field sign calculations (not shown) are shown as lines of small circles and thick black dashes. Generalized visual and multisensory area names (V1, V2, V3, VP/V3v, V6, V6A, V8/VO1, V3A, MT/V5, FST, LIP (multiple), VIP (multiple), PrCu [pre-cuneus visual area]) were drawn judiciously from the existing conflicting literature.

Figure 4.

Cross-subject surface average quantitative R1 (=1/T1) maps. Spherical morph average maps of quantitative R1 values sampled at 50% of cortical thickness are projected back to the same subject as in Figure 3, and posed and annotated identically. Quantitative R1 is illustrated as variation from the mean (ΔR1) after removing variance due to local cortical curvature. A histogram of ΔR1 values from 3 right hemisphere visual areas (V6, MT, and VIP, defined by retinotopy) and 2 right hemisphere non-visual areas (inferior pre-cuneus “default mode” area and non-visual lateral occipital cortex) is shown in the lower middle inset. The non-visual regions of interest are shown as dashed lines in matching blue and green (upper left and middle left). The maxima shown are 3–4% higher than the average R1.

Figure 4.

Cross-subject surface average quantitative R1 (=1/T1) maps. Spherical morph average maps of quantitative R1 values sampled at 50% of cortical thickness are projected back to the same subject as in Figure 3, and posed and annotated identically. Quantitative R1 is illustrated as variation from the mean (ΔR1) after removing variance due to local cortical curvature. A histogram of ΔR1 values from 3 right hemisphere visual areas (V6, MT, and VIP, defined by retinotopy) and 2 right hemisphere non-visual areas (inferior pre-cuneus “default mode” area and non-visual lateral occipital cortex) is shown in the lower middle inset. The non-visual regions of interest are shown as dashed lines in matching blue and green (upper left and middle left). The maxima shown are 3–4% higher than the average R1.

V1 (see medial views) was characterized by high R1R1 ∼0.031s−1), and its retinotopically defined borders in both hemispheres corresponded almost exactly with a sharp drop in R1. Superior to V1, V6 was bilaterally identified by a characteristic upper-to-lower field transition (moving superiorly and posteriorly) on the posterior bank of the parieto-occipital sulcus (Pitzalis et al. 2006) with R1 values similar to V1 (ΔR1 ∼0.035s−1) with a shallow R1 valley between it and V1 (see contrast ramp movie in Supplementary Information). Interestingly, Flechsig (1920) identified an early myelinating medial area in an almost identical location (Fig. 2, top inset, medial view). There was an additional shallower maximum of myelination (ΔR1 ∼0.012s−1) just anterior and superior to V6, containing mostly lower visual fields that we have tentatively labeled V6A (see Gamberini et al. 2011, for finer subdivisions).

In lateral occipital cortex, there was a prominent oval of high R1R1 ∼0.029s−1). We initially expected this might correspond to MT/V5. However, both polar angle retinotopy (reversal at an upper visual field vertical meridian) and the posterior bound of ipsilateral responses were situated near the middle of the oval. This region was located within a sulcus (Supplementary Fig. 2, showing individual subject data on each folded surface; the average was also resampled to each subject and shown as a green outline). We labeled the posterior part MT. To avoid name proliferation, we labeled the anterior and inferior part of this area, fundus of the superior temporal sulcus area (FST), on the basis of a similar retinotopic layout in non-human primates [e.g. macaque monkey, owl monkeys; see also Kolster et al. (2010) and Amano et al. (2009) who have subdivided this region into a larger number of areas with somewhat different names].

In dorsolateral occipital cortex near the midline, we identified a patch of high R1R1 ∼0.028s−1) that corresponded well with V3A (defined by the first appearance of upper field here moving anteriorly) in the left hemisphere and roughly so in the right. Inferiorly, beginning with the center-of-gaze representation of V1 (see inferior view), there is a large anteriorly extending region of high R1 values (ΔR1 ∼0.031s−1) that also covers the adjoining center-of-gaze representations V2, V3/VP, and hV4 in both hemispheres. This maximum extends across parts of at least 4 different visual areas and may be artifactual (see Discussion). Still on the inferior surface, moving medially along the horizontal meridian of hV4 (Wandell et al. 2007), an anterior right-angle bend is visible in a region originally identified as V8 by Hadjikhani et al. (1998). This precisely corresponds to a small maximum of R1R1 ∼0.021s−1) in both hemispheres.

Back on the lateral surface, there is a prominent maximum of R1R1 ∼0.020s−1) in both hemispheres at the expected location of the ventral intraparietal area (VIP) (Sereno and Huang 2006). This is connected to the large high R1 strip in somatosensory and motor cortex via a small isthmus. In our polar angle mapping, this region contained lower visual field and horizontal meridian responses. The lower field bias may partly reflect head-centered remapping of receptive fields as a result of the lowered gaze position common with direct-view mapping.

Just posterior to VIP, best visible in the posterior view, there is an elongated region of moderately high R1 values (ΔR1 ∼0.013s−1) in both hemispheres that extends posteriorly through the region originally identified as the human intraparietal area (LIP) (Sereno et al. 2001) eventually joining up with the prominent V3A maximum. Subsequently, this region has been subdivided into a number of LIP's (or IPS's) and contains several somewhat variable visual maps (e.g. Swisher et al. 2007; Silver and Kastner 2009). Finally, in frontal cortex, there are shallow R1 maxima in the frontal eye fields, in an area identified as the polysensory zone (Graziano and Gandhi 2000; Sereno and Huang 2006; Huang and Sereno 2007), which appears as an extension off the motor strip, and there is a shallow R1 maximum in the retinotopic part of the dorsolateral prefrontal cortex (Hagler and Sereno 2006).

Discussion

Cortical-surface-based analysis of quantitative longitudinal relaxation rate, R1 (=1/T1), showed reliable small (1–4%) regional differences in cortical gray matter likely due to differences in myelination, as verified by test–retest and by statistical analysis in selected regions defined in previous studies (postmortem-defined V1 and MT, angular gyrus). We compared this in vivo architectonic measure with retinotopic mapping data obtained in separate sessions from the same subjects. This revealed that a number of boundaries visible in the R1 maps correspond to previously recognized retinotopic borders, demonstrating the feasibility of using R1 maps to help with cortical parcellation. Though not all adjoining areas have distinguishable R1 values, a number can now be identified structurally without elaborate functional scans.

The cortex-wide average of R1 as a function of laminar depth resembled postmortem myeloarchitectonic measures, systematically decreasing from the gray–white matter border to the pial surface with a moderate plateau in middle layers. These results are consistent with findings from localized high-resolution imaging and postmortem imaging (Walters et al. 2003; Eickhoff et al. 2005). The R1 ranking of different ROIs—from most to least myelinated—is also remarkably consistent over sampling depths.

R1 was significantly related to local curvature at intermediate depths—with convex parts of the cortex more strongly myelinated—quantitatively confirming observations in myelin-stained sections (Smart and McSherry 1986; Annese et al. 2004). Local curvature predicted R1 better than a measure of whether a cortical region was on a sulcus or gyrus (the primary measure used for cross-subject surface alignment; Fischl et al. 1999). The strength and generality of this finding across the cortex suggests that there is a fundamental asymmetry between concave and convex bends in the cortex having to do with differences in the way deeper and more superficial layers can be deformed (Xu et al. 2010).

At a macroscopic level, the pattern of R1 across the cortex showed a remarkable resemblance to postmortem studies of the order of myelogenesis (Flechsig 1920). By comparing these R1 maps with retinotopic maps, we found that the borders of a number of visual areas defined in the same subjects—including V1, V3A, V6, V6A, V8, VIP, retinotopic dorsolateral prefrontal cortex, and the frontal eye fields—were marked by abrupt increases in R1 (0.02–0.04s−1). Perhaps our most surprising finding concerned the heavily myelinated oval in the lateral occipital cortex. Previously identified as MT/V5 in a number of studies, retinotopic mapping (polar angle, posterior boundary of ipsilateral responses) showed that MT proper only accounted for the posterior 1/3 to 1/2 of that oval. This suggests that previous in vivo and postmortem studies in humans may have overestimated the extent of MT. Given that there are several areas with above average myelination anterior and inferior to MT in monkeys (MST and FST—e.g. Bock et al. 2009), this observation is not unprecedented; however, those other areas may be relatively larger or contain additional subdivisions in humans. Somatomotor and auditory cortex were also prominently recognizable (treated elsewhere).

Potential Applications

In contrast to T1w imaging, quantitative R1 imaging directly estimates a basic physical property, the longitudinal relaxation rate of protons in a given tissue. Thus, the resulting maps are directly comparable without intensity normalization across individuals, scanners, and time—which especially lends itself to studying cortical areal differences, development, health versus disease, and disease progression. In this respect, quantitative neuroimaging is superior to standard postmortem myelin stains, which, despite their much higher resolution, are more subject to the uncontrollable vagaries of silver impregnation.

As an example, R1 values in middle cortical layers for the high-probability MT label differ significantly from immediately adjacent regions, with an average R1 value of ∼0.62s−1 at a fractional cortical depth 0.5 (Fig. 1B), and an average de-curved ΔR1 of 0.023s−1. To investigate whether this highly myelinated region could be identified on a single subject basis, we thresholded each subject's surface-smoothed (4-mm FWHM) de-curved R1 map at the average ΔR1 value of 0.023s−1 and searched for a disconnected supra-threshold R1 patch within a lateral occipital MT+ search space based on previous human studies (Annese et al. 2005). Such patches were identified in all 6 subjects, and in 11 of 12 hemispheres (Supplementary Fig. 1).

Methodological Issues and Prospects

Automated in vivo mapping of architectonics in single subjects requires detecting local changes in R1 on the order of 1% at a high isotropic resolution of ∼800um, which is made possible by advanced quantitative methods. In particular, we employed high-quality mapping of the B1+ transmit field (Lutti et al. 2010, 2012) and highly sensitive multi-echo 3D FLASH R1 mapping to achieve maximal resolution, accuracy, and precision in the shortest time possible (Weiskopf et al. 2011). We optimized the TR, flip angle and RF spoiling phase increment for highest accuracy and precision (Helms et al. 2011) and applied corrections for imperfect spoiling of magnetization coherence pathways (Preibisch and Deichmann 2009). No post hoc brightness intensity normalization was used. The excellent scan-rescan reproducibility (Fig. 2) demonstrates that our R1 estimates are robust, showing that quantitative in vivo maps can be obtained in only 25min of scan time (see Supplementary Information for further discussion and description of R1 mapping).

In recent work done in parallel with ours, cortical myelination has been estimated using a ratio between T1w and T2w scans (Glasser and Van Essen 2011). The results were broadly similar to ours but focused on large group results. Also, an additional series of post hoc normalization steps were required to visualize this non-quantitative data, which contain uncorrected artifacts in tissue contrast and brightness due to non-uniform B1+ transmit fields and to the effects of dividing pixel values from the 2 different image types that have contrasts sensitive to different aspects of microanatomy as well as different artifacts. The normalization steps are more difficult to reproduce at other centers than would be a hard threshold at cortical depth = 0.5 of R1= 0.62s−1.

Although local increases in cortical R1 corresponded well with known patterns of myelination and with retinotopically defined areas, there were some unexpected local increases in R1—for example, at the nearly adjoined center-of-gaze regions of V1, V2, V3, and hV4. Because R1 differences between cortical areas are comparable with R1 differences between different cortical depths (Fig. 1), detection of interareal differences requires uniformly accurate estimation of the gray–white and pial surfaces just as much as it requires R1 uniformity. To improve those estimates, the more uniform R1 images rather than T1w images could be used to reconstruct the pial surface. However, this will require extending the current FreeSurfer pial-surface-finding algorithm, which is finely adapted to T1w image intensity relations among tissue types. For example, in quantitative R1 images, there is no change in image intensity between the superior margin of the cerebellum and inferior margin of the inferior occipital lobe—which leads to a ballooning of the pial surface there—whereas in T1w images, the cerebellum and inferior occipital lobe are separated by a thin but distinct hypointensity corresponding to the dura—likely a mixed contrast artifact, but one that FreeSurfer relies on. Another region where the FreeSurfer pial-surface reconstruction using R1 images fails is near the thin temporal bone. It will be necessary to strategically incorporate quantitative proton density images into the pipeline to resolve this surface-reconstruction problem.

Another problem arises in areas such as S-I and V1 where parts of the cortical ribbon are so thin (∼1.5 mm) that they are spanned by less than two 0.8 mm3 voxels. The problem of misestimating the cortical depth fraction and hence misestimating R1 is particularly acute in this case, which might explain the unexpected reduction in R1 in parafoveal V1 (see medial surface views in Fig. 3).

Quantitative R1 mapping is fitted to answering a number of questions beyond those introduced here. For example, we are currently using it to characterize auditory and somatomotor fields (Dick et al. 2011) and are combining it with high-angular-resolution diffusion measurements in the cortex. This method could also illuminate structure/function mapping within cortical areas (e.g. in V2 stripes). Finally, the significant individual differences in local quantitative R1 we have uncovered may have functional correlates.

Supplementary Material

Supplementary material can be found at: http://www.cercor.oxfordjournals.org/.

Funding

This work was supported by the Wellcome Trust (the Wellcome Trust Centre for Neuroimaging is supported by core funding from the Wellcome Trust 091593/Z/10/Z), NIH R01 MH 081990 (M.I.S.), Royal Society Wolfson Research Merit Award (M.I.S.), Royal Society Research Grant (F.D.), and Capital Investment Funds from Birkbeck College and University College London (M.I.S. and F.D.). Funding to pay Open Access publication charges for this article was provided by the Wellcome Trust.

Notes

All authors participated in data collection, data analysis, and preparation of this paper. Conflict of Interest: None declared.

References

Allman
JM
Kaas
JH
A representation of the visual field in the caudal third of the middle temporal gyrus of the owl monkey (Aotus trivirgatus)
Brain Res
 , 
1971
, vol. 
31
 (pg. 
85
-
105
)
Amano
K
Wandell
BA
Dumoulin
SO
Visual field maps, population receptive field sizes, and visual field coverage in the human MT+ complex
J Neurophysiol
 , 
2009
, vol. 
102
 (pg. 
2704
-
2718
)
Annese
J
Gazzaniga
MS
Toga
AW
Localization of the human cortical visual area MT based on computer aided histological analysis
Cereb Cortex
 , 
2005
, vol. 
15
 (pg. 
1044
-
1053
)
Annese
J
Pitiot
A
Dinov
ID
Toga
AW
A myelo-architectonic method for the structural classification of cortical areas
Neuroimage
 , 
2004
, vol. 
21
 (pg. 
15
-
26
)
Bailey
P
von Bonin
G
The isocortex of man
 , 
1951
Urbana
University of Illinois Press
Bock
NA
Kocharyan
A
Liu
JV
Silva
AC
Visualizing the entire cortical myelination pattern in marmosets with magnetic resonance imaging
J Neurosci Meth
 , 
2009
, vol. 
185
 (pg. 
15
-
22
)
Braak
H
Architectonics of the human telencephalic cortex
 , 
1980
New York
Springer-Verlag
Caspers
S
Schleicher
A
Bacha-Trams
M
Palomero-Gallagher
N
Amunts
K
Zilles
K
Organization of the human inferior parietal lobule based on receptor architectonics
Cereb Cortex
 , 
2012
Clark
VP
Courchesne
E
Grafe
M
In vivo myeloarchitectonic analysis of human striate and extrastriate cortex using magnetic resonance imaging
Cereb Cortex
 , 
1992
, vol. 
2
 (pg. 
417
-
424
)
Dale
AM
Fischl
B
Sereno
MI
Cortical surface-based analysis. I. Segmentation and surface reconstruction
Neuroimage
 , 
1999
, vol. 
9
 (pg. 
179
-
194
)
DeYoe
EA
Carman
G
Bandettini
P
Glickman
S
Cox
R
Miller
D
Neitz
J
Mapping striate and extrastriate visual areas in human cerebral cortex
Proc Natl Acad Sci USA
 , 
1996
, vol. 
93
 (pg. 
2382
-
2386
)
Dick
F
Tierney
AT
Weiskopf
N
Lutti
A
Sereno
MI
Combined tonotopic and myeloarchitectonic mapping of human auditory regions
Neuroscience meeting planner
 , 
2011
Washington (DC)
 
Program No. 171.26
Draganski
B
Ashburner
J
Hutton
C
Kherif
F
Frackowiak
RSJ
Helms
G
Weiskopf
N
Regional specificity of MRI contrast parameter changes in normal ageing revealed by voxel-based quantification (VBQ)
Neuroimage
 , 
2011
, vol. 
55
 (pg. 
1423
-
1434
)
Eickhoff
S
Walters
NB
Schleicher
A
Kril
J
Egan
GF
Zilles
K
Watson
JDG
Amunts
K
High-resolution MRI reflects myeloarchitecture and cytoarchitecture of human cerebral cortex
Hum Brain Mapp
 , 
2005
, vol. 
24
 (pg. 
206
-
215
)
Engel
SA
Rumelhart
DE
Wandell
BA
Lee
AT
Glover
GH
Chichilnisky
EJ
Shadlen
MN
fMRI of human visual cortex
Nature
 , 
1994
, vol. 
36
 pg. 
525
 
Fischl
B
Dale
AM
Measuring the thickness of the human cerebral cortex from magnetic resonance images
Proc Natl Acad Sci USA
 , 
2000
, vol. 
97
 (pg. 
11050
-
11055
)
Fischl
B
Rajendran
N
Busa
E
Augustinack
J
Hinds
O
Yeo
BTT
Mohlberg
H
Amunts
K
Zilles
K
Cortical folding patterns and predicting cytoarchitecture
Cereb Cortex
 , 
2008
, vol. 
18
 (pg. 
1973
-
1980
)
Fischl
B
Sereno
MI
Tootell
RB
Dale
AM
High-resolution intersubject averaging and a coordinate system for the cortical surface
Hum Brain Mapp
 , 
1999
, vol. 
8
 (pg. 
272
-
284
)
Flechsig
P
Antomie des menschlichen Gehirns und Rückenmarks auf myelogenetischer Grundlage
 , 
1920
Leipzig
Georg Thieme
Gamberini
M
Galletti
C
Bosco
A
Breveglieri
R
Fattori
P
Is the medial posterior parietal area V6A a single functional area?
J Neurosci
 , 
2011
, vol. 
31
 (pg. 
5145
-
5157
)
Geyer
S
Weiss
M
Reimann
K
Lohmann
G
Turner
R
Microstructural parcellation of the human cerebral cortex – from Brodmann's post-mortem map to in vivo mapping with high-field magnetic resonance imaging
Front Hum Neurosci
 , 
2011
, vol. 
5
 
Glasser
MF
Van Essen
DC
Mapping human cortical areas in vivo based on myelin content as revealed by t1- and t2-weighted MRI
J Neurosci
 , 
2011
, vol. 
31
 (pg. 
11597
-
11616
)
Graziano
MS
Gandhi
S
Location of the polysensory zone in the precentral gyrus of anesthetized monkeys
Exp Brain Res
 , 
2000
, vol. 
135
 (pg. 
259
-
266
)
Hadjikhani
N
Liu
AK
Dale
AM
Cavanagh
P
Tootell
RB
Retinotopy and color sensitivity in human visual cortical area V8
Nat Neurosci
 , 
1998
, vol. 
1
 (pg. 
235
-
241
)
Hagler
DJ
Riecke
L
Sereno
MI
Parietal and superior frontal visuospatial maps activated by pointing and saccades
Neuroimage
 , 
2007
, vol. 
35
 (pg. 
1562
-
1577
)
Hagler
DJ
Sereno
MI
Spatial maps in frontal and prefrontal cortex
Neuroimage
 , 
2006
, vol. 
29
 (pg. 
567
-
577
)
Helms
G
Dathe
H
Dechent
P
Quantitative FLASH MRI at 3T using a rational approximation of the Ernst equation
Magn Reson Med
 , 
2008
, vol. 
59
 (pg. 
667
-
672
)
Helms
G
Dathe
H
Weiskopf
N
Dechent
P
Identification of signal bias in the variable flip angle method by linear display of the algebraic Ernst equation
Magn Reson Med
 , 
2011
, vol. 
66
 (pg. 
669
-
677
)
Helms
G
Draganski
B
Frackowiak
R
Ashburner
J
Weiskopf
N
Improved segmentation of deep brain grey matter structures using magnetization transfer (MT) parameter maps
Neuroimage
 , 
2009
, vol. 
47
 (pg. 
194
-
198
)
Helms
G
Finsterbusch
J
Weiskopf
N
Dechent
P
Rapid radiofrequency field mapping in vivo using single-shot STEAM MRI
Magn Reson Med
 , 
2008
, vol. 
60
 (pg. 
739
-
743
)
Huang
R-S
Sereno
MI
Dodecapus: an MR-compatible system for somatosensory stimulation
Neuroimage
 , 
2007
, vol. 
34
 (pg. 
1060
-
1073
)
Huk
AC
Dougherty
RF
Heeger
DJ
Retinotopy and functional subdivision of human areas MT and MST
J Neurosci
 , 
2002
, vol. 
22
 (pg. 
7195
-
7205
)
Kolster
H
Peeters
R
Orban
GA
The retinotopic organization of the human middle temporal area MT/V5 and its cortical neighbors
J Neurosci
 , 
2010
, vol. 
30
 (pg. 
9801
-
9820
)
Lashley
KS
Clark
G
The cytoarchitecture of the cerebral cortex of Ateles; a critical examination of architectonic studies
J Comp Neurol
 , 
1946
, vol. 
85
 (pg. 
223
-
305
)
Lutti
A
Hutton
C
Finsterbusch
J
Helms
G
Weiskopf
N
Optimization and validation of methods for mapping of the radiofrequency transmit field at 3T
Magn Reson Med
 , 
2010
, vol. 
64
 (pg. 
229
-
238
)
Lutti
A
Stadler
J
Josephs
O
Windischberger
C
Speck
O
Bernarding
J
Hutton
C
Weiskopf
N
Robust and fast whole brain mapping of the RF transmit field B1 at 7T
PLoS One
 , 
2012
, vol. 
7
 pg. 
e32379
 
Malikovic
A
Amunts
K
Schleicher
A
Mohlberg
H
Eickhoff
SB
Wilms
M
Palomero-Gallagher
N
Armstrong
E
Zilles
K
Cytoarchitectonic analysis of the human extrastriate cortex in the region of V5/MT+: a probabilistic, stereotaxic map of area hOc5
Cereb Cortex
 , 
2007
, vol. 
17
 (pg. 
562
-
574
)
Merzenich
MM
Kaas
JH
Principles of organization of sensory-perceptual systems in mammals
Prog Psychobiol Physiol Psych
 , 
1980
, vol. 
9
 (pg. 
1
-
42
)
Meynert
T.
Der Bau der Grosshirnrinde und seine örtlichen Verschiedenheiten nebst einem pathologisch-anatomischen Corollarium
Vjschr Psychiat
 , 
1867
, vol. 
1
 (pg. 
198
-
217
)
Pitzalis
S
Galletti
C
Huang
R-S
Patria
F
Committeri
G
Galati
G
Fattori
P
Sereno
MI
Wide-field retinotopy defines human cortical visual area V6
J Neurosci
 , 
2006
, vol. 
26
 (pg. 
7962
-
7973
)
Preibisch
C
Deichmann
R
Influence of RF spoiling on the stability and accuracy of T1 mapping based on spoiled FLASH with varying flip angles
Magn Reson Med
 , 
2009
, vol. 
61
 (pg. 
125
-
135
)
Schmierer
K
Scaravilli
F
Altmann
DR
Barker
GJ
Miller
DH
Magnetization transfer ratio and myelin in postmortem multiple sclerosis brain
Ann Neurol
 , 
2004
, vol. 
56
 (pg. 
407
-
415
)
Sereno
MI
Dale
AM
Reppas
JB
Kwong
KK
Belliveau
JW
Brady
TJ
Rosen
BR
Tootell
RB
Borders of multiple visual areas in humans revealed by functional magnetic resonance imaging
Science
 , 
1995
, vol. 
268
 (pg. 
889
-
893
)
Sereno
MI
Huang
R-S
A human parietal face area contains aligned head-centered visual and tactile maps
Nat Neurosci
 , 
2006
, vol. 
9
 (pg. 
1337
-
1343
)
Sereno
MI
Pitzalis
S
Martinez
A
Mapping of contralateral space in retinotopic coordinates by a parietal cortical area in humans
Science
 , 
2001
, vol. 
294
 (pg. 
1350
-
1354
)
Sigalovsky
IS
Fischl
B
Melcher
JR
Mapping an intrinsic MR property of gray matter in auditory cortex of living humans: a possible marker for primary cortex and hemispheric differences
Neuroimage
 , 
2006
, vol. 
32
 (pg. 
1524
-
1537
)
Silver
MA
Kastner
S
Topographic maps in human frontal and parietal cortex
Trends Cogn Sci
 , 
2009
, vol. 
13
 (pg. 
488
-
495
)
Smart
IH
McSherry
GM
Gyrus formation in the cerebral cortex of the ferret. II. Description of the internal histological changes
J Anat
 , 
1986
, vol. 
147
 (pg. 
27
-
43
)
Swisher
JD
Halko
MA
Merabet
LB
McMains
SA
Somers
DC
Visual topography of human intraparietal sulcus
J Neurosci
 , 
2007
, vol. 
27
 (pg. 
5326
-
5337
)
Walters
NB
Egan
GF
Kril
JJ
Kean
M
Waley
P
Jenkinson
M
Watson
JDG
In vivo identification of human cortical areas using high-resolution MRI: an approach to cerebral structure-function correlation
Proc Natl Acad Sci USA
 , 
2003
, vol. 
100
 (pg. 
2981
-
2986
)
Wandell
BA
Dumoulin
SO
Brewer
AA
Visual field maps in human cortex
Neuron
 , 
2007
, vol. 
56
 (pg. 
366
-
383
)
Weiskopf
N
Lutti
A
Helms
G
Novak
M
Ashburner
J
Hutton
C
Unified segmentation based correction of R1 brain maps for RF transmit field inhomogeneities (UNICORT)
Neuroimage
 , 
2011
, vol. 
54
 (pg. 
2116
-
2124
)
Xu
G
Knutsen
AK
Dikranian
K
Kroenke
CD
Bayly
PV
Taber
LA
Axons pull on the brain, but tension does not drive cortical folding
J Biomech Eng.
 , 
2010
, vol. 
132
 (pg. 
071013–1
-
071013–8
)
Zilles
K
Amunts
K
Centenary of Brodmann's map–conception and fate
Nat Rev Neurosci
 , 
2010
, vol. 
11
 (pg. 
139
-
145
)

Author notes

Some of these results were presented as poster 580.18/SS17 at the November 2010 Society for Neuroscience meetings.
N.W. and F.D. shared senior authors.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/), which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.