Measuring gyrification (cortical folding) at any age represents a window into early brain development. Hence, we previously developed an algorithm to measure local gyrification at thousands of points over the hemisphere1. In this paper, we detail the computation of this local gyrification index.
Date Published: 1/02/2012, Issue 59; doi: 10.3791/3417
Keywords: Medicine, Issue 59, neuroimaging, brain, cortical complexity, cortical development
Schaer, M., Cuadra, M. B., Schmansky, N., Fischl, B., Thiran, J. P., Eliez, S. How to Measure Cortical Folding from MR Images: a Step-by-Step Tutorial to Compute Local Gyrification Index. J. Vis. Exp. (59), e3417, doi:10.3791/3417 (2012).
Cortical folding (gyrification) is determined during the first months of life, so that adverse events occurring during this period leave traces that will be identifiable at any age. As recently reviewed by Mangin and colleagues2, several methods exist to quantify different characteristics of gyrification. For instance, sulcal morphometry can be used to measure shape descriptors such as the depth, length or indices of inter-hemispheric asymmetry3. These geometrical properties have the advantage of being easy to interpret. However, sulcal morphometry tightly relies on the accurate identification of a given set of sulci and hence provides a fragmented description of gyrification. A more fine-grained quantification of gyrification can be achieved with curvature-based measurements, where smoothed absolute mean curvature is typically computed at thousands of points over the cortical surface4. The curvature is however not straightforward to comprehend, as it remains unclear if there is any direct relationship between the curvedness and a biologically meaningful correlate such as cortical volume or surface. To address the diverse issues raised by the measurement of cortical folding, we previously developed an algorithm to quantify local gyrification with an exquisite spatial resolution and of simple interpretation. Our method is inspired of the Gyrification Index5, a method originally used in comparative neuroanatomy to evaluate the cortical folding differences across species. In our implementation, which we name local Gyrification Index (lGI1), we measure the amount of cortex buried within the sulcal folds as compared with the amount of visible cortex in circular regions of interest. Given that the cortex grows primarily through radial expansion6, our method was specifically designed to identify early defects of cortical development.
In this article, we detail the computation of local Gyrification Index, which is now freely distributed as a part of the FreeSurfer Software (http://surfer.nmr.mgh.harvard.edu/, Martinos Center for Biomedical Imaging, Massachusetts General Hospital). FreeSurfer provides a set of automated reconstruction tools of the brain's cortical surface from structural MRI data. The cortical surface extracted in the native space of the images with sub-millimeter accuracy is then further used for the creation of an outer surface, which will serve as a basis for the lGI calculation. A circular region of interest is then delineated on the outer surface, and its corresponding region of interest on the cortical surface is identified using a matching algorithm as described in our validation study1. This process is repeatedly iterated with largely overlapping regions of interest, resulting in cortical maps of gyrification for subsequent statistical comparisons (Fig. 1). Of note, another measurement of local gyrification with a similar inspiration was proposed by Toro and colleagues7, where the folding index at each point is computed as the ratio of the cortical area contained in a sphere divided by the area of a disc with the same radius. The two implementations differ in that the one by Toro et al. is based on Euclidian distances and thus considers discontinuous patches of cortical area, whereas ours uses a strict geodesic algorithm and include only the continuous patch of cortical area opening at the brain surface in a circular region of interest.
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.
- 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
- Create the three-dimensional cortical mesh models8,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.
- 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.
- 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.
- 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.
- 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:
Input Patient_Male 20
Input Control_Female 23
- 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
- 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
- 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
- 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 rate10.
- 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
Alternatively, the statistical analyses could eventually be computed at the level of the cortical parcellation integrated in FreeSurfer11. 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 group1,12 or by others13-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 maps1, 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 6mm17, 13mm18-21 and 22mm22,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 adulthood24. In a schematic simplified view, measuring these complementary properties offers the possibility to advance our understanding of the pathogenesis of neurodevelopmental disorders25.
The authors have nothing to disclose.
This research was supported by the National Center of Competence in Research (NCCR) "SYNAPSY - The Synaptic Bases of Mental Diseases" financed by the Swiss National Science Foundation (n° 51AU40_125759). Development of the local Gyrification Index was supported by grants from the Swiss National Research Fund to Dr. Marie Schaer (323500-111165) and to Dr. Stephan Eliez (3200-063135.00/1, 3232-063134.00/1, PP0033-102864 and 32473B-121996) and by the Center for Biomedical Imaging (CIBM) of the Geneva-Lausanne Universities and the EPFL, as well as the foundations Leenaards and Louis-Jeantet. Support for the development of FreeSurfer software was provided in part by the National Center for Research Resources (P41-RR14075, and the NCRR BIRN Morphometric Project BIRN002, U24 RR021382), the National Institute for Biomedical Imaging and Bioengineering (R01 EB001550, R01EB006758), the National Institute for Neurological Disorders and Stroke (R01 NS052585-01) as well as the Mental Illness and Neuroscience Discovery (MIND) Institute, and is part of the National Alliance for Medical Image Computing (NAMIC), funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 EB005149. Additional support was provided by The Autism & Dyslexia Project funded by the Ellison Medical Foundation.