In vivo measurements of lamination patterns in the human cortex.

The laminar composition of the cerebral cortex is tightly connected to the development and connectivity of the brain, as well as to function and pathology. Although most of the research on the cortical layers is done with the aid of ex vivo histology, there have been recent attempts to use magnetic resonance imaging (MRI) with potential in vivo applications. However, the high-resolution MRI technology and protocols required for such studies are neither common nor practical. In this article, we present a clinically feasible method for assessing the laminar properties of the human cortex using standard pulse sequence available on any common MRI scanner. Using a series of low-resolution inversion recovery (IR) MRI scans allows us to calculate multiple T1 relaxation time constants for each voxel. Based on the whole-brain T1 -distribution, we identify six different gray matter T1 populations and their variation across the cortex. Based on this, we show age-related differences in these population and demonstrate that this method is able to capture the difference in laminar composition across varying brain areas. We also provide comparison to ex vivo high-resolution MRI scans. We show that this method is feasible for the estimation of layer variability across large population cohorts, which can lead to research into the links between the cortical layers and function, behavior and pathologies that was heretofore unexplorable.

INTRODUCTION

In vivo measurement of the cortical layers with magnetic resonance imaging (MRI) in humans has been a particular challenge in recent years (Trampel, Bazin, Pine, & Weiskopf, 2017). The different lamination patterns of the cortex are connected to the function of different brain areas (Brodmann & Garey, 1999; Von Economo & Koskinas, 1925). Notably, the granularity of the cortex—that is, the relative volume of the inner granular layer (Layer IV)—is related to the type of functionality, as sensory cortices are more heavily granular, while the motor cortex, for example, is agranular (Beul & Hilgetag, 2015). Measuring the lamination patterns of the cortex is a crucial link between microstructure and function (Glasser & Van Essen, 2011; Lifshits et al., 2018; Trampel et al., 2017). However, in vivo MRI measurements of the cortical layers present several major challenges (Edwards, Kirilina, Mohammadi, & Weiskopf, 2018; Trampel et al., 2017). Currently, most common MRI scanners designated for human brain imaging have resolution limits that prevent accurate lamination measurements in practical acquisition times (Lifshits et al., 2018; Shamir et al., 2019; Trampel et al., 2017). Moreover, the curvature of the human cortex further exacerbates the partial volume effect (PVE), which limits the isolation of each unique layer (Leprince et al., 2015; Lifshits et al., 2018; Shamir et al., 2019).

Despite the challenges, there are various attempts to measure cortical lamination patterns with MRI, both in vivo and ex vivo. These are usually focus on brain areas where there are certain distinct layers, such as the heavily myelinated Layer IV in the primary visual cortex (Barbier et al., 2002; Clark, Courchesne, & Grafe, 1992; Geyer, Weiss, Reimann, Lohmann, & Turner, 2011). It has been shown that T1‐weighted contrasts can distinguish between different cortical microstructures, as T1 is affected by various microstructural properties, such as myelination, iron content, and morphology (Clark et al., 1992; Stüber et al., 2014; Barazany and Assaf, 2012). Furthermore, Barazany and Assaf (2012) showed that it is possible to identify several different clusters of gray matter voxels in the brain based on quantitative T1‐mapping. Several other methods, such as T2‐weighted and diffusion‐weighted MRI, have also been used to visualize the layers, with varying degrees of success (Assaf, 2018; Edwards et al., 2018; Trampel et al., 2017).

Recently, Lifshits et al. (2018) attempted to overcome the resolution issues by estimating the PVEs of different tissues based on their T1 by performing multicomponent analysis that infers subvoxel information. Acquiring multiple inversion recovery (IR) scans, they have shown that subvoxel T1 modeling can be used to identify different populations of cortical tissue (Lifshits et al., 2018; Shamir et al., 2019).

