Contents lists available at ScienceDirect NeuroImage journal homepage: www.elsevier.com/locate/neuroimage Resolution considerations in imaging of the cortical layers Shlomi Lifshitsa , Omri Tomerb , Ittai Shamirc , Daniel Barazanyd , Galia Tsarfatye , Saharon Rosseta , Yaniv Assafb,f,⁎ a Department of Statistics and Operations Research, Faculty of Exact Sciences, Tel-Aviv University, Tel-Aviv, Israel b Sagol School of Neuroscience, Tel Aviv University, Tel Aviv, Israel c Department of mechanical engineering, Faculty of Engineering, Tel Aviv University, Tel Aviv, Israel d The Strauss Center for Computational Neuroimaging, Tel Aviv University, Tel Aviv, Israel e Department of Diagnostic Imaging, Sheba Medical Center, Tel Hashomer, Ramat Gan, Israel f Department of Neurobiology, Faculty of Life Sciences, Tel Aviv University, Tel Aviv, Israel A R T I C L E I N F O Keywords: Cortical layers T1 relaxation Partial volume effect MRI resolution Brain parcellation A B S T R A C T The cortical layers are a finger print of brain development, function, connectivity and pathology. Obviously, the formation of the layers and their composition is essential to cognition and behavior. The layers were traditionally measured by histological means but recent studies utilizing MRI suggested that T1 relaxation imaging consist of enough contrast to separate the layers. Indeed extreme resolution, post mortem, studies demonstrated this phenomenon. Yet, one of the limiting factors of using T1 MRI to visualize the layers in neuroimaging research is partial volume effect. This happen when the image resolution is not high enough and two or more layers resides within the same voxel. In this paper we demonstrate that due to the physical small thickness of the layers it is highly unlikely that high resolution imaging could resolve the layers. By contrast, we suggest that low resolution multi T1 mapping conjugate with composition analysis could provide practical means for measuring the T1 layers. We suggest an acquisition platform that is clinically feasible and could quantify measures of the layers. The key feature of the suggested platform is that separation of the layers is better achieved in the T1 relaxation domain rather than in the spatial image domain. Introduction The cortical layers are assumed to play an important role in brain function, development and pathology. The layout of the layers not only defines the anatomical location of different brain regions but also affects their functions (Kandel et al., 2000). The layers are formed through brain development while neurons migrate to form the cortex (Kandel et al., 2000). Through that migration process, cortical connections are formed, hence the formation of the layers also affects brain connectivity. The layers are defined by the arrangement and density of neurons (cyto-architecture) and myelin (myelo-architecture) in the cortex (Brodmann and Garey, 1999). With that definition, the common quantitative measure of the layers is their thickness and its variation across the cortex (Brodmann and Garey, 1999; Vogt, 1911; von Economo and Koskinas, 1925). Yet the cyto- and myelo-architectonic layer thickness can be measured only with histological measures, which were used to demonstrate that the layout of the layers (both composition and thickness) overlaps with the functional representation of cortical regions. It was also demonstrated that abnormal layer composition leads to functional deficits and cognitive impairments (Hof et al., 1996; Kalus et al., 1997; Kordower et al., 2001; Masliah et al., 1991). Despite the assumed great importance of the layers, as no state-ofthe-art in-vivo imaging bio-markers of the layers exist, large population studies on the layers role in cognition, behavior, brain physiology and neurodegeneration are limited. While fMRI and Voxel-based-morphometry (VBM) provide information on cortical function and thickness across groups of subjects (typically where the cortex width is represented by 2–4 voxels), none of these methods are measured in a resolution that allows visualization of the layers. Recently, meso-scale resolution MRI of brain function or anatomy, typically using high field MRI, allowed rough localization (of various parameters from BOLD response to susceptibility and diffusion properties) within the cortical stripe from inner (close to the WM) to outer (close to the CSF) parts of the cortex (Barazany and Assaf, 2011; Barbier et al., 2002; Duyn et al., 2007; Koopmans et al., 2010; Leuze et al., 2012; Olman et al., 2012). Of specific interest was the use of T1 contrast to estimate myelin content along the cortical surface as a neuro-anatomical and functional indicator for differences between brain regions and cortical sub- http://dx.doi.org/10.1016/j.neuroimage.2017.02.086 Accepted 27 February 2017 ⁎ Corresponding author at: Department of Neurobiology, Faculty of Life Science, Tel Aviv University, Tel Aviv, Israel. E-mail address: assafyan@gmail.com (Y. Assaf). NeuroImage 164 (2018) 112–120 Available online 06 March 2017 1053-8119/ © 2017 Elsevier Inc. All rights reserved. MARK structures (Dinse et al., 2015; Fracasso et al., 2016; Goncalves et al., 2015; Lutti et al., 2014; Sereno et al., 2013). The understanding that investigation of the layers might reveal more detailed information on brain function triggered, from the early days of MRI, research aimed to find which MRI contrast demonstrates best the layers (Barbier et al., 2002; Clark et al., 1992; Eickhoff et al., 2005; Walters et al., 2003). These kinds of studies focused on excised tissue where extreme image resolution can be obtained at hours of scanning and on striate cortex where the heavily myelinated layer 4 (stripe of Gennari) can be easily detected (Barazany and Assaf, 2011; Barbier et al., 2002; Clare and Bridge, 2005; Turner et al., 2008). It was found that T1 and T2 weighted contrasts enable demonstration of the layers probably due to differences in the myelin content across the layers. Over the years T1 contrast was preferred over T2 as the differences in T1 spread over a much large scale and the contrast in T1 (especially with inversion recovery (IR) sequence) is much higher. Despite this promising observation, one of the challenges in cortical layer imaging is to go beyond the level of the striate cortex and provide a robust framework that will enable characterization of all the layers for the entire brain and in-vivo (Barazany and Assaf, 2011). Additional studies assumed that T1 values within the cortex are correlated with the layers myelin content (Dinse et al., 2015; Fracasso et al., 2016; Lutti et al., 2014; Sereno et al., 2013). These studies showed that high resolution quantitative T1 MRI can provide additional insight into cortical anatomy and provide new means to investigate the cortex. Yet, it should be noted that even the highest resolution MRI can achieve (submillimeter) does not allow visualization of the layers but rather estimation of various gross factors related to the layers (e.g. myelination content). In recent years, with the advances in high magnetic field scanners, it was speculated that in-vivo sub-millimeter imaging of the human brain will enable cortical layer characterization. Several studies used T1 contrast to explore the cortical layers providing additional proof that T1 contrast is highly sensitive to the layer composition, implying that even higher resolution will enable better characterization of the layers (Barazany and Assaf, 2011; Geyer et al., 2011; Turner et al., 2008). In this paper we argue that extreme resolution MRI will not allow adequate and reasonable characterization of the layers. We note that as the cortex is only 2–4 mm thick (and in many regions some layers’ width is even less than 200 μm), the minimum MRI resolution to adequately resolve the layers should be in the order of tens of microns for the entire brain. Such acquisition resolution cannot be currently achieved and therefore partial volume effects (PVE) is the main obstacle for comprehensive and robust cortical layer imaging. In this paper we demonstrate an alternative approach to resolve PVE in cortical layer imaging using composition analysis of multi T1 components within a voxel. We show, across species, that in vivo imaging of the layers can be achieved with good separation between the layers, for the whole brain and with reasonable scanning time. We suggest that separation between the layers cannot be achieved by means of increased resolution, but rather by composition analysis in the T1 time domain. Methods In this work we have compared 4 experiments: high vs. low resolution on human and rat samples. A comparison between all experimental conditions is given in details in Table 1 and in the text below. Human experiments Subjects 33 healthy subjects (age range 25–35) recruited for this study, 15 (8 males, 7 females) of them were scanned in the high-resolution experiment (see details below) and 18 (10 males, 8 females) for the low-resolution experiment (see details below). Subjects were neurologically and radiologically healthy with no history of neurological diseases. The imaging protocol was approved by the institutional review boards of Tel Aviv Sourasky Medical Center, Sheba Medical Center and Tel Aviv University. All subject signed informed consent before enrollment in the study. Acquisition We have used 2 different protocols to demonstrate the ability of T1 weighted MRI to characterize the layers: 1. A high resolution protocol that included a series of inversion recovery prepared fast spin echo images acquired at resolution of 0.43×0.43×1.5 mm3 (matrix size of 512×384 reconstructed to 512×512) covering the left hemisphere in the sagittal plane (data taken from:(Barazany and Assaf, 2011)). These experiments were scanned on a 3T Signa general electric scanner (GE Healthcare, Milwaukee, WI) with an 8-channel RF coil. TR/TE were 11,000/ 14 ms and 7 inversion times (TIs) (230, 432, 575, 760, 920, 1080 and 1380 ms). The acquisition time for the inversion recovery protocol lasted 45 min. In addition we have scanned a spoiled gradient recalled echo (SPGR) (TR/TE/TI/Flip angle=9.3 ms/ 3.8 ms /450 ms/13°) at resolution of 1x1×1 mm3 as an anatomical reference with high gray/white matter contrast. 2. A low resolution protocol that included a series of inversion recovery prepared echo planar images acquired at resolution of 3x3×3 mm3 in the axial plane covering the entire brain. These experiments were scanned on a 3T Magnetom Siemens Prisma (Siemens Medical Solutions, Erlangen, Germany) scanner with a 64-channel RF coil and the following parameters: TR/TE=10000/30 ms and 107 inversion times spread between 50 ms up to 3000 ms), GRAPPA factor of 2 with matrix size of 68×68. The acquisition time for the inversion recovery data set was approximately 20 min. In addition we have scanned MPRAGE sequence (TR/TE=1750/2.6 ms, TI=900 ms) at resolution of 1x1×1 mm3 as an anatomical reference with high gray/ white matter contrast. Rat experiments An excised male rat brain aged 20 weeks (Wistar, Harlan labs), formalin fixated was scanned on a 7T/30 Bruker Biospec scanner (Bruker, Karlsruhe, Germany) with a transmit volume coil and quadrature surface coil as receiver. As in the human experiment we used two protocols: 1. A high resolution a 3D inversion recovery prepared fast spin echo with the following parameters: TR=4000 ms, TE=5.92 ms (with 16 echoes train length resulting in effective TE of 29 ms), TIs varied from 100 ms to 3000 ms in 18 steps, 2 averages with cubic image resolution of 110 μm3 . The acquisition time for each TI was 42 min and the total scanning time was about 16 h. 2. Low-resolution 3D inversion recovery prepared fast spin echo with 60 inversion times spread between 100 and 3000 ms with high sampling rate around the TIs that nulls the GM layers. The resolution was 510 μm cubic with TR of 4000 ms and TE of 4.9 ms (effective TE of 41.4 ms). The acquisition time for each TI was 80 s and the total acquisition time was 80 min. Image analysis The IR datasets (either human or rats) were fitted to the conventional inversion recovery decay function (Barral et al., 2010) (using an in-house code written in Matlab, (Mathworks, Natick, MA)) with possible multiple T1 components, on a voxel-by-voxel basis: S. Lifshits et al. NeuroImage 164 (2018) 112–120 113 ∑TI M eM( ) = 0 ∙(1 − 2 )i j j TI T− /i j1 (1) where TIM( )i is the magnetization at the i-th inversion recovery image, M0j is the predicted magnetization at TI=0 ms for each T1 component (j) in the voxel, and T j1 is the longitudinal relaxation time for each T1 component. j was set to 1 for the high resolution experiments and to 7 for the low resolution experiments. In the low-resolution experiment, j was set to 7 following the assumptions that there are 7 T1 components in the tissue – 1 for CSF, 1 for WM and heavily myelinated layer of the cortex and additional 5 cortical layers. Calculation of the T1component probability maps A histogram of the whole brain T1 values was used to compute T1component probability maps. A Gaussian mixture model was used to fit the whole brain T1 values (here we used the gmdistribution fit function in Matlab). Using a Bayesian index criterion (BIC) we found that the best fit the T1-histogram provides 7 or 8 distinct T1 components (presumably 1 for CSF, 1 for WM and 5–6 for GM). Fig. 1 shows such analysis for the entire dataset: high resolution (A and C) vs lowresolution (B and D) and human (A and B) compared with rat (C and D). In general we found 7 Gaussian classes: 1 for white matter and layer 6 (with T1 values lower than ~600 ms), 5 for additional gray matter layers (noted L 1–5) and one for CSF (not shown on graph). Based on the Gaussian mixture analysis we could compute the probability of each voxel to be assigned as each of the classes. Eventually it was possible to create class specific probability maps. Therefore, each image voxel, regardless of the acquisition resolution, is represented by a vector of the abovementioned T1 class's probabilities. To visualize this information we projected the T1 class probabilities onto a surface representation of the brain. Another way we used to visualize the data is by clustering. K-means clusters employed to segment the cortex to the different layers based on the class probability and location along the cortex (k was set to 7 or 8 based on the number Table 1 Experiment Low-res rat High-res rat Low-res human High-res human Magnet 7T Biospec Bruker 3T Siemens Prisma TR 4,000 ms 4,000 ms 11,000 ms 10,000 ms TE (eff. TE) 4.9 (41.4) ms 5.92 (29) ms 30 ms 14ms TI rangea 50–3000 ms 100–3000 ms 50–3000 ms 230-1380 ms No. of TI measurements 60 18 107 7 Resolution 0.51×0.51×0.51 mm3 0.11×0.11×0.11 mm3 3x3×3 mm3 0.43×0.43×1.5 mm3 a The full list of TI's is given in the Supplementary information. Fig. 1. T1 histograms of the different experiments performed in this study. (A) Histogram of T1 values for the entire brain of the high resolution human experiments. The histogram demonstrate a mixture of Gaussians with the gray lines indicates the histograms for each subject and the black line – the average of all subjects. The arrows are used for visual guidance on the centers of the different Gaussian classes (1–6 for gray matter and their layer classification). It is interesting to note that there is increased variability across subject especially at the low-T1 range (lower layers of the cortex and white matter). Similar analysis on the T1 data of the low resolution experiments (B) shows that his variability reduces (probably since we model the partial volume effect). (C) and (D) shows the same but for the rat brain high resolution (C) and low-resolution (D). Note that the centers of the Gaussians remained the same both in the high resolution (A and C) compared with the low-resolution (B and D) indicating that the low-resolution T1 mapping is adequate for separation between the layers. S. Lifshits et al. NeuroImage 164 (2018) 112–120 114 of Gaussians in the T1 histogram). Results In this paper we wish to explore PVE effects in cortical layers imaging through T1 relaxation mapping. At first step we acquired a high resolution (110 μm3 ) on excised rat brain to investigate the T1 layers with minimal PVE. Then we demonstrate composition based analysis on two in-vivo human data sets (one with meso-scale resolution and one with low resolution). High resolution rat data High cubic resolution MRI (110 μm3 ) of a fixated rat brain was acquired to demonstrate the utility of inversion recovery and T1 MRI to visualize and quantify the layers. The purpose of this experiment was to demonstrate that: a) Different layers have different and distinct T1 values; b) The layer distribution has neuro-anatomical meaning. Fig. 2 shows analysis of the high resolution IR data set. Fig. 2A depicts the T1 map of one representative slice with an enlargement of the S1BF cortex demonstrating the variation in T1 perpendicular to the cortex (numerically demonstrated in Fig. 2B). It is obvious that T1 changes across the cortex are not just a ‘gradual’ change from CSF's high T1 to white matter's low T1. The T1 perpendicular to the cortical stripe changes in a stepwise manner with a trend of reduction in T1 towards the white matter probably due to increased myelination in deeper layers (Fig. 2B). The T1 class probability maps were used to segment the cortex to the different layers based on the values and location along the cortex (using k-means clustering). Following clustering, the thickness of each layer perpendicular to the cortical surface could be calculated. Fig. 2C shows the projection of each layer thickness onto a surface representation of cortex. The layer thickness surfaces show different regional distribution across the cortex (and hence demonstrate the sensitivity of the T1 layers). By comparing the layer probability maps with the cytoarchitectonic atlas of brain regions (Fig. 2D), it is possible to observe that different cyto-architectonic areas also differ in their T1 layer thickness. Noteworthy are the missing layer 4 from the motor cortex (arrow 1), the thick layer 4 in visual cortex and sensory cortex (arrows 2 and 3) and thick layer 2 in the cingulate cortex (arrow 4). Meso-scale human data Cubic resolution of 110 μm in the human brain is in-achievable, yet the cortical width in humans is about twice that of the rat. Therefore comparable resolution to Fig. 2 in the human brain would be in the order of 200–300 μm – which cannot be achieved in-vivo for this protocol even with high-field magnets (7T and above). Alternatively, we have scanned human subjects with the best achievable resolution in an acceptable scanning time (~1 h) of 430×430 μm in-plane resolution and slice thickness of 1.5 mm. Analysis of T1-layers in such setup were reported previously (Barazany and Assaf, 2011). Yet, it is well accepted that at such resolution significant PVE between the layers and adjacent tissues is apparent and a more sophisticated model should be used to infer sub-voxel layers distribution. In the following analysis, instead of computing single T1 value for each voxel and then assigning each T1/voxel to a layer, we performed a multi T1 analysis revealing the sub-voxel T1 composition. Conventionally such multi-component analysis can be achieved by acquiring multiple inversion recovery images (with large number of inversion times) to allow multi T1 analysis (as in Eq. (1) when j is set to 2 or higher). Yet, as only 7 TIs were acquired in the meso-scale resolution human data, we could not fit more than T1 value per voxel. To overcome this obstacle we used the T1 classes as shown in Fig. 1. Here, we used the T1 classes as a-priori knowledge and simulated, per class, the IR signal dependency (sampled several time over the Gaussian distribution of each class). That signal dependency was used as a training set. Once this classifier was trained, we could compute the Fig. 2. High resolution T1 analysis for the rat brain. (A) T1 map of one representative slice along with enlargement of the S1BF area. (B) A profile of T1 values perpendicular to the s1BF area showing the borders between the layers. (C) Projection of the layer class probability maps (L1-L6) onto a cortical surface. Arrows indicate areas of interest: 1 – M1 cortex, 2 – V1 cortex, 3 – S1 cortex, 4 – Cingulate cortex. (D) shows the orientation surface projection for a cyto-architectonic atlas of the rat brain indicating the entorhinal and perirhinal cortices (dark blue), visual cortex (blue), cingulate cortex (cyan), temporal cortex (mainly auditory cortex, green), sensory cortex (all parts, orange), motor cortex (all parts, red), piriform cortex and adjacent regions (olive green), and insular cortex (yellow). S. Lifshits et al. NeuroImage 164 (2018) 112–120 115 class's probability values for each voxel. Such analysis is shown in Fig. 3 where 5 T1 classes were identified. In this analysis, the T1 component in Fig. 1 with the longest T1 values has the highest probability at the superficial parts of the cortex and the probability of the components with the shorter T1 values increases at voxels that are closer to the white matter surface. The T1 classes' probability maps represent the sub-voxel composition of each class (i.e., layer). The sub-voxel composition can be used to enhance the image resolution (similar to super-resolution) by solving a regularized optimization problem. In this analysis each voxel was divided into 4-sub voxels and each one of them was assigned to a class based on the composition analysis of the parent voxel. The class of each sub-voxel was assigned by keeping the mean value of all sub-voxels with minimum deviation from the parent voxel. This process provided new classified layer images at interpolated in-plane resolution of 215 μm. In this analysis we used a total variation regularization function, which preserves the edges in the image and ensure continuity of the classes in the spatial space. The optimization problem was solved by using the FISTA algorithm (Beck and Teboulle, 2009) providing the enhanced resolution images (Fig. 4). Fig. 4A shows a segmentation of the cortex into 5 classes based on majority vote of computed T1 values (i.e. the class that received the highest probability). The white circles in Fig. 4A represent areas where the continuity of some layer classes was disturbed as a consequence of PVE. In these voxels, a class was assigned by majority vote but the probability for other classes was not insignificant. The enhanced resolution image (Fig. 4B) shows that even with subtle division into 4 sub-voxel – continuity of the layers can be obtained providing better Fig. 3. Layer class probability maps. Following training a classifier based on the Gaussian mixture model applied on the data shown in Fig. 1, the probability for belonging to each class was computed. Class 1 shows high probability and superficial parts of the cortex while the probability of the other classes increases approaching the inner surface of the cortex. Fig. 4. Enhanced cortical layer resolution. (A) Majority votes of the 5 Gy matter classes in the original resolution. Each color represent a class who had the highest probability in each voxel. (B) Enhanced resolution majority vote image where each voxel in (A) was sub-divided into 4 voxels and a class was assigned based on the probabilities (composition) of the classes. Circles represents areas where low probability classes could have been visualized only in the enhanced resolution image (see class 4 in cyan). S. Lifshits et al. NeuroImage 164 (2018) 112–120 116 mean to visualize the layers. Fig. 5 shows the same analysis along with FreeSurfer segmentation into neuro-anatomical regions zoomed at the bottom panel. Here we focused on an area where transition between two neuro-anatomically distinct regions occurs. The transition, that is also reflected by a change in the layer thickness distribution is better observed in the enhanced image (Fig. 5B) compared with the original resolution image (Fig. 5A). Low resolution human data Figs. 1–5 show that T1 can be used to characterize layers in the cortex yet PVE remains a major obstacle before establishing this method for whole brain layer visualization. In addition, the issue of scanning time provides another drawback for this method. While the high resolution rat experiment nicely demonstrated the ability of IRMRI and T1 mapping to visualize layer and identify the neuroanatomical brain regions, similar high resolution data, in the human brain (ex-vivo or in-vivo) will be extremely difficult to achieve, scanning time wise. Even the combination of meso-scale resolution and composition analysis (Figs. 4 and 5) does not provide optimal framework for layer analysis. Therefore, in the following analysis we focused on extracting more accurately the T1 layers from composition analysis. To have enough data points to allow robust and accurate composition analysis, image resolution had to be reduced to remain within acceptable scanning time-scale. Here we have sampled 107 inversion times per voxel, while each voxel was scanned at 3 mm3 resolution obviously containing multiple layers (and even all of them in some cases) within a single voxel. Multi T1 component function was fitted to each voxel with j=7 (in Eq. (1)). We choose 7 as in a preliminary test we didn’t find voxels with more than 7 distinct T1 components. Eventually the T1 maps consisted of up to 7 T1s per voxels - ranging from one T1 class in homogeneous voxel such as in the CSF up to 7 classes in voxels that include the entire cortical stripe. In each voxel, aside from the T1 values, we also computed the relative fraction of each class (M0j in Eq. (1)). Using this data we computed a whole brain T1 histogram (similar to the one shown in Fig. 1) where at least 8 Gaussian classes were revealed − 1 for white matter, 1 for CSF (very broad distribution) and 6 different Gaussians in the gray matter T1 regime. As in the previous analyses given in Figs. 2–5, we used the Gaussian mixture analysis to assign the T1 values to the different classes. Once assigned, we could compute class specific maps where the values in each voxel represent the volume fraction of this class from the multi T1 analysis. These class weighted maps were then projected onto a surface representation of the cortex (Fig. 6). The layer surface representation reveals interesting features on the T1 layers – For example the high intensity of layer 4 around the occipital cortex (arrow 1) and primary visual areas (Fig. 6A), the missing layer 4 from the primary motor area (arrow 2), the high intensity of layer 6 in temporal regions (arrow 3), the small fraction of layers 1 and 2 compared to the others and the almost even distribution of layer 3 across the cortex. Quantification of the observations shown in Fig. 6 was performed by segmenting the cortex into 6 types: Allocortex, agranular cortex, and additional 4 types with increasing amount of granularity of the cortex as suggested previously (Barbas and Rempel-Clower, 1997; Brodmann and Garey, 1999; Triarhou, 2007; von Economo and Koskinas, 1925). Fig. 7 depicts the layer composition in each of the 6 cortical types. In this analysis we extracted, per cortical type (see Figure caption for the list of regions), the relative volume fraction of each layer out of the total gray matter volume area in that region (i.e. relative layer width). The total volume fraction of a layer was computed by summing the probability of a certain layer for each region (from the layer probability maps as shown in Fig. 6) and multiplying it by the voxel volume. This number was then divided by the total volume of gray matter for each region to produce the relative weights of the 6 layers in each cortical region. This kind of analysis avoids complicated computation of layer Fig. 5. Comparison of layer class distribution with neuro-anatomical atlas. (A) Shows the FreeSurfer atlas segmentation into neuro-anatomical regions with a section of the frontal cortex enlarged below. (B) Majority vote of the different classes in the original and (C) enhanced resolution. The transition between the two areas (blue and red in (A)) is marked in the white arrows where class 4 almost disappear in the red region and class 1 is dramatically narrowed in the blue region. S. Lifshits et al. NeuroImage 164 (2018) 112–120 117 variation along the folded cortex and also minimize the bias of variation in cortical thickness across different regions. This figure nicely demonstrate the ability to distinguish between the types based of the differences in layer 2 (decreasing in relative width with increase in cortical granularity) and layer 4 (increase in relative width with the increase in granularity). In order to investigate the information embedded in the layer composition analysis, we have compared the layer class probabilities in different Brodmann areas (BAs). This was achieved by registering both to the same volume space and averaging each layer class probability over each BA (see bottom panel in Fig. 6). Analysis of variance (ANOVA) of the 6 layer class probabilities x 37 Brodmann areas (BA) revealed that most of them are distinguishable from the others (red square in Fig. 8). More importantly, out of the 76 regions pairs (blue outlined squares) that are adjacent, we could separate between 42 pairs (p < 0.05) (red squares with blue outline) with additional 11 showing a trend of difference (orange square with blue outline). Analysis of the variance across subjects revealed that the highest layer variability exists in the following regions: prefrontal cortex, orbitofrontal cortex, and cingulate cortex. The most homogenous regions (layer probability wise) were: angular gyrus, infer frontal gyrus, and associative visual cortex. Discussion The cortical layers can be resolved by T1 MRI – a phenomenon that was already demonstrated several times in numerous papers, most of them on a limited region of the brain (Barazany and Assaf, 2011; Barbier et al., 2002; Clare and Bridge, 2005; Geyer et al., 2011; Turner et al., 2008). Yet the spatial resolution is the main limitation for cortical layer imaging since accurate visualization and quantification of the layers requires acquisition resolution that is far beyond the limits of current MR technology. Indeed, high-resolution MRI does not allow adequate characterization of the layers, yet it does allow the extraction of laminar features (e.g. myelination content) as demonstrated in several works (Dinse et al., 2015; Fracasso et al., 2016; Lutti et al., 2014; Sereno et al., 2013; Waehnert et al., 2016, 2014) underscoring the utility of cortical sub-component characterization. In this work we wish to push the limits of MRI even further and provide better and potentially more accurate quantification of the layers themselves but with a different approach. Reverting to the conventional definition of resolution might help to understand the concept behind the approach used in this work: how close objects can be to one another and still be resolved or the capability of making parts of an object distinguishable Fig. 6. Layer specific surface projection maps of the human brain. (A) Occipital view of the layer class specific probability maps showing high probability values in the expected regions for layer 4 (especially in V1, see Brodmann atlas surface projection at the bottom). (B) Right side view of the layer class specific probability maps. Arrows indicate areas of interest: high probability for layer 4 in the visual cortex (arrow 1), low probability for layer 4 in the motor cortex (arrow 2), and high probability for layer 6 in the temporal cortex (layer 3). Fig. 7. Quantification of layer composition at different types of cortices. In this figure we quantified the normalized layer fraction out of the total gray matter volume (relative layer width) of the Piriform cortex (as a representative of the allo-cortex, BA27), Pre-motor cortex (BA6 representing the agranular cortex) and in addition BA20 (Inferior temporal gyrus, S. Granular 1), BA11 (Orbitofrontal cortex, S. Granular 2), BA19 (Associative visual corte, S. Granular 3) and BA17 (Primary visual cortex, Granular cortex) representing increasing degrees of cortical granularity). The graph demonstrate the ability the proposed analysis and specifically the quantification of layers 2 and 4 to distinguish between the different cortical types. The brain surface projection on the top right side represent different brain regions that correspond to each of the cortex types (see color index on the left). This map is based on similar segmentation shown in Beul and Hilgetag (2014) based on Von Economo (2009). S. Lifshits et al. NeuroImage 164 (2018) 112–120 118 (and quantifiable). In this paper we argue that distinguishing and quantifying the cortical layers (via their T1 differences) is feasible in the relaxation time domain rather than in the spatial image domain. While the best demonstration of high spatial resolution in-vivo was resolving some properties of the layers (e.g. the stripe of Gennari), accurate T1 measurements are able to distinguish between all layers regardless the spatial resolution, once measured in the time domain and observed in the spatial domain. With the technological advanced in MRI (i.e. high magnetic field scanners, efficient r.f. excitation and receiving schemes, strong gradient capacity) signal to noise ratio (SNR) increased dramatically. Trivial consequence of SNR improvement is increasing image resolution leading to visualization of finer details in the image. Indeed better resolution is a key feature in imaging as summarized in the well-known phrase: ‘seeing’ is ‘believing’. Yet, let's consider one of the most popular MRI methods: diffusion imaging – which measure a spatial parameter (translational displacement) whose scale is orders of magnitude smaller than high-resolution imaging (μm vs. mm). This is achieved by measuring the signal in q-space (reciprocal spatial domain) rather than in the k-space (the spatial frequency domain). The micron scale resolution of diffusion imaging is achieved not by making an image in the micron scale but rather characterizing a phenomenon that happen on the micron-scale. We suggest the same solution for cortical layer imaging - measuring the signal in poor spatial resolution but with high resolution in the T1 relaxation domain and achieving good separation between the layers there. Continuing the analogy to diffusion imaging – one of the main difficulties in diffusion imaging is providing biological interpretation to the diffusion measures as the relation between brain anatomy/physiology and biophysical properties of water motion is indirect. The same happen to some extent with T1 measures of the layers. The cortical layers are defined, traditionally, as a measure of cyto-architecture and/ or myelo-architecture. The relation between T1 and these histological properties is, as expected, complicated. Yet, there are many evidences in the literature that T1 is highly influenced by the degree of myelination (Barazany and Assaf, 2011; Barbier et al., 2002; Turner et al., 2008) therefore the observed T1 layers might be highly correlated with myelo-architecture. Indeed previous studies (Dinse et al., 2015; Fracasso et al., 2016; Lutti et al., 2014; Sereno et al., 2013), as well as the comparison with the cyto-architectonic atlases shown in this paper (Figs. 2, 5–7) indicate that the T1 layers resemble the properties of the cyto and myelo-architecture. An excellent demonstration of that was provided very recently by Waehnert et al. (2016) indicating the T1 mapping provide invaluable, subject specific, information of cortical architecture. However, one cannot claim that T1 is a measure of any of these histological properties as it is not – it measure the relaxation rate of excited spins to the surrounding. While the ‘surrounding’ indeed represents the cyto and myelo architecture, it is also represents many other factors. Therefore, we suggest to refer to cortical layer imaging by T1 as T1 layers bearing in mind that they resemble the traditional cortical layers but not a direct measure of them. The vague biological interpretation is one limitation of using lowresolution multi T1 analysis to quantify the layers. Yet, the biggest limitation of this approach is in its concept: the need to ‘see’ the layers. The use of low resolution which does not allow direct imaging of the layers may raise doubts about the validity of this observations. This is a conceptual limitation, yet there are some technical issues: 1) the need for dozens of inversion time data points to allow robust estimation of the T1 values. 2) Errors in the estimation might lead to misclassification of a T1 component and “noise” in the measured layers. Finally, we argue that the layers can be better separated in the T1 relaxation domain rather than the spatial domain. This is true, yet even in the relaxation domain the separation is not optimal and there is some overlap between the T1 components. Shifting to higher magnetic fields (7T and higher) might reduce this overlap as higher field spreads the T1 values. Fig. 8. ANOVA of layer distribution x Brodmann areas. Analysis of variance between the averaged probability for each layer class x 37 Brodmann areas revealed significant differences (p < 0.05) in the red squares and trend towards significance in the orange squares (0.1 < p < 0.05). Blue squares represent region that adjacent to one another. S. Lifshits et al. NeuroImage 164 (2018) 112–120 119 The need to ‘see’ the layers is an obstacle that neuroimagers need to overcome in order to be able to measure the layers. Even when ultrahigh magnetic fields will be constructed, imaging the human brain at 50 μm resolution (which might allow direct visualization of the layers) will lead to additional computational and statistical challenges, as each brain will be represented by an enormous amount of data. Alternatively, combined with surface projection (as shown in Figs. 2 and 6), low-resolution T1 bears great potential for large population studies with moderate computational needs and low statistical risks (in terms of multiple comparisons). The high sensitivity of the T1 layers to the known neuroanatomy and even some kind of across species homology in the layer distribution (compare Figs. 2 and 6) strengthen the need to implement this methodology to further study its sensitivity to various conditions (and its relation to cognition and pathological conditions). Acknowledgements The authors wish to thank the NOFAR program of the Israel innovation authority, the Israeli ministry of Economy and Industry (50647) who funded part of the project (High resolution composition analysis). Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at doi:10.1016/j.neuroimage.2017.02.086. References Barazany, D., Assaf, Y., 2011. Visualization of cortical lamination patterns with magnetic resonance imaging. Cereb. Cortex. Barbas, H., Rempel-Clower, N., 1997. Cortical structure predicts the pattern of corticocortical connections. Cereb. Cortex 7, 635–646. Barbier, E.L., Marrett, S., Danek, A., Vortmeyer, A., van Gelderen, P., Duyn, J., Bandettini, P., Grafman, J., Koretsky, A.P., 2002. Imaging cortical anatomy by highresolution MR at 3.0T: detection of the stripe of Gennari in visual area 17. Magn. Reson. Med. 48, 735–738. Barral, J.K., Gudmundson, E., Stikov, N., Etezadi-Amoli, M., Stoica, P., Nishimura, D.G., 2010. A robust methodology for in vivo T1 mapping. Magn. Reson. Med. 64, 1057–1067. Beck, A., Teboulle, M., 2009. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems. IEEE Trans. Image Process. 18, 2419–2434. Beul, S.F., Hilgetag, C.C., 2014. Towards a "canonical" agranular cortical microcircuit. Front. Neuroanat. 8, 165. Brodmann, K., Garey, L., 1999. Brodmann's Localisation in the cerebral cortex. Imperial College Press, London, River Edge, NJ. Clare, S., Bridge, H., 2005. Methodological issues relating to in vivo cortical myelography using MRI. Hum. Brain Mapp. 26, 240–250. Clark, V.P., Courchesne, E., Grafe, M., 1992. In vivo myeloarchitectonic analysis of human striate and extrastriate cortex using magnetic resonance imaging. Cereb. Cortex 2, 417–424. Dinse, J., Hartwich, N., Waehnert, M.D., Tardif, C.L., Schafer, A., Geyer, S., Preim, B., Turner, R., Bazin, P.L., 2015. A cytoarchitecture-driven myelin model reveals areaspecific signatures in human primary and secondary areas using ultra-high resolution in-vivo brain MRI. Neuroimage 114, 71–87. Duyn, J.H., van Gelderen, P., Li, T.Q., de Zwart, J.A., Koretsky, A.P., Fukunaga, M., 2007. High-field MRI of brain cortical substructure based on signal phase. Proc. Natl. Acad. Sci. USA 104, 11796–11801. Eickhoff, S., Walters, N.B., Schleicher, A., Kril, J., Egan, G.F., Zilles, K., Watson, J.D., Amunts, K., 2005. High-resolution MRI reflects myeloarchitecture and cytoarchitecture of human cerebral cortex. Hum. Brain Mapp. 24, 206–215. Fracasso, A., van Veluw, S.J., Visser, F., Luijten, P.R., Spliet, W., Zwanenburg, J.J., Dumoulin, S.O., Petridou, N., 2016. Lines of Baillarger in vivo and ex vivo: myelin contrast across lamina at 7T MRI and histology. Neuroimage 133, 163–175. Geyer, S., Weiss, M., Reimann, K., Lohmann, G., Turner, R., 2011. Microstructural parcellation of the human cerebral cortex - from Brodmann's post-mortem map to in vivo mapping with high-field magnetic resonance imaging. Front Hum. Neurosci. 5, 19. Goncalves, N.R., Ban, H., Sanchez-Panchuelo, R.M., Francis, S.T., Schluppeck, D., Welchman, A.E., 2015. 7 T FMRI reveals systematic functional organization for binocular disparity in dorsal visual cortex. J Neurosci. 35, 3056–3072. Hof, P.R., Glannakopoulos, P., Bouras, C., 1996. The neuropathological changes associated with normal brain aging. Histol. Histopathol. 11, 1075–1088. Kalus, P., Senitz, D., Beckmann, H., 1997. Cortical layer I changes in schizophrenia: a marker for impaired brain development? J Neural Transm. 104, 549–559. Kandel, E.R., Schwartz, J.H., Jessell, T.M., 2000. Principles of Neural Science 4th ed.. McGraw-Hill, Health Professions Division, New York. Koopmans, P.J., Barth, M., Norris, D.G., 2010. Layer-specific BOLD activation in human V1. Hum. Brain Mapp. 31, 1297–1304. Kordower, J.H., Chu, Y., Stebbins, G.T., DeKosky, S.T., Cochran, E.J., Bennett, D., Mufson, E.J., 2001. Loss and atrophy of layer II entorhinal cortex neurons in elderly people with mild cognitive impairment. Ann. Neurol. 49, 202–213. Leuze, C.W., Anwander, A., Bazin, P.L., Dhital, B., Stuber, C., Reimann, K., Geyer, S., Turner, R., 2012. Layer-specific intracortical connectivity revealed with diffusion MRI. Cereb. Cortex.. Lutti, A., Dick, F., Sereno, M.I., Weiskopf, N., 2014. Using high-resolution quantitative mapping of R1 as an index of cortical myelination. Neuroimage 93 (Pt 2), 176–188. Masliah, E., Mallory, M., Hansen, L., Alford, M., Albright, T., DeTeresa, R., Terry, R., Baudier, J., Saitoh, T., 1991. Patterns of aberrant sprouting in Alzheimer's disease. Neuron 6, 729–739. Olman, C.A., Harel, N., Feinberg, D.A., He, S., Zhang, P., Ugurbil, K., Yacoub, E., 2012. Layer-specific fMRI reflects different neuronal computations at different depths in human V1. PLoS One 7, e32536. Sereno, M.I., Lutti, A., Weiskopf, N., Dick, F., 2013. Mapping the human cortical surface by combining quantitative T(1) with retinotopy. Cereb. Cortex 23, 2261–2268. Triarhou, L.C., 2007. A proposed number system for the 107 cortical areas of Economo and Koskinas, and Brodmann area correlations. Stereotact. Funct. Neurosurg. 85, 204–215. Turner, R., Oros-Peusquens, A.M., Romanzetti, S., Zilles, K., Shah, N.J., 2008. Optimised in vivo visualisation of cortical structures in the human brain at 3 T using IR-TSE. Magn. Reson. Imaging 26, 935–942. Vogt, O., 1911. Die myeloarchitektonik des isocortex parietalis. J. Psychol. Neurol. 18, 379–390. Von Economo, C., 2009. Cellular Structure of the Human Cerebral Cortex. KARGER, Berlin. von Economo, C.B., Koskinas, G.N., 1925. The Cytoarchitectonics Of The Adult Human Cortex. Julius Springer Verlag, Vienna and Berlin. Waehnert, M.D., Dinse, J., Schafer, A., Geyer, S., Bazin, P.L., Turner, R., Tardif, C.L., 2016. A subject-specific framework for in vivo myeloarchitectonic analysis using high resolution quantitative MRI. Neuroimage 125, 94–107. Waehnert, M.D., Dinse, J., Weiss, M., Streicher, M.N., Waehnert, P., Geyer, S., Turner, R., Bazin, P.L., 2014. Anatomically motivated modeling of cortical laminae. Neuroimage 93 (Pt 2), 210–220. Walters, N.B., Egan, G.F., Kril, J.J., Kean, M., Waley, P., Jenkinson, M., Watson, J.D., 2003. In vivo identification of human cortical areas using high-resolution MRI: an approach to cerebral structure-function correlation. Proc. Natl. Acad. Sci. USA 100, 2981–2986. S. Lifshits et al. NeuroImage 164 (2018) 112–120 120