Freesurfer Retinotopic Analysis Essay

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.

MT, myelination, parcellation, surface reconstruction, visual areas

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.

top



  • The purpose of this FAQ is to provide a Wiki space where users can add the most frequent questions to avoid asking already answered questions in the support list.

Contents

  1. General
    1. Q. How can I help this FAQ?
    2. Q. What are the advantages of FreeSurfer over VBM?
    3. Q. How do I know the version of FreeSurfer I'm running?
    4. Q. Is the FreeSurfer source code available?
  2. Computing
    1. Q. How long does it take to finish a reconstruction?
    2. Q. How can I reduce the time of recon-all in a group of patients?
    3. Q. I want to use freesurfer to quantify cortical thickness, and I want to use the software under Linux. What equipment for computer do you recommend?
    4. Q. Is it possible to run FreeSurfer in Ubuntu Linux?
    5. Q. Could I run multiple instances of Freesurfer (on my virtual box)?
  3. Acquisition
    1. Q. Is it recommended that people use memprage? How are they analyzed? Just sqrt sum sqr of the echoes? Or is there something more elaborate?
    2. Q. Are there suggested scan sequences which you recommend I use with Freesurfer?
    3. Q. Is it possible to analyze clincial structural MRI exams with help of Freesurfer?
  4. Processing/Re-processing Data in FreeSurfer
    1. Q. Can I run some cases in my dataset using one version of FreeSurfer and others using a different version of FreeSurfer?
    2. Q. How to run recon-all in a ssh terminal
    3. Q. I had a [Power Outage|Computer Failure|Spilled coffee on computer|etc], how can I resume the recon-all?
    4. Q. I have a subject which has been running for a really long time in recon-all, how can I tell if there is something wrong?
    5. Q. I have already skull-stripped data. Can I submit it to recon-all?
    6. Q. I made white matter and pial edits to my volume. Do I need to run -autorecon2-wm and when that is finished run -autorecon2-pial?
    7. Q. One of my cases doesn't have the cerebellum so I adjusted the watershed threshold to fix it. However, at one threshold, the cerebellum is still not included, and at the next the entire skull is there. What should I do?
    8. Q. One or more of the files in a subject's directory was accidentally deleted or has become corrupted. How do I recreate the missing file(s)?
  5. Common Error messages
    1. Q. Help! I got this error message: "mri_watershed error: GLOBAL region of the brain empty!" - what should I do?
    2. Q. I get an error message from the Talairach Failure Detection. What does this mean?
    3. Q. I get an error during the mri_ca_label step while running recon-all. The last thing it says is: "saving intensity scales to aseg.auto_noCCseg.label_intensities.txt". What's wrong?
    4. Q. I get an error during the Talairach transform step when running recon-all on ANALYZE images. Why is this happening?
  6. FreeSurfer Output Questions
    1. Q. The surfaces near the medial wall, hippocampus, and amygdala aren't accurately following the structures there. How can I fix this?
    2. Q: It seems that there is a 5mm thickness upper limit in my volumes. Is it normal? How can I change this?
    3. Q. Why in many subjects the insular cortical surface seems so thick? Is the convoluted nature of the Insula that causes that?
    4. Q. Why would the orientation of my scans look wrong after processing my Siemens DICOM files with 'mri_convert'?
    5. Q. What are the sulc and curv overlays (in QDEC) showing us?
    6. Q. How would I put a FreeSurfer output volume and/or segmentation back into the same space as my original anatomical input (native space)?
    7. Q. How do I transform coordinates in one space to those in another space (eg, a point on the surface to MNI305 space or to the col, row, slice of a functional volume)?
    8. Q. How can I measure the distance along the cortical surface between two points located on the cortical surface?
    9. Q. The pial surface includes some cerebellum. How can I fix this?
    10. Q. How do I get the cortical thickness maps for several subjects into the same template space (fsaverage)?
    11. Q. How can I get a high resolution atlas of the cortex in Freesurfer?
    12. Q. Where can I find V1 labels?
  7. Analysis of FreeSurfer Data
    1. Q. I am trying to measure the cortical thickness of a specific ROI. How can I do this?
    2. Q. I am using QDEC to examine the anatomical differences between two groups of subjects. The surface-based measures I can select are thickness, area, area.pial, sulc, curv, and jacobian_white. Could anybody tell me what anatomical features the later three (sulc, curv, and jacobian_white) actually measure?
    3. Q. I am trying to measure the global mean cortical thickness (i.e combined across hemispheres). How would I do this?
    4. Q: How can I get measurements of the lobes (parietal, temporal, frontal, & occipital)?
    5. Q. How would I calculate the total CSF volume?
    6. Q. What is the unit measure of mean curvature and Gaussian curvature?
    7. Q. I want to run make_average_subject in order to prepare my data for GLM analysis, which file should I specify for the -xform flag, and what are the differences between my options?
    8. Q. What goes into the calculation of subcortical volume? We have tried adding volumes of individual subcortical areas but the total does not equal the subcortical volume provided by aseg. Is the subcortical volume calculation accurate?
    9. Q. How would you perform a power analysis for a whole brain group comparison QDEC analysis?
    10. Q. How can I get the volume of the different lobes of the brain (occipital, parietal, temporal, etc) ? White matter and gray matter? Labels ?
    11. Q. Why do I get different results when scanning the same person twice?
  8. FreeSurfer & Matlab
    1. Q. How can I use the Freesurfer Matlab commands if I have a copy of Matlab installed on a Windows machine?
    2. Q. Can I load FreeSurfer output in Matlab?
    3. Q. How can I make a histogram of cortical thickness?
    4. Q. How can I obtain the thickness of each vertex, and how can you identify which structure each vertex belongs to?
  9. FreeSurfer GUI
    1. Q. How can I troubleshoot rendering problems when using TKsurfer?