In this article, we demonstrate the effectiveness and validity of such a method and show that it can be used to identify microstructural properties of the cortex beyond what is possible with standard/high‐resolution imaging techniques, and, as a demonstration of the method, show how these properties correlate with age. We also demonstrate that how cortical parcellation is possible with this method, as well as examine how this method compares to high‐resolution ex vivo imaging analysis. We suggest that this method can lead to new avenues of research into cortical parcellation and the link between microstructure and function. In addition, this approach can be used in whole brain, large population studies which allows, for the first time, to explore the variance of layer composition over wide populations. We anticipate that despite the low‐resolution imaging protocol used, this method will allow for measuring the different microstructural patterns of the cortex, reveal correlations between cortical microstructure and age.

METHODS

Subjects

Two‐hundred subjects were recruited for this study (98 females and 102 males, aged 30.09 ± 11.82 years; see Figure S1 for subject age and sex distribution). Subjects were neurologically and radiologically healthy, with no history of neurological diseases, and normal appearance of clinical MRI protocol. The imaging protocol was approved by the institutional review boards of Sheba Medical Centers and Tel Aviv University, where the MRI investigations were performed. All subjects provided signed informed consent before enrollment in the study.

MRI acquisition

All subjects were scanned at the Alfredo Federico Strauss Center for Computational Neuroimaging at Tel Aviv University, with a 3 T Siemens MAGNETOM Prisma MRI scanner (Siemens Medical Solutions, Erlangen, Germany) with a 64‐channel RF head coil. The protocol consisted of a series of 44 IR prepared spin echo echo‐planar‐images (IRSE‐EPI) with inversion times (TI) between 50 and 2,500 ms. The scans were acquired at a 3 × 3 × 3 mm3 resolution in the axial plane, covering the entire brain, with the following parameters: TR/TE = 12,000/30 ms, GRAPPA factor of 2 with a matrix size of 68 × 68 and 42 slices (no gap). The subjects were also scanned with an MPRAGE sequence (TR/TE/TI = 2400/2.78/1000 ms) at a 1 × 1 × 1 mm3 resolution, used as an anatomical reference with a good contrast between gray and white matter. The duration of the protocol was roughly 15 min. All subjects were scanned as part of the Tel Aviv University (TAU) Brain Brank Initiative, and the imaging protocol included additional sequences that were not used in this study.

MRI data analysis

The series of IR EPI scans were registered to the first IR scan using MATLAB's (Mathworks, Natick, MA) Image Processing Toolbox and SPM 12 (Penny, Friston, Ashburner, Kiebel, & Nichols, 2006). Following that, the IR data were fitted using nonlinear least‐squares optimization on a voxel‐by‐voxel basis to the multicomponent IR function (Lifshits et al., 2018) using an in house MATLAB code: where and are the inversion time and the magnetization of the ith IR image, respectively, is the voxel's magnetization at and is proportional to proton density, is the longitudinal relaxation time for each component with being the component's volume fraction. is the number of components fitted for a specific voxel. This method generates a multicomponent T1‐map for each subject, with a matching partial volume map for each subvoxel component. The number of computed components for each voxel was determined by fitting a range of possible components and comparing the root mean squared error (RMSE) of their fit.

In order to identify different tissue types in the brain and specifically different distributions of gray matter tissue, we fitted a mixture model to each subject's whole‐brain T1‐distribution (calculated from the subject's multicomponent T1‐map). We used a mixture of Student's t‐distributions model (Peel & McLachlan, 2000) instead of the more common Gaussian mixture model, as it is more robust to noise and outliers in the data (McLachlan & Peel, 2004; Peel & McLachlan, 2000). For each subject, we attempted to fit up to 25 mixtures, where the optimal number was determined by Bayesian information criteria (BIC; Schwarz, 1978) of each subject's distribution to the mean model across subjects. We found that on average, the 18‐mixture fit had the best fit (Figure S2), and thus it was used for the following analysis steps.

