How to Measure Cortical Folding from MR Images: a Step-by-Step Tutorial to Compute Local Gyrification Index

1. Reconstruct the 3D cortical surfaces

This first part of the protocol uses the standard FreeSurfer pipeline as described in the Wiki (http://surfer.nmr.mgh.harvard.edu/fswiki). Note that the commands detailed here describe one way of achieving the cortical surface reconstructions, but equivalent commands may also be used.

1. Import the raw MRI DICOM into FreeSurfer and verify the quality of the image (e.g. that the orientation is correct, the contrast sufficient and the images not moved). This process uses the following commands (replace the text between <...> (inclusive) with values appropriate to a specific instance, and "#" denotes comments):
   
   mksubjdirs # create the folder architecture used by FreeSurfer
   cd /mri # go to the mri folder of your subject
   mri_convert -cm 001.mgz # convert the raw MRI in the # FreeSurfer format
   tkmedit 001.mgz # visualize the converted volume
   
   
2. Create the three-dimensional cortical mesh models^8,9. In order to cope with the problem of buried sulci, FreeSurfer first creates a unitary white matter volume, which is used as a starting point for the initial gray-white surface. This surface is then optimized according to the local gradient of intensity and further expanded to the gray-CSF interface.
   
   recon-all -s # launch the cortical surface reconstruction
   
   At the end of the reconstruction process, you will get two mesh models composed of about 150,000 points for each hemisphere: a white (gray-white interface) and a pial (gray-CSF interface) surface. It is important to note that all the surfaces and volumes remain in the native space, allowing measurement such as volume, surface area, thickness or gyrification index to be measured without deformation.
   
3. Check for the accuracy of these reconstructed surfaces:
   
   tkmedit T1.mgz ?h.pial # the white surface is overlaid in green and the pial surface in red
   
   where ?h denotes the hemisphere: lh.pial for the left hemisphere and rh.pial for the right hemisphere. The Figure 2 (in 2 versions: an animated gif image to be included in the movie and a static one for the website) shows an example of correct white and pial surfaces reconstructions for the "bert" subject distributed along with the FreeSurfer package. If you have to manually correct the result of the reconstruction process, you will find a tutorial on the FreeSurfer Wiki (http://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/WhiteMatterEdits , http://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/PialEdits).

2. Compute the local Gyrification Index

When you are satisfied with your surfaces, compute the local Gyrification Index (lGI) using the command:

recon-all -lgi -s

This command usually runs for about 3 hours for the two hemispheres of one study participant, depending on the power of your workstation. The different steps of lGI process are overviewed in Fig. 1. The computing begins with the creation of an outer surface using morphological closing operation. This outer surface, denoted ?h.pial_outer_smoothed, is further illustrated in Fig. 3. Then, about 800 overlapping circular regions of interest are created on the outer surface. For each one of these regions, a corresponding region of interest is defined on the pial surface. The whole computation ends up with the creation of an individual map containing one lGI value for each point of the cortical surface (i.e. ~150,000 values per hemisphere).

3. Check the result of the lGI calculation for each hemisphere

tksurfer ?h pial -overlay /surf/?h.pial_lgi -fthresh 1

The lGI values are overlaid over the cortical surface. As correct lGI values are typically comprised between 1 and 5, setting the minimum threshold at 1 (with the option fthresh) allows a rapid check: you should not see any gray cortical area. An example of correct individual result is shown in Fig. 4.

4. Statistical group comparisons

The purpose is to quantify the effect of group at each vertex over the cortical surface while controlling for the effect of gender and age. You will need to follow the same process as if you would like to compare cortical thickness at each vertex, but giving ?h.pial_lgi instead of ?h.thickness. Two options are possible to compute the statistical group comparisons: the classical commands are listed first, and the graphical interface (Qdec) is briefly mentioned thereafter.

1. The first option to compare lGI results between groups use the commands listed below; further details can be obtained at https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/GroupAnalysis.
   1. First you will need to create a study specific template giving all your subjects in input:
      
      make_average_subject --subjects ...
      
      The command above will create a subject named "average". Alternatively, you can use the subject "fsaverage" distributed as a part of the FreeSurfer distribution.
   2. Then, create the text file containing the description of the subjects involved in your study (the "FreeSurfer Group Descriptor File"). Your FSGD.txt should look like this:
      
      GroupDescriptorFile 1
      Group Control_Male
      Group Control_Female
      Group Patient_Male
      Group Patient_Female
      Variables Age
      Input Patient_Male 20
      Input Control_Female 23
      <…>
   3. Resample the lGI data in the space of the average subject using the following command for each hemisphere:
      
      mris_preproc --fsgd FSGD.txt --target average --hemi ?h --meas pial_lgi --out ?h.lgi.mgh
      
   4. Smooth the data on the cortical surface to reduce the signal to noise:
      
      mri_surf2surf --hemi ?h --s average --sval ?h.lgi.mgh --fwhm 10 --tval ?h.10.lgi.mgh
      
   5. Compute the group comparison at the level of each vertex. For that you will need to create a contrast text file (e.g. in the case of the FSGD.txt described above, the "contrast.txt" will contain the values "1 1 -1 -1 0" to compute the difference between controls and patients while controlling for age and gender). Finally run the comparison:
      
      mri_glmfit --y ?h.10.lgi.mgh --fsgd FSGD.txt doss --glmdir ?h.lgi.glmdir --surf average ?h --C contrast.txt
      
   6. Visualize the results on your average subject using tksurfer:
      
      tksurfer average ?h inflated
      
      Then load as overlay the sig.mgh file located in the folder ?h.lgi.glmdir/contrast.txt/sig.mgh. Using the option "configure overlay" you can further modify the p threshold as well as correct for multiple comparisons using false discovery rate¹⁰.
      
2. The alternative option for group comparison is to use Qdec, a graphical user interface implemented in FreeSurfer. The use of Qdec with local Gyrification Index implies to pre-smooth the lGI data:
   
   recon-all -qcache -measure pial_lgi -average -s
   
   With Qdec, the FreeSurfer Group Descriptor File is replaced by a slightly different version, the Data Table (qdec.table.dat) that includes the description of the different groups and other confounding variables such as age. A detailed description of the utilization of Qdec is provided at http://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/QdecGroupAnalysis.
   Of note, if the lGI is not available in the list of dependent variables in Qdec, you must add the following line to the .Qdecrc file located in your home directory:
   
   MEASURE1 = pial_lgi

5. Analysis

Alternatively, the statistical analyses could eventually be computed at the level of the cortical parcellation integrated in FreeSurfer¹¹. For that purpose, average lGI values can be extracted for the 34 gyral regions of interest for each hemisphere, and these measurements can be further compared between the different groups. This parcel-wise analysis (as opposed to the vertex-wise analysis described above) could be attractive as it limits the amount of statistical comparisons. However, the lGI at each point quantifies the gyrification in the surrounding circular area, so that the average lGI in a gyral region of interest also reflects to some extent the gyrification in the neighboring regions of interest.

Finally, although the most important issues were described in this protocol, a solution to the other problems that may be encountered during the FreeSurfer or the lGI processing can be found in the archives of the FreeSurfer mailing list (http://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSupport).

6. Representative Results

As described in section 1c of the protocol, you should always carefully check the accuracy of the reconstruction of the cortical surfaces prior to lGI computation. While scrolling between the frontal and the occipital lobe, pay particular attention that vessels and membrane are not included in the pial surface. Check also that the white surface accurately follows the gray-white interface. An example of correct reconstruction is provided in Figure 2 (see the animated gif figure for the whole volume).

At the end of the lGI computation, you will also have to check the result for both hemispheres of each subject. There shouldn't be any cortical area with a lGI result smaller than 1. The section 3 of the protocol and the figure 4 show how to check correct if the output of the lGI computation is correct.

Figure 1. Overview of the lGI computation. First, three-dimensional cortical mesh models are reconstructed from the raw images using the standard FreeSurfer pipeline. These reconstruction algorithms use a binary white matter volume as starting point to overcome the issue of buried sulci. The cortical mesh models typically comprises about 150,000 vertices and are classically used to compute cortical thickness at each point. Similarly, local Gyrification Index (lGI) will be computed at each vertex. For that purpose, an outer surface is created. Then corresponding circular regions of interest are identified on the outer and cortical surface using matching algorithm. After about 800 of generating overlapping regions of interest, the process results in the creation of individual maps of lGI. These maps can be easily interpreted: an index of 5 means that there is 5 time more cortical surface invaginated within the sulci in the surrounding area that the amount of visible cortical surface; an index of 1 means that the cortex is flat in the surrounding area. Finally, statistical group comparisons are computed at the level of each vertex, similarly to cortical thickness comparisons.

Figure 1B. Individual cortical map of lGI. This small movie shows a 360-degrees rotation of an individual cortical lGI map as shown in Fig. 1. It is striking to note that the cortical regions with higher lGI values correspond to the first fold to be created during the in utero life: the sylvian fissure, the superior temporal sulcus and the intraparietal sulcus on the lateral view of the brain, and the parieto-occipital sulcus on the medial view of the brain. View movie

Figure 2. Example of adequate cortical surface reconstruction (one coronal section). After the end of the reconstruction process, the cortical surfaces should be accurately verified across the entire cerebral volume. The inner cortical surface (denoted white surface, in green on the image) should precisely follow the gray-white interface. The outer cortical surface (i.e. gray-CSF interface, denoted pial surface, here in red) should not include any piece of vessel or membrane. Of note, the example presented here uses the "bert" subject distributed along with the FreeSurfer package.

Figure 2B. Example of adequate cortical surface reconstruction (full volume). This animated gif image shows the cortical surface of the left hemisphere of the "bert" subject on each coronal section, as seen by scrolling from the most frontal to the most occipital coronal sections with FreeSurfer. View movie

Figure 3. Example of outer surface computed as a part of the lGI process (one coronal section). The first step in the lGI computation is the creation of an outer surface enveloping the hemisphere. This surface (denoted ?h.pial_outer_smoothed in FreeSurfer) can be checked using tkmedit. Here, the "bert" subject distributed with FreeSurfer is used as an example.

Figure 3B. Example of outer surface computed as a part of the lGI process (full volume). This animated gif image shows the outer surface of the left hemisphere on each coronal section, as seen by scrolling from the most frontal to the most occipital coronal sections with tkmedit in FreeSurfer. View movie

Figure 4. Example of correct lGI output as viewed with FreeSurfer. Different orientations of the cortical surface of the "bert" subject with lGI values overlaid. The color-code is the default "heat" overlay as seen with tksurfer in FreeSurfer. Using a minimum threshold of 1, all the vertices must be colored and no cortical area should appear in gray. Of note, the color overlay can be modified using the option "Configure Overlay" in tksurfer, where the minimum and maximum values, as well as the histogram of the overall distribution of the lGI can also be checked.

The protocol above describes how to measure local Gyrification Index based on cerebral T1-weighted MRI and conduct statistical group comparisons. Our method has been specifically designed to localize early disruption in the cortical expansion process and as such is of particular interest in many neurodevelopmental or psychiatric conditions. Examples of group comparisons in clinical samples can be found in publications by our group^1,12 or by others^13-16. The process is fully automated and requires only command to be executed, although two parameters can be modified.

The first modifiable parameter is at the level of lGI computation: the radius of the circular region of interest. The radius is by default set to 25 mm, which was chosen in order to include more than one sulcus at a time while retaining an adequate resolution. Our validation paper includes an experiment of the effect of the radius on the cortical gyrification maps¹, showing that large radii tend to smooth the cortical maps with a dilution of the local maxima. For clinical studies, we would recommend a radius between 20 and 25 mm.

The second adjustable parameter is the amount of smoothing at the level of the statistical analyses. In order to increase signal-to-noise ratio, data are smoothed on the cortical mesh employing an iterative nearest-neighbor averaging procedure. Cortical thickness studies with similar conditions as ours (i.e. data measured in the native space, same data distribution over the cortical surface, and same sulcal-based registration technique with a study specific template) commonly use a full-width at half-maximum (FWHM) of 10 mm (For reference, cortical thickness studies co-authored by the developers of FreeSurfer used a kernel of 6mm¹⁷, 13mm^18-21 and 22mm^22,23). In the above protocol, we propose to use a FWHM of 10 mm to keep in line with most of the cortical thickness literature. However, as the cortical lGI maps are already relatively smooth, the result of the comparisons will barely change according to the presence or absence of smoothing.

Although the implementation of the methods for measuring gyrification and thickness share common steps, we would like to emphasize that both measures reflect different properties of the cortical morphology. As already stressed above, gyrification is mostly influenced by early events. Contrarily, cortical thickness is largely sensitive to maturational changes during childhood, adolescence and early adulthood²⁴. In a schematic simplified view, measuring these complementary properties offers the possibility to advance our understanding of the pathogenesis of neurodevelopmental disorders²⁵.