General

Q. How can I help this FAQ?

A: If you are able to edit pages, go ahead! If you don't have permission to write on the Wiki, send an e-mail to ppj at netfilter dot com dot br

Q. What are the advantages of FreeSurfer over VBM?

A:

  1. FS uses geometry to do inter-subject registration, which experience has shown results in a much better matching of homologous cortical regions than volumetric techniques.
  2. FS allows you to look at the two components of volume separately (thickness and surface area). It has been found that these two do not necessarily track one another, and in the worst case where one is increasing and the other decreasing the volume change can be 0.
  3. The target that FS uses for registration (the white matter surface geometry) is completely invariant to gm atrophy, so gm changes won't result in a different registration.

Q. How do I know the version of FreeSurfer I'm running?

A: Run:

Q. Is the FreeSurfer source code available?

A: Yes, the source code is accessible via CVS access using the instructions described here.

It is also recommended that you subscribe to the freesurfer mailing list to post questions and monitor solutions from other users.


Computing

Q. How long does it take to finish a reconstruction?

A: It depends on your processor speed and machine performance (notice that dual/quad/hex core or hyperthead won't speed up one analysis significantly), give a look at the link below in the section "Step-wise directives": http://surfer.nmr.mgh.harvard.edu/fswiki/recon-all

Also you can contribute to our running time statistics: https://surfer.nmr.mgh.harvard.edu/fswiki/ReconAllRunTimes

Q. How can I reduce the time of recon-all in a group of patients?

A: FreeSurfer run its process in a non-parallel environment, so you won't have benefit from a dual/quad/hex core machine for a single case analysis. However if you have many cases you can start two FreeSurfer recon-all process in the same machine and theoretically you can reduce by half the time to analyze your group of cases. A similar procedure can also be used in quad/hex-core environment. Note that to take benefit of a multi-core environment you need to use a SMP kernel in your OS.

Q. I want to use freesurfer to quantify cortical thickness, and I want to use the software under Linux. What equipment for computer do you recommend?

This answer requires constant updates:

For Best Performances:

Use a multi-Core processor Intel (a i7-900 or newer series processor or Xeon 6500/7500 Series)

Install at Least 8GB of Memory

If you prefer AMD you can use:

AMD Opteron 12-Core or Phenon II X6

Other alternative much more expensive is: Mac Pro with 12 cores and 16GB RAM.

The number of cores is roughly the number of studies you can process simultaneously. Notice that each process will take 20-24 hrs.

Keep in mind that you need to use the fastest memory in order to achieve maximum benefit from multi-core architecture.

Q. Is it possible to run FreeSurfer in Ubuntu Linux?

A: Yes. Ubuntu Linux is basically a Debian distro, so you should use FreeSurfer RH9 version. Depending on your video card you should disable the DRI using the option NoDRI in the Device section of your X configuration file. Notice that in older Ubuntu version there's a bug that prevents NoDRI from working.

Q. Could I run multiple instances of Freesurfer (on my virtual box)?

A: Yes, it is certainly possible (even advisable) to run more than one copy of Freesurfer on a machine. It depends on how many processor cores you have. You'll need at least 2Gb per process if not 3Gb, one core per process. It may be possible to run one process using 1Gb if you use the 'no-gcut' flag. Open one terminal window for each instance of FreeSurfer, using serperate commands in each terminal window. When running four individual recon-all jobs on four cores, we have found that they will all complete in about 110% of the time of a single run.

When using a virtual box, assign N-1 cores to your virtual box, where N is the total number of processor cores you have, along with at least (N-1) x 2Gb of RAM.


Acquisition

Q. Is it recommended that people use memprage? How are they analyzed? Just sqrt sum sqr of the echoes? Or is there something more elaborate?

A: Yes, particularly for longitudinal. The increased bandwidth makes a huge difference in being able to register across time. And at the moment we do use the rms. An optimal combo would maybe be a tiny bit better, but probably not much difference.

Q. Are there suggested scan sequences which you recommend I use with Freesurfer?

A: Yes, Custom multiecho sequences for Siemens scanners are available from the Martinos Center. The multiecho FLASH is available as a standard sequence on the Siemens platform and we distribute a slightly modified version that makes it more convenient to set up the echo timing and allows a little more flexibility w.r.t. geometry and spoiling but this isn't critical. We also distribute a multiecho version of the MPRAGE sequence that isn't a product sequence yet.

We are happy to send you binaries but MGH and Siemens requires that you sign a "C2P" agreement, which is an indemnification document, before we can provide you with the sequences and protocols for free.

Q. Is it possible to analyze clincial structural MRI exams with help of Freesurfer?

A: Yes, if the resolution is around 1mm, it should be possible. If the resolution is greater than 1.3mm, then surfaces will probably not be valid, but subcortical structures might be valid.


Processing/Re-processing Data in FreeSurfer

Q. Can I run some cases in my dataset using one version of FreeSurfer and others using a different version of FreeSurfer?

A: No, mixing versions is never a good idea as results are expected to differ, if only slightly. The same version of FreeSurfer should be used to process all cases within a dataset. Another consideration is that someone editing data run with an older version of FreeSurfer will probably have more edits to do than if the same case was run with the most recent version.

If you are sourcing the stable or dev version of FreeSurfer from within the NMR center, you are strongly encouraged to created a frozen copy of FreeSurfer for your study. To do that, just download the public distribution from our site and install it in your SUBJECTS_DIR or some other project-appropriate place.

Q. How to run recon-all in a ssh terminal

A: This is tricky, because recon-all spawns many process and redirect I/O to the terminal that opened it. Run

  • recon-all -s subjid -all > /dev/null 2>&1 </dev/null & disown -a

Q. I had a [Power Outage|Computer Failure|Spilled coffee on computer|etc], how can I resume the recon-all?

A: You can try

recon-all -make all -s subject

Also, depending on your version, there might be a file in the scripts directory called "IsRunning". You should delete this if it is there.

Q. I have a subject which has been running for a really long time in recon-all, how can I tell if there is something wrong?

A: Assuming you have processed other subjects in considerably less time, this probably means that there is a giant defect that will not be corrected properly (eg. cerebellum or skull attached). You should check the ?h.inflated.nofix or ?h.orig.nofix surfaces for any large defects. If you load the filled.mgz you should be able to see whether the cerebellum is still attached, and if so you would want to edit the wm.mgz. Also, verify that the pons and corpus collosum were properly detected.

Q. I have already skull-stripped data. Can I submit it to recon-all?

A: If your skull-stripped volume does not have the cerebellum, then no. If it does, then yes, however you will have to run the data a bit differently.

First you must run only -autorecon1 like this:
recon-all -autorecon1 -noskullstrip -s <subjid>

Then you will have to make a symbolic link or copy T1.mgz to brainmask.auto.mgz and a link from brainmask.auto.mgz to brainmask.mgz. Finally, open this brainmask.mgz file and check that it looks okay (there is no skull, cerebellum is intact; use the sample subject bert that comes with your FreeSurfer installation to make sure it looks comparable). From there you can run the final stages of recon-all:
recon-all -autrecon2 -autorecon3 -s <subjid>

Q. I made white matter and pial edits to my volume. Do I need to run -autorecon2-wm and when that is finished run -autorecon2-pial?

A: No. If you made both white matter and pial edits, you only need to run -autorecon2-wm which is higher up in the processing stream (before the steps of -autorecon2-pial). By running -autorecon2-wm, it will perform the steps necessary to fix a white matter edit and then all the remaining steps in -autorecon2 including those steps that -autorecon2-pial would run. The OtherUsefulFlags wiki shows the hierarchy of the -autorecon2 flag shortcuts (-autorecon2-cp starts before -autorecon2-wm which starts before -autorecon2-pial). In general, you should start the processing stream at the earliest step where you made an manual intervention and the rest will be taken care of. Don't forget to run -autorecon3!

Q. One of my cases doesn't have the cerebellum so I adjusted the watershed threshold to fix it. However, at one threshold, the cerebellum is still not included, and at the next the entire skull is there. What should I do?

A: Removal of cerebellum is always bad and must be fixed (cerebellum is included in the subcortical segmentation portion and is needed for a proper atlas alignment). Including skull is not good, but leaving skull is acceptable. The problem being sometimes the wm seg catches skull. So make adjustments till you get cerebellum, then just see what happens with the surfaces given remaining skull, and edit out skull in brainmask.mgz if it hurts it.

For the command line below, try it with and without the -no-wgcaatlas flag. Generally the -wsthresh only needs to be varied.

  • recon-all -skullstrip -wsthresh 35 -clean-bm -no-wsgcaatlas -subjid subjid

If that still doesn't work, try using the multistrip command with the watershed method. This will give an output of 4 sets of different thresholds each for the orig, nu, and T1 volume. You can choose the best one from the set and copy it to brainmask.auto.mgz and run the last step of -skullstrip. More info can be found here: http://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/SkullStripFix

Q. One or more of the files in a subject's directory was accidentally deleted or has become corrupted. How do I recreate the missing file(s)?

A: You will need to rerun the recon-all step in which that file is first created. For example, to recreate the norm.mgz volume you would need to run recon-all -canorm -s <subjid>. Click here to find out in which step a certain file is created in the version 5.3 recon-all process flow.


Common Error messages

Q. Help! I got this error message: "mri_watershed error: GLOBAL region of the brain empty!" - what should I do?

A: Run this:
recon-all -skullstrip -no-wsgcaatlas -s <subjid>

If that goes through without error, you can finish processing your data by running:
recon-all -autorecon2 -autorecon3 -s <subjid>

Q. I get an error message from the Talairach Failure Detection. What does this mean?

A: In certain versions of FreeSurfer, the talairach failure detection is too conservative so this may not be a problem. First, you will want to check the talairach transform to make sure it looks okay. Directions on how to do that are here: https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/Talairach

Essentially, you want to use the command:

tkregister2 --mgz --s <subjid> --fstal

and see if the green lines line up with the anatomy of the blurry brain image. If so, then you can run this subject as you did before, just add -notal-check to your command. For example, recon-all -all -s subjid -notal-check. If the green lines do not line up well with the anatomy, you will want to follow directions on the wiki page mentioned above to fix it. If the images you passed to recon-all are in Analyze format, verify that the orientation is correct.

Q. I get an error during the mri_ca_label step while running recon-all. The last thing it says is: "saving intensity scales to aseg.auto_noCCseg.label_intensities.txt". What's wrong?

A: This is a bug in v.4.5. A fixed mri_ca_label can be downloaded from here:
ftp://surfer.nmr.mgh.harvard.edu/pub/dist/freesurfer/misc/linux-centos4_x86_64/

copy it your $FREESURFER_HOME/bin

We would advise rerunning all subjects with this fixed version.

Q. I get an error during the Talairach transform step when running recon-all on ANALYZE images. Why is this happening?

A: Images in the ANALYZE format do not retain orientation information (left/right). You can bypass the Talairach check by including the '-notal-check' flag in your recon-all command. Alternatively, you can specify the orientation (neurological or radiological) using mri_convert. Run 'mri_convert --help' to get info on specifying volume orientation in the description section. Lastly, you can manually register your brain volume in Talairach space by following this tutorial.


FreeSurfer Output Questions

Q. The surfaces near the medial wall, hippocampus, and amygdala aren't accurately following the structures there. How can I fix this?

A: The good news is you don't have to! These areas are generally unreliable for a thickness study so you would want to exclude them from your analysis. If you use glmfit, it calls the ?h.cortex.label file in order to determine thickness. This file automatically has set these areas to zero so they will not be included in your analysis. To see what area ?h.cortex.label excludes you can load it in tksurfer on top of the inflated surface. If you have tkmedit open at the same time for the same subject, you can use the Go To Saved Point function to make sure areas of concern are in this excluded section. *Note: If you are using a FreeSurfer version older than 4.3.0, the excluded region includes the insula. If you are using an analysis program other than glmfit, just be sure to use the ?h.cortex.label file to get your thickness measurements.

Several reasons make it difficult to generate pial/wm surfaces in the medial temporal lobe area. This region tends to be furthest away from the coils receiving the MR signal, potentially late to myelinate in individuals under 20 years old, very thin white matter in general, and is in close proximity to regions of susceptibility.

The picture below shows the ?h.cortex.label (outlined in yellow) over the ?h.aparc.annot. Notice how the ?h.cortex.label includes some of the ?h.unknown.label:

As always, use the example subject bert that comes with the FreeSurfer installation to see what we define as acceptable surfaces.

Q: It seems that there is a 5mm thickness upper limit in my volumes. Is it normal? How can I change this?

A: We did this to prevent noncortical regions such as the basal ganglia from corrupting the thickness measure through averaging. With the ?h.cortex.label it is probably no longer needed, however you can use mris_thickness -max <max thick> to generate a thickness with a different max. Take a look in the example below:

mris_thickness -max 10 bert lh newlh.thickness

The file newlh.thickness will be created inside the surf directory of your subject. But we don't think the true thickness is ever that much except in pathological cases like dysplasia and other disorders of cortical development.

Q. Why in many subjects the insular cortical surface seems so thick? Is the convoluted nature of the Insula that causes that?

A: It's not the convoluted nature of the Insula. It's the fact that extreme capsule is so thin that it frequently isn't very apparent on MR, and so there appears to be continuous gray matter from the basal ganglia into the cortex. We think version 4.0 fixes this.

Q. Why would the orientation of my scans look wrong after processing my Siemens DICOM files with 'mri_convert'?

A: It can be misleading to check orientation using FSLVIEW because it orients volumes based on the way they are on disk so you often get things looking pretty strange (e.g., upside-down). FSLVIEW does dispaly little letters to indicate what it thinks the orientation is, so if you stick to those, then you can properly judge whether the orientation is correct.

Q. What are the sulc and curv overlays (in QDEC) showing us?

A: The 'sulc' conveys information on how far removed a particular vertex point on a surface is from a hypothetical "mid-surface" that exists between the gyri and sulci. This surface is chosen so that the "mean" of all these displacements is zero. The 'sulc' gives a indication then of linear distance and displacements: how "deep" and how "high" are brain folds

In FreeSurfer, gyri have negative 'sulc' values, are colored green, and indicate how far "down" a point has to travel to reach this "mid-surface". Sulci have positive 'sulc' values, are colored red, and indicate how far "up" a point needs to travel to reach the mid-surface.

The 'curv' conveys information on the curvature (not distance) at a specific vertex point. The color conveys the sign, and is just an arbitrary choice. The sharper the curve, the higher the value (positive or negative). Areas with positive curvature, are colored red, and correspond to curvatures in sulci, i.e. curving "up". Areas with negative curvature are colored green, and correspond to curves pointing "down", i.e. gyri.

So, in a nutshell, the difference is that the 'curv' files contain information about curvatures, and the 'sulc' files contain information about displacement.

Q. How would I put a FreeSurfer output volume and/or segmentation back into the same space as my original anatomical input (native space)?

A: Please see directions on how to do this here, and also on this page http://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems

Q. How do I transform coordinates in one space to those in another space (eg, a point on the surface to MNI305 space or to the col, row, slice of a functional volume)?

A: Directions on how to do this are on this page: http://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems

Q. How can I measure the distance along the cortical surface between two points located on the cortical surface?

A: There is a tool which allows you to do this called 'mris_pmake'. Please read the associated help file for details on how to use it ('mris_pmake --help').

Q. The pial surface includes some cerebellum. How can I fix this?

A: This fix only works for v4.3 and up. If you are using an earlier version, please contact us for help.

In order to fix this problem, you will have to edit the brain.finalsurfs.mgz (and not the brainmask.mgz). Remove the parts of the cerebellum that are affecting the surfaces. Save your changes and then run:
recon-all -make all -s subjid

Q. How do I get the cortical thickness maps for several subjects into the same template space (fsaverage)?

A: Include the '-qcache -measure thickness' flags in your recon-all command for each subject. This will create files in the /surf directory which sample the thickness data at different smoothing levels onto the fsaverage subject space. If the '-measure' flag is not included, all cortical measurements will be sampled onto fsaverage and smoothed (thickness, sulc, area, curv, etc.).

Q. How can I get a high resolution atlas of the cortex in Freesurfer?

A: This can be accomplished using 'mris_divide_parcellation'. You can pass your preferred max area of a parcellation and it will keep subdividing along the primary eigen-axis until no units are above that surface area.

Q. Where can I find V1 labels?

A: V1 labeling is produced by default in recon-all as part of the Brodmann area set (located in ../fsaverage/label/?h.V1.label). Those labels are described here. The Hind's V1 labeling method must be run seperately by adding the '-label-v1' flag to your recon-all command:

recon-all -all -s bert -label-v1

The accuracy of the Hinds and Fischl V1 labeling method compared to retinotopy depends under what conditions the retinotopy is carried out (field strength, # of coils elements, voxel size, etc.). In our comparision, the folding patterns predict the border location to about 2.5mm, which is probably better than you can get with retinotopy. (click here for tips on creating retinotopy stimuli)


Analysis of FreeSurfer Data

Q. I am trying to measure the cortical thickness of a specific ROI. How can I do this?

A: You can save the ROI to a label file and then use:

  • mris_anatomical_stats -l <label file>

Q. I am using QDEC to examine the anatomical differences between two groups of subjects. The surface-based measures I can select are thickness, area, area.pial, sulc, curv, and jacobian_white. Could anybody tell me what anatomical features the later three (sulc, curv, and jacobian_white) actually measure?


sulc = "average convexity" from our 1999 reconII paper(https://surfer.nmr.mgh.harvard.edu/ftp/articles/fischl99b-recon2.pdf).
Essentially measures the depth/height of each point above the average surface.
curv = smoothed mean curvature.
jacobian_white = the jacobian of the spherical transform. Measures the amount of distortion needed to warp a subject into register with the atlas.
You might want to look at the slides downloadable from the top of this page:

http://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial
which contain pictures showing the meaning of 'curv' and 'sulc'.

Q. I am trying to measure the global mean cortical thickness (i.e combined across hemispheres). How would I do this?

A: One suggestion is to use the surface area of each hemisphere as the weighting factor. In which case the global mean thickness including both hemispheres would be given by:

bh.thickness = ( (lh.thickness * lh.surfarea) + (rh.thickness * rh.surfarea) ) / (lh.surfarea + rh.surfarea)

If you use the values in the ?h.aparc.stats, it already factors out the 'unknown' region, so you don't have to do it yourself.

Q: How can I get measurements of the lobes (parietal, temporal, frontal, & occipital)?

A: You can either use the PALS_B12 atlas (mapped to fsaverage in v5.0) or create your own following the directions on the CorticalParcellation wiki.

Q. How would I calculate the total CSF volume?

A: You can add up the various ventricular structures to get total ventricular volume, but we don't segment sulcal CSF, since it's not distinguishible from bone on a T1-weighted MRI.

Q. What is the unit measure of mean curvature and Gaussian curvature?

A: See these wiki pages for more info: Mean curvature, Gaussian curvature

Q. I want to run make_average_subject in order to prepare my data for GLM analysis, which file should I specify for the -xform flag, and what are the differences between my options?

A: Your typical options are talairach.lta, talairach.xfm, and talairach.m3z. The .lta and .xfm are linear transforms to a Talairach coordinate system. The .m3z is a non-linear morph and will give you a much higer anatomical resolution.

Q. What goes into the calculation of subcortical volume? We have tried adding volumes of individual subcortical areas but the total does not equal the subcortical volume provided by aseg. Is the subcortical volume calculation accurate?

A: The aseg.stats file takes into account partial voluming which is not taken into account when you simply sum up the subcortical structures because some of those structures will be partially volumed with white matter.

Q. How would you perform a power analysis for a whole brain group comparison QDEC analysis?

A: For the power analysis, you need four things:

  1. Effect Size
  2. Number of Subjects
  3. Target False Positive Rate (alpha)
  4. Target False Negative Rate (beta)

Given any 3, you can compute the 4th. You can get the Effect Size from the output of the QDEC analysis. Each analysis creates a GLM directory, and there is a directory for each contrast in the GLM dir. In the contrast dir, you will find several files, but the important ones for this are the gamma.mgh and the gammavar.mgh. The gammavar is the square of the error bar (ie, t=gamma/sqrt(gammavar)). The gamma will not change as you add subjects (at least in expectation). The gammavar will drop linearly with the number of subjects (again in expectation). The effect size will then be gamma/sqrt(gammavar*Npilot), where Npilot is the number subjects in the pilot study. In a new study with Nnew subjects, the expected t will be tnew = gamma/sqrt(gammavar*Npilot/Nnew).

Q. How can I get the volume of the different lobes of the brain (occipital, parietal, temporal, etc) ? White matter and gray matter? Labels ?

A: For V5.0 and later, you can run mri_annotation2label with --lobesStrict to get a lobe annotation. If that definition of "lobes" is good for you, then you can run mris_anatomical_stats to get the volume for each lobe. For versions prior to 5.0, you can use mri_annotation2label to break the labels apart, then use mri_mergelabels to combine the individual labels into lobe labels, then use mris_label2annot to create a lobe annotation, then use mris_anatomical_stats.

Q. Why do I get different results when scanning the same person twice?

A: There are a variety of sources of variance. First there are some reasons that affect the images. Different looking images will of course affect the structural estimates. The second set of reasons is related to random noise and processing bias.

a) Acquisition

a) 1. Head Motion during acquisition changes the image and results in smaller GM estimates (even for images that pass QC)
http://reuter.mit.edu/publications/pid/reuter-motion14
http://reuter.mit.edu/publications/pid/tisdall15

a) 2. Hydration levels affect the image and e.g. GM estimates
http://reuter.mit.edu/publications/pid/biller15