We then identified tissue distributions whose mean T1‐value was in the expected gray matter T1 regime (Figure S2), which is roughly between 800 ms and 1800 ms (Clare & Jezzard, 2001; Deoni, 2007; Ethofer et al., 2003; Gelman, Ewing, Gorell, Spickler, & Solomon, 2001; Lu et al., 2005; Stanisz et al., 2005; Wansapura, Holland, Dunn, & Ball, 1999; Weiskopf et al., 2013; Wright et al., 2008). For each of these gray matter mixtures (T1 classes), we calculated the probability of the mixture in each voxel using Bayes' formula: where represents a mixture, is the general whole‐brain probability of a ‐value, is the mixture probability, and is the probability of the ‐value in the mixture .

These tissue probability maps, which are calculated in the IR scan space, were then registered and normalized to the MNI152 space (Evans et al., 1993; Fonov, et al., 2011; Fonov, Evans, McKinstry, Almli, & Collins, 2009) using SPM 12 (Penny et al., 2006) as following: first registering the shortest TI IR scan to the subject's MPRAGE, then registering the MPRAGE to the MNI152 space, and finally applying the transformation to the probability maps. We then calculated the relative laminar proportions for each parcel in three different brain atlases: the 100 and 1,000 parcels version of the 17‐network atlas from Schaefer et al.'s (2018) local–global parcellation; and the 214 cortical areas from Julich‐Brain's maximum probability map (Amunts, Mohlberg, Bludau, & Zilles, 2020). These resulted in a p‐by‐k matrix for each subject (where p is the number of parcels in the atlas and k is the number of gray matter mixtures), in which each row is the T1‐class probability vector and calculated the mean probability vectors across the subjects. Furthermore, we calculated the cortical thickness estimation for the 100 parcels version of Schaefer et al.'s (2018) atlas using CAT12 (Gaser & Dahnke, 2016; Luders et al., 2006). Then, for each subject, we multiplied each area's T1‐class probability vector by the area's cortical thickness, resulting in T1‐class thickness vectors, which were also averaged across all subjects.

Based on these mean probabilities, we used the fuzzy C‐means (FCM) method (Bezdek, 1981; Dunn, 1973) to cluster the 1,000 parcels of the local–global parcellation into six fuzzy clusters, mimicking the six cortical types identified by Beul and Hilgetag (2015). As this version of Schaefer's atlas is composed of 1,000 small areas, it serves as a practical basis for further clustering, and as it is a function‐based atlas, it inserts less structural bias in our analysis than an anatomy‐based atlas would.

Using the T1‐class thickness vectors of the 100‐parcels Schaefer et al. (2018) atlas, we calculated the correlations between age and the thickness of each T1‐class of each area and tested for significance following false discovery rate (FDR) correction (Benjamini & Hochberg, 1995). We used the 100‐parcel atlas as it is composed of larger areas that are consistent across subjects.

Finally, we calculated the class covariance matrices within and between T1‐classes, as well as cortical volume covariance matrices, based on the T1‐class probabilities of the areas in the Julich‐Brain atlas. The latter was done using CAT12 (Gaser & Dahnke, 2016). This process was done in order to assess whole‐brain interclass and intraclass correlations, whether there is interhemispheric correlation within and between classes. The Julich‐Brain atlas was used as it is a symmetrical cytoarchitectonic atlas with homologous areas in each hemisphere, allowing for analysis of interhemispheric correlations.

Ex vivo acquisition and analysis

From a fixed human brain (age: 40 years), auditory (A1), motor (M1), somatosensory (S1) and visual (V1) cortices were dissected. The dissected sections were preserved in 4% paraformaldehyde in polybutylene succinate (PBS). Twenty‐four hours prior to MRI, the sections were soaked in PBS at room temperature and then placed in Fluorinert (FC‐770; 3M, St. Paul, MN, USA) tube for scanning. MRI was performed on a 7 T/30 Bruker Biospec (Bruker, Germany) at the Alfredo Federico Strauss Center for Computational Neuroimaging at Tel Aviv University. The protocol consisted of seven 3D IR rapid acquisition scans with relaxation enhancement (IR‐RARE) with the following TIs: 250, 350, 380, 420, 450, 480, and 650 ms at a resolution of 150 × 150 × 225 μm3 (TR/TE = 5000/8 ms, RARE factor = 8). The TIs were selected based on Barazany and Assaf's (2012) protocol, with the higher TIs discarded as there is no cerebrospinal fluid with higher T1‐value to account for. Only seven TIs were used as due to the high resolution of the scans, modeling the PVE was less relevant.