a) 3. Different head positions can cause large differences if images are not gradient unwarped (depending on your hardware these effects can be small or large). Also there is of course noise in the images which can lead to different results.

a) 4. Different person: sometimes errors in the de-identification (annonymization) pipeline can causes images with the same subject id to come from different subjects. Also people sometimes send a sibling or friend to a follow-up visit (in order to complete the study and get paid, if they cannot come themselves).

b) Processing

b) 1. Different noise or small changes in the image can trigger processing to do different things (e.g. small skull strip failure, different surface placement, e.g. one time including dura the other time not). Small changes can accumulate. This can produce arbitrary large changes in a single subject. It can be reduced/fixed via manual edits.

b) 2. Processing bias, for example, mapping the follow-up to baseline (which some people do as a pre-processing step), is problematic as follow-ups get interpolated (which looks like smoothing) and so you introduce a change to some of your images (the follow-up only), biasing your measurements of change. See
http://reuter.mit.edu/publications/pid/reuter-long12
http://reuter.mit.edu/publications/pid/reuter-bias11
This problem (and some of b.1) can be improved by using the longitudinal pipeline in FreeSurfer to analyze your data:
https://surfer.nmr.mgh.harvard.edu/fswiki/LongitudinalProcessing

Finally note, that you cannot trust individual cortical thickness values. People usually run studies comparing 10 or 15 subjects per group. There can always be individual outliers. Manual checking and editing can fix some of that, but not all. Especially effects from motion or hydration cannot be removed.