For each region, the series of IR images were coregistered and resliced using SPM8 (UCL, London), and the gray matter was manually segmented. We then calculated quantitative T1‐maps according to the inversion recovery spin echo (IR‐SE) equation (Barral et al., 2010): Following that, we fitted a mixture of Student's t‐distributions with six components to the T1‐maps, to account for six cortical layers. We then calculated the relative ratio of each to assess the relative volume of the T1‐layer it represents, and performed hard clustering on the gray matter voxels based on the mixture model. These results were compared qualitatively compared to the in vivo T1‐layers extracted in these specific areas, as described in the previous section, by calculating the mean laminar differences between the two measurements.

RESULTS

In this article, we explore the microstructural properties of the human cortex as measured by T1‐relaxation time. We used multicomponent T1 analysis to identify different T1 populations of gray matter and found six populations in the T1 gray matter regime (ranging 800–1800 ms). It should be noted that we did not predetermine the number of populations in this range, and that they are by no means a direct measures of the histological layers (see Section 4). We then calculated the correlations between these populations and age and examined the relationship between the different populations across the brain. We also performed ex vivo analysis using our method on four different excised cortical sections.

T‐layers analysis

The mean probabilities for each of these gray matter T1‐class across 200 subjects can be seen in Figure 1 (for SDs, see Figure S2). It can be seen, for example, that Class 5 has a higher probability in tempo‐occipital areas, while Class 6 has a higher probability in visual areas in the occipital cortex. Class 3 has a higher probability in inferior parts of the cortex in the general area of the hippocampal formation. Class 4 has a higher probability in parietal and temporal areas. Classes 1 and 2 have a higher probability in frontal and parietal regions, compared to the rest of the cortex. It should be noted while the gray matter classes were identified based solely on their T1, they do not necessarily correspond to the cytoarchitectonic cortical layers, and that there was no predetermined specific number of such mixtures that we were looking for.

FIGURE 1

Mean probabilities of the T1‐classes for the 1,000 areas in Schaefer et al.'s (2018) 17‐network parcellation of the human cortex, across the 200 subjects. All the data are presented on the MNI152 brain (Evans et al., 1993; Fonov et al., 2009; Fonov, et al., 2011). For SDs, see Figure S3. The same data with dynamic color range for each T1‐classes are presented in Figure S4. For these results presented on the FsAverage surface, see Figure S8

Six fuzzy clusters were identified using the FCM algorithm (Figure 2). It can be seen that Cluster 1 includes mostly temporal and parietal regions, while Cluster 2 includes areas in all lobes. Cluster 3 is located almost exclusively in the visual cortex, and Cluster 4 similarly is located in the inferior parts of the cortex around the uncus and the orbital gyri. Cluster 5 is dispersed heavily in upper frontal and parietal areas, and Cluster 6 appears to include mostly primary and supplementary motor areas.

FIGURE 2

Maximum fuzzy clusters using the T1‐class probabilities for the areas in Schaefer et al.'s (2018) 17‐network, 1,000‐areas parcellation of the human cortex. All the data are presented on the MNI152 brain (Evans et al., 1993; Fonov et al., 2009; Fonov, et al., 2011). For these results presented on the FsAverage surface, see Figure S10

Correlation with age

Correlations between age and T1‐class thickness measurements, as well as between age and cortical thickness, were calculated and corrected for multiple comparisons (Figure 3). The results show major negative correlations between age and T1‐classes 4 and 5 in most of the cortex, as well as positive correlations between age and T1‐class 1. Some minor correlations can be seen between the other T1‐classes and age. Cortical thickness is shown to be negatively correlated with age across most of the cortex.

FIGURE 3

Significant correlations between T1‐class thickness and age across the six T1‐classes, as well as correlations between age and overall cortical thickness (top row). All the data is presented on the FreeSurfer FsAverage brain surface. For these results presented on the MNI152 brain, see Figure S11

Cross‐population analysis

T1‐classes covariance matrices (n = 200) show positive correlations and high interhemispheric correlation between homologous areas (Figure 4a), which can also be seen in the top right and bottom left of the matrices. It can also be seen that there are positive interclass correlations between T1‐classes 1, 2, and 3 (Figure 4a), as well as between T1‐classes 4 and 5 (Figure 4b), suggesting that the variability in T1‐class composition across the cortex has a certain pattern, as opposed to simply being random. The interhemispheric patterns are consistent in these interclass matrices as well.

FIGURE 4

Correlation matrices. (a) Correlation and covariance matrices between the T1‐classes,, where correlations between homologous areas can be seen; positive intra‐class correlations can be seen, especially between homologous areas in both hemispheres, as well as similar interclass correlations between T1‐classes 1, 2 and 3, between T1‐classes 3 and 4, and between T1‐classes 4 and 5; negative interclass correlations between T1‐class 6 and all other classes, and weaker negative correlation between T1‐classes 1 and 5; (b) Cortical gray matter volume covariance matrix. In all matrices, the upper half of the y‐axis and left half of the x‐axis represent left hemisphere areas, and the lower half of the y‐axis and right half of the x‐axis represent the right hemisphere

Ex vivo comparison

Ex‐vivo analysis found T1‐layer patterns similar to those found in our in vivo method (Figure 5) in the four cortical sections examined— the auditory (A1), motor (M1), somatosensory (S1), and visual (V1) primary cortices. These similarities mean layer ratios difference between in vivo and ex vivo measurements were similar across all four sections (A1: 0.04 ± 0.04; M1: 0.04 ± 0.03; S1: 0.06 ± 0.05; V1: 0.04 ± 0.03). In areas A1 and V1, it can be seen that the sixth T1‐layer is notably the thickest in both in vivo and ex vivo measurements. It can also be seen that T1‐layers 1, 2, and 3 are relatively smaller both in the in vivo and ex vivo analysis in all four areas. These similarities are qualitative only as no statistical analysis could be performed on single samples, yet they provide some high‐resolution validation for our methods.

FIGURE 5

T1‐lamination patterns in four ex vivo cortical sections: the primary (a) auditory, (b) motor, (c) somatosensory and (d) visual cortices, showing comparison between in vivo and ex vivo measurements (left), a quantitative T1 image of a slice of the ex vivo section (center), and the same slice overlaid by the six different gray matter T1‐clusters

DISCUSSION

In this article, we present a method for in vivo measurement of whole‐brain cortical lamination. This method is based on multicomponent T1 analysis of inversion recovery imaging and identification of multiple gray matter clusters in the human cortex, which allows us to overcome the common resolution and PVE limitations of conventional MRI protocols, which hinder measurements of cortical microstructure. We show that T1‐MRI properties provide a unique measure of the cortex, and that these measures of cortical microstructure can be used to distinguish between different cortical areas. We also show how it can be used to perform research on various groups and populations. As the cortical layers are the fingerprint of human function (Barazany and Assaf, 2012), there is a critical need for tools that enable exploration and understanding of this phenomenon in large populations. Such studies may also help unravel the laminar patterns' role in neurodegenerative diseases and other pathologies.

Much of the research into in vivo imaging of cortical layers has focused on high‐resolution imaging of small parts of the brain, usually requiring a high magnetic field or long scanning times (Barbier et al., 2002; Clark et al., 1992; Geyer et al., 2011; Lifshits et al., 2018). These methods' applicability to large population studies and whole‐brain analysis is therefore limited. The approach proposed in this article, however, is to use fast, low‐resolution T1 imaging to model the PVE rather than to resolve it with a high enough resolution.

By performing subvoxel analysis on multiple IR datasets with varying contrasts, we are able to characterize the microstructure of the cortex voxel without using the high‐resolution techniques required to actually image it. The main advantages of this approach are its simple protocol, which is applicable on any MRI scanner, its short acquisition time (roughly 15 min), and its ability to perform this analysis on a whole‐brain basis. We show here that this combination makes our method suitable to general population studies and comparisons with large sample sizes.

As we demonstrated, our T1‐layers measurements can capture some of the known microstructural differences between various cortical regions. An example of this is the unique features of the visual cortex or the paleocortex (see Figures 1 and 2). While these are not exact measurements of the cortical layers as they are defined by histology, comparing them to high‐resolution imaging of ex vivo sections shows a relatively good correspondence. It should be stressed that though we found six distinct populations of gray matter T1 tissues, they do not necessarily correspond fully to the cytoarchitectonic layers. However, it is not unreasonable to assume that differences in cell arrangement density, and hence in tissue density, will be reflected in the T1 properties of the tissue (Fatterpekar et al., 2002; Lifshits et al., 2018). Some of the results of our clustering are also in line with histological lamination measurements (Beul & Hilgetag, 2015). Less distinct, yet still valid comparisons can also be seen in the frontal and temporal cortices, which provides some validity for our method.

Assessing the covariance matrices provides support for the reliability of our measurements, as well as additional validation. The high interhemispheric correlations between homologous areas show the robustness of our method. Meanwhile, interlayer correlations may be related to the intrinsic connections between the layers, though further research is needed for conclusive statements on the matter. In addition, we have shown that the negative correlation between cortical thickness and age is expressed mostly in some of the T1‐classes and not in others. We also show that some of the classes actually show somewhat of an increase in thickness. While the underlying cause for this is difficult to identify, similar phenomena in lamination patters had been previously observed with histological methods (Zhang, Hua, Zhu, & Luo, 2006). Another possible explanation for this is myelination of some cortices that occurs with age and is related to the thinning of the cortex (Natu et al., 2019). Finally, high‐resolution ex vivo comparisons show similarity to our method in areas that are considered both granular, such as the primary visual and somatosensory cortices, as well as to agranular cortices such as the primary motor cortex (Beul & Hilgetag, 2015).

It should be noted that our approach is not without its limitations. First and most obviously, we do not directly visualize the cortical layers themselves, but only measure microstructural features indirectly in low resolution. For this reason, we use the terms “T1‐layers” or “T1‐classes” instead, to differentiate from histological terms. The other major limitation of this method is that the low‐resolution imaging technique used prevents us from capturing the variation in laminar thickness where cortical curvature affects it. The cortical layers vary not only between cortical areas, but also within areas according to the folding of the cortex (Waehnert et al., 2014), variations that require a significant increase in image resolution. This would require not only a much longer scanning protocol but also more advanced MRI technology than is commonly available for in vivo acquisitions. Therefore, our method effectively only measure the average microstructural composition within small areas. Our method is also sensitive to noise and outliers in the data, which may be mitigated by increasing time resolution of our protocol (i.e., the number of different IR scans), which will, again, come at the expense of the short, convenient nature of the current protocol used. With regards to the study cohort, it can be seen that it trends toward younger subjects, mostly due to the recruiting logistics of the TAU Brain Bank Initiatives, from which this study draws its participants.

This method and its relative simplicity may enable population‐wide, whole‐brain research into cortical lamination patterns that heretofore was only studied on a much smaller scale and prespecified regions. It has already been shown that similar lamination measurements identify unique characteristics in stroke (Lotan et al., 2019) and cortical dysplasia (Lotan et al., 2021) patients. Further use of our method could allow for further investigation of the link between cortical microstructure and various functions and pathologies more extensively than has yet been done with brain imaging.

Appendix S1: Supporting Information

Click here for additional data file.