FreeSurfer & Matlab

(In the directory '$FREESURFER_HOME/matlab' you will find several Matlab-based scripts for reading and writing surface- and volume-based data generated in freesurfer.)

Q. How can I use the Freesurfer Matlab commands if I have a copy of Matlab installed on a Windows machine?

A: You can try to transport files from Linux to Windows which will be tedious and things may not work properly. Alternatively, you could install Matlab for Linux or you can try using Octave (a kind of GNU Matlab for Linux).

Q. Can I load FreeSurfer output in Matlab?

A: If you write it as .mgz format then load_mgh (or MRIread) will read it

Q. How can I make a histogram of cortical thickness?

A: In Matlab you can use the read_curv() function, as in the example below.

  • thick = read_curv ('/usr/local/freesurfer/subjects/bert/surfer/lh.thickness');
    hist (thick,100);


    Notice that there are many zero values that refer to areas where there's no cortical surface. The best procedure is to create a new array without the zero values and then make the histogram.

Q. How can I obtain the thickness of each vertex, and how can you identify which structure each vertex belongs to?

A: You will need to use the read_curv() function to get the thickness estimates, and use the read_annotation() function to assign each index in the surface vector to a specific structure. You can then match up corresponding indices in the the thickness vector and the structure vector to determine the thickness and structure label for each vertex.


FreeSurfer GUI

Q. How can I troubleshoot rendering problems when using TKsurfer?

A: Please refer to TksurferDisplayProblems

UserContributions/FAQ (last edited 2018-01-05 10:52:27 by MorganFogarty)

One thought on “Freesurfer Retinotopic Analysis Essay

Leave a Reply

Your email address will not be published. Required fields are marked *