Parcellation of the Human Cerebral Cortex Based on Molecular Targets in the Serotonin System Quantified by Positron Emission Tomography In vivo

Abstract Parcellation of distinct areas in the cerebral cortex has a long history in neuroscience and is of great value for the study of brain function, specialization, and alterations in neuropsychiatric disorders. Analysis of cytoarchitectonical features has revealed their close association with molecular profiles based on protein density. This provides a rationale for the use of in vivo molecular imaging data for parcellation of the cortex with the advantage of whole-brain coverage. In the current work, parcellation was based on expression of key players of the serotonin neurotransmitter system. Positron emission tomography was carried out for the quantification of serotonin 1A (5-HT1A, n = 30) and 5-HT2A receptors (n = 22), the serotonin-degrading enzyme monoamine oxidase A (MAO-A, n = 32) and the serotonin transporter (5-HTT, n = 24) in healthy participants. Cortical protein distribution maps were obtained using surface-based quantification. Based on k-means clustering, silhouette criterion and bootstrapping, five distinct clusters were identified as the optimal solution. The defined clusters proved of high explanatory value for the effects of psychotropic drugs acting on the serotonin system, such as antidepressants and psychedelics. Therefore, the proposed method constitutes a sensible approach towards integration of multimodal imaging data for research and development in neuropharmacology and psychiatry.


Introduction
The classification of the cerebral cortex into discrete regions bears several advantages for the study of brain function. Ideally, the distinction of parcels is reproducible across individuals and laboratories, and accompanied by a nomenclature, such that specific scrutiny and scientific discourse is facilitated. Historically, Brodmann areas constitute the most prominent attempt in this direction and are based on Brodmann's observations of cytoarchitectonical features of the cerebral cortex (Zilles and Amunts 2010). The advance of new technologies that allow for the measurement of different properties of human brain tissue has spurred attempts to improve parcellation. Especially in the field of brain imaging, a reliable reference system is indispensable (Amunts et al. 2014), as exemplified by the effort of Talairach-Tournoux to translate Brodmann areas into a standardized coordinate system (Talairach and Tournoux 1988). Based on macroanatomical landmarks, this atlas was transformed into the widely used MNI template created from an average of a multitude of MRI scans, which resulted in a more accurate representation of the population's neuroanatomy (Chau and McIntosh 2005). Atlases based on averages of multiple subjects do not entirely accommodate inter-individual variation, especially if their alignment is based on macro-anatomical landmarks, such as gyrification of the cortex (Desikan et al. 2006;Destrieux et al. 2010), as these do not necessarily coincide with cytoarchitectonical borders (Uylings et al. 2005). One promising solution has been proposed and validated using data from the human connectome project, which entails the acquisition of standardized structural and functional MR sequences for each subject and enables an individualized delineation of cortical areas using a machine-learning classifier (Glasser et al. 2016). Another invaluable approach is the construction of probabilistic cytoarchitectonical atlases derived from computational analysis of histological sections in a set of several brains (Eickhoff et al. 2005). The finding that interindividual variations in functional MRI activation could be reduced by representing the cerebral cortex in a two-dimensional coordinate system has initiated the trend of surface-based analysis Fischl, Sereno, Tootell, et al. 1999). The benefits of this methodology have led to its extensive application in processing of structural and functional MRI, as well as molecular imaging data (Greve et al. 2014) and the investigation of structural and functional connectivity.
In the domain of connectivity analysis, the identification of parcels using resting-state fMRI signal and graph theoretical (Power et al. 2011), independent component (Yeo et al. 2014, Markov random field (Schaefer et al. 2017), or spectral clustering analysis (Craddock et al. 2012) has yielded appealing results with regard to reproducibility and elucidation of information processing in the cerebral cortex. On the other hand, hierarchical clustering approaches to parcellation may parallel the development of specialization in the cerebral cortex, as demonstrated in structural connectivity data (Moreno-Dominguez et al. 2014). A MRI study in twins which delineated clusters based solely on the heritability of cortical surface area similarly revealed a hierarchical organization of cortical areas (Chen et al. 2013). Other invaluable approaches are derived from the analysis of functional connectivity networks in comparison to gene expression patterns (Ganglberger et al. 2017).
It could be argued that the field of parcellation of the cortex has undergone a gradual paradigm shift from dissecting maximally distinct areas towards the identification of the topology of coexisting sub-systems. Interestingly, borders between regions defined by cytoarchitectonical features were found to coincide to a great extent with those defined by receptor profiles using autoradiography (Zilles et al. 2002), suggesting that functional differentiation of the cerebral cortex may be paralleled by receptor compositions. Not only does this constitute complementary information and render multimodal support for the parcellation derived from cytoarchitectonics, but it provides a rationale for the utilization of molecular imaging data in vivo to study brain ontology. A linear relationship between protein densities measured using autoradiography and outcome measures derived from molecular imaging methods such as positron emission tomography (PET) has been demonstrated for several proteins Beliveau et al. 2017). A critical advantage of molecular imaging is that acquisition of comprehensive whole-brain data is possible within a single scan, as opposed to effortful sectioning and staining of selected regions. We therefore propose clustering using PET imaging data of molecular targets as a complementary approach to parcellation of the cortical surface. We expect that clusters defined on the expression of molecular targets will be of special relevance for the analysis of pharmacological imaging data and therefore facilitate research and development in this field. The presented analysis is constrained to PET data on four proteins involved in serotonergic neurotransmission for which a representative sample is available. Inarguably the serotonergic system exerts a strong influence during brain development, modulation of a large set of brain functions and behavior, and is one of the most important targets of therapeutic intervention in neuropsychopharmacology. Extensive application of the proposed parcellation methods to different neurotransmitter systems can be expected with the emergence of data sharing initiatives in the molecular imaging field (Knudsen et al. 2016). The resulting clusters were compared to established cytoarchitectonic and functional brain parcellations and interpreted in light of pharmacological effects targeting the investigated binding proteins.

Participants and Imaging Procedures
The dataset used in this analysis comprises PET and structural 3T MR measurements of 108 healthy participants in total. Each subject underwent one PET scan for one target. For details on demographics and data acquisition see Table 1, or previous publications covering the data used here (Elmenhorst et al. 2012;Lanzenberger et al. 2007Lanzenberger et al. , 2012Gryglewski et al. 2017;Vanicek et al. 2017). Subjects were recruited by advertisement and provided written informed consent according to the study procedures approved by the Ethics Committees at the Medical University of Vienna, Austria, and the Medical Faculty of the University of Düsseldorf, Germany. Participants underwent medical examinations to assure their physical and mental health. Subjects had no history of psychopharmacological drug treatment, substance abuse or diagnosis of psychiatric conditions according to the Diagnostic and Statistical Manual of Mental Disorders IV (DSM-IV). Scans were performed regardless of menstrual cycle of female participants.
PET imaging was carried out using optimized protocols with application of highly specific tracers for quantification of the 5-HT 1A receptor ([carbonyl-11 C]WAY-100635), 5-HT 2A receptor ([ 18 F] altanserin), MAO-A ([ 11 C]harmine) and serotonin transporter 5-HTT ([ 11 C]DASB). [carbonyl-11 C]WAY-100635 has a tenfold lower affinity for dopamine D4 receptors than for 5-HT 1A receptors (Chemel et al. 2006;Martel et al. 2007). Based on the low expression of D4 receptor mRNA in the cortex, we do not expect a relevant influence on binding potentials in this region (Matsumoto et al. 1996). Changes in [ 11 C]harmine binding up to 17% were observed in one study after altering endogenous monoamine concentrations pharmacologically (Sacher et al. 2012). The other tracers used were shown to be insensitive to changes in endogenous ligand concentrations, with changes in [ 11 C]DASB binding being observed only after extreme pharmacological challenges performed in animal studies (Paterson et al. 2010).

Image Preprocessing
PET images were motion corrected in SPM12 (Wellcome Trust Centre for Neuroimaging, London, UK; http://www.fil.ion.ucl.ac. uk/spm) by rigid realignment of every frame image to a median image, which consisted of a motion-free period, selected by visual inspection. Duration of frames were between 5 and 300 s, depending on the protocol and time point (frame length was increasing towards the end of the measurements to improve signal to noise ratio). The median PET images were coregistered to the individual MR images. The cortical surface was reconstructed using FreeSurfer 6.0 (Harvard Medical School, Boston, USA; http://www.surfer.nmr.mgh.harvard.edu) with individual T1-weighted MR images serving as input. In addition, every reconstruction was visually inspected and manually corrected for possible inaccuracies in segmentation of the cortical surface by setting control points or erasing falsely segmented areas and re-initiating reconstruction according to published procedures (McCarthy et al. 2015). Volume to surface projection was performed using the FreeSurfer function mri_vol2surf. The obtained registration parameters were subsequently applied to the motion-corrected dynamic PET images resulting in dynamic PET data in fsaverage surface space.

Quantification of Molecular Targets
The dynamic PET data registered to the cortical surface served as input for the quantification of protein distributions, which was performed in MATLAB 8.2 (https://www.mathworks.com) using models specified below. PET data was smoothed prior to quantification using an 8-mm kernel in surface space. The resulting outcome measures (binding potential, BP ND , BP P , V T ) are proportional to absolute density of available receptors (Innis et al. 2007).

5-HT 1A Receptor
The cortical 5-HT 1A receptor binding was quantified by applying the multilinear reference tissue model (MRTM2) (Ichise et al. 2003). The insular cortex was used as high-uptake and the cerebellar white matter as reference region, due to the low 5-HT 1A concentration in this area (Parsey et al. 2005). The choice of quantification method for 5-HT 1A receptor binding potentials may affect results of group comparisons based on variability of receptors in the reference region (Parsey et al. 2010). This issue is not relevant to the current analysis, as a potentially different scaling of data is removed by normalization prior to clustering.

5-HT 2A
Receptor 5-HT 2A receptor binding potential (BP P ) was quantified using radioactivity in brain and blood plasma at equilibrium during tracer infusion (Elmenhorst et al. 2012). Activity in the reference region (cerebellum) was subtracted from cortical activity and divided by plasma activity.

Serotonin Transporter
Modeling of 5-HTT binding was similar to that in previous publications . In brief, the MRTM2 was applied with thalamus serving as high-uptake and cerebellar gray matter as reference region to obtain the binding potential (BP ND ; Parsey et al. 2006).
Monoamine Oxidase A Arterial blood sampling was performed to obtain arterial input functions after correction for the presence of radioactive metabolites using high-performance liquid chromatography. Arterial input functions were derived from the kinetic modeling tool implemented in PMOD 3.509 (PMOD Technologies Ltd., Zurich, Switzerland; http://www.pmod.com). The arterial input functions and cortical activity were processed using Logan plots to obtain volume of distribution (V T ) as the outcome measure for MAO-A binding (Logan et al. 1990;Ginovart et al. 2006).

Clustering Method
For each protein, data was averaged across subjects resulting in population-based maps for the distribution of 5-HT 1A receptors, 5-HT 2A receptors, MAO-A and 5-HTT on the cortical surface. In order to avoid weighting in the parcellation procedure owing to the differently scaled outcome measures, binding data for each protein was z-scored, which resulted in values with a mean of 0 and a standard deviation of 1. The standardized mean maps for the four proteins were used as input for the k-means clustering algorithm. The initial starting values were chosen according to the K-means++ algorithm (Arthur and Vassilvitskii 2007). Clustering was performed with a squared Euclidean distance measure. The solution with the lowest within-cluster sums of point-to-centroid distances was selected from 50 replicates. To obtain a robust estimate of the optimal number of clusters, bootstrapping was performed using subsets of the PET data. One hundred subsets were generated by randomly selecting five subjects for each molecular target. For each subset, the optimal number of clusters was determined by the silhouette criterion using squared Euclidean distance (Rousseeuw 1987). Possible K were evaluated from 4 (corresponding to the number of proteins used for clustering) to a limit of 20, due to the tendency of the criterion values to decrease with an increasing number of clusters. The K indicated for the majority of subsets in bootstrapping was selected as the optimal solution for the entire dataset. While selection of the optimal number of clusters was performed on a combined dataset of both hemispheres for computational reasons, clustering of the entire dataset was performed both for each hemisphere separately and on a combined dataset of both hemispheres.
To evaluate the resulting clusters, we examined the overlap of parcellation based on serotonergic molecules with those of the population-average, landmark-and surface-based version of the Brodmann atlas (Van Essen 2005), as well as a cortical parcellation based on major functional connectivity networks (Yeo et al. 2011).

Results
Cortical protein maps generated for 5-HT 1A receptors, 5-HT 2A receptors, MAO-A and 5-HTT were in good agreement with distributions reported in the literature (Figure 1, Table 2) (Varnas et al. 2004;Savli et al. 2012;Beliveau et al. 2017). The silhouette criterion indicated K = 5 as the optimal solution in 83 out of 100 subsets during bootstrapping. As results for combined and separate analysis of hemispheres were highly similar ( Supplementary  Figure 1), the subsequent paragraphs only refer to results of the combined analysis. The topology and molecular profiles of clusters are shown in Figure 2 and Table 2. The comparison with other cortical parcellation methods is shown in Figure 3 and Supplementary Table 1. Supplementary Figure 2 displays multiple perspectives of clustering results. Clustering results are available for download in NIfTI and FreeSurfer (.annot) file formats at http://www.meduniwien.ac.at/neuroimaging/parcellation.html. Properties of clusters are summarized in the following paragraphs. Numbering of the clusters is arbitrary.

Cluster 1
Cluster 1 has the lowest binding of MAO-A, 5-HT 1A , and 5-HT 2A receptors. The majority of its area is located on the pre-and postcentral gyri and the paracentral lobule and extends to the superior parietal gyrus. Another part can be found adjacent to the corpus callosum extending from the caudal anterior cingulate cortex to the retrosplenial isthmus cingulate cortex and parahippocampal gyrus (Brodmann areas (BA) 26-29 and 35). Almost two-thirds of this cluster correspond to the somatomotor network and comprise the majority of BA 1-5, 6, and 7.

Cluster 2
Cluster 2 is characterized by the highest binding of excitatory 5-HT 2A receptors among all clusters and second highest expression of inhibitory 5-HT 1A receptors. The cluster spans 24% of the cortical surface. Two-thirds of the cluster are located in the superior, middle and inferior temporal gyri (BA 20-22 which contain Wernicke's speech area), and extend to adjacent lateral occipital, fusiform (BA 37), angular, supramarginal, inferior parietal gyri (BA 39, 40). The remaining part of the cluster is located in the frontal cortex and comprises the majority of the rostral middle frontal gyrus (BA 9, 10, 46). Approximately 10% of the cluster are located in the medial and lateral orbitofrontal cortex (BA 11) and extend into the medial aspect of the superior frontal gyrus. Cluster 2 covers 30% of the default and 43% of the limbic networks.

Cluster 3
Cluster 3 is the largest cluster and includes 37% of the cortical surface. Binding of 5-HTT is the lowest among all clusters, 5-HT 1A receptors and MAO-A are below average. One half of the cluster is located in the superior and inferior parietal cortex and extends into the precuneus (BA 7), lateral occipital (BA 19), supramarginal, and postcentral gyri. The remainder of the cluster can be found in the superior frontal gyrus (BA 8,9) and covers the majority of the middle and inferior frontal (BA 45-47) gyri, as well as parts of the lateral orbitofrontal cortex and precentral gyrus (BA 6). Cluster 3 contributes 60% each to the dorsal attention and frontoparietal networks, and 40% each to the default mode and visual networks.

Cluster 4
Cluster 4 has the highest binding of MAO-A and the second highest binding of 5-HTT. Binding of 5-HT 2A and 5-HT 1A receptors is slightly above and below average, respectively. The majority of Cluster 4 is located on the medial aspect of the cortex, spans from the rostral (BA 24, 33) to the posterior cingulate cortex (BA 23, 30, 31) and reaches into visual areas including the pre-transverse temporal gyri (BA 41,42). Cluster 4 covers one-third of the visual and one-fifth of the default mode networks.

Cluster 5
The highest binding of 5-HT 1A receptors and 5-HTT can be found in Cluster 5, and binding of MAO-A is ranked second. Outcome measures (binding potentials (BP ND or BP P ) or volume of distribution (V T )) are proportional to absolute density of available protein. A total of 108 subjects were investigated, such that 22-32 individual scans were averaged to obtain each protein map. This cluster is the smallest and encompasses most of the insula (BA 13-16), temporal pole (BA 38) and entorhinal cortex and reaches out to the parahippocampal, fusiform (BA 36), and superior temporal gyri. Furthermore, Cluster 5 has an extension into the posterior parts of the orbitofrontal cortex and includes the subgenual area (BA 25). Cluster 5 contributes to 42% of the limbic network.

Relevance of Molecular Clusters in the Serotonergic System
The current results represent a proof of principle for the parcellation of the cerebral cortex based on the topology of serotonergic molecular targets measured using PET. The convergence of clustering results in separate and combined analysis of left and right hemispheres provides internal validation for the method and indicates that clusters identified are not subject to a pronounced influence of lateralization. While k = 5 is a relatively low number of clusters compared to the number of regions identified by cytoarchitectonical features and the MR based methods discussed above, many of the allocations of cortical areas are consistent in synopsis with other methods. It must be kept in mind that the current clustering approach is merely based on four parameters, i.e., four different proteins involved in the same neurotransmitter system, such that differentiation of regions is achieved with respect to the distribution of these particular molecules. It is likely that the inclusion of PET data on other proteins will result in a finer parcellation of cortical regions. Nevertheless, the unique characteristics of each cluster provide a novel perspective on the nature of serotonergic neuromodulation in the cerebral cortex and hint at the potential utility of the proposed parcellation for the investigation of pharmacological effects in this system. In fact, it may be a sensible approach to tailor regions of interests for the analysis of pharmacological imaging data using clustering based on expression of molecular targets of the studied drugs (Gryglewski et al. 2018). If brain areas within clusters defined in this manner react uniformly to pharmacological intervention, a substantial increase in power can be expected. Multiple comparisons can be avoided in contrast to analyses based on atlases obtained from other parcellation methods which may suggest distinctions between brain areas that are not paralleled by a differential response to the effects of the studied drug. In detail, Cluster 1 is almost exclusively covered by the somatomotor network located at the pre-and postcentral gyrus and it may be argued that the proposed parcellation method failed to identify the indisputable border between primary motor and somatosensory cortices. However, if we consider that Cluster 1 is characterized by the lowest expression of all proteins studied, we may hypothesize that the effect of serotonergic neuromodulation is relatively low in these areas after completion of brain development (D'Amato et al. 1987). This is in line with the lack of effects on brain activation of 5-HT 2A receptor agonists in this area (Kraehenmann et al. 2015). The adjacent Cluster 3 is characterized by the lowest binding of serotonin transporters, the expression of which has been found to be highly associated with extracellular serotonin levels (Dewar et al. 1991). Higher binding of serotonin receptors suggest that serotonergic neuromodulation may be effectuated in secondary and associative regions within Cluster 3, i.e. parts of the dorsal attention, frontoparietal, default mode and visual networks, at lower levels of serotonin or serotonin receptor binding drugs than in Cluster 1. Based on the low expression of 5-HTT, no pronounced effects of selective serotonin reuptake inhibitors (SSRIs) can be expected in clusters 1, 2, and 3. In contrast, Cluster 2 is characterized by relatively higher binding of MAO-A compared to 5-HTT and 30% and 43% of this cluster overlap with the default mode and limbic networks, respectively. These areas include the medial prefrontal cortex, temporoparietal junction, orbitofrontal cortex (Raichle and Snyder 2007;Sheline et al. 2009;Zhu et al. 2012), as well as the dorsolateral prefrontal cortex which is targeted with transcranial magnetic stimulation for the treatment of depression atlas misses the assignment of one area, which has a considerable overlap with the insular cortex. Thus, we summarized the unassigned region with Brodmann area 13-16 according to Brodmann and Garey (2006). (Kolbinger et al. 1995;Pascual-Leone et al. 1996). This may underlie the response of patients to MAO inhibitors who have failed antidepressant treatment with SSRIs and justifies the independent role of the success of a trial with MAO inhibitors in staging of treatment resistant depression (Ruhe et al. 2012). The high expression of MAO-A in Cluster 4, 50% of which belongs to the visual network, may further underlie the increase of visual hallucinations with the addition of MAO inhibitors in the form of harmine (radiolabeled for the quantification of MAO-A in our study) and other harmala alkaloids to dimethyltryptamine in ayahuasca preparations, next to their influence on first-pass metabolism (McKenna et al. 1984;De Araujo et al. 2012). Given their high expression of 5-HT 2A receptors, Clusters 2 and 4 are implicated in the effects of drugs acting as agonists on these receptors, such as psilocybin and LSD, which is in line with psychedelic effects on visual and auditory perception and corresponding imaging studies on glucose metabolism, fMRI parameters and blood flow in these areas which were found to correlate with subjective drug effects (Vollenweider et al. 1997;Kraehenmann et al. 2015;Carhart-Harris et al. 2012. Lastly, with their high expression of 5-HTT Clusters 4 and 5 are most likely to be affected by the application of SSRIs which constitute the first line of treatment for depression and anxiety disorders. The extent of Cluster 4 in the cingulate and Cluster 5 in the parahippocampal and insular cortices includes regions that were historically assigned to the limbic system and are continuously under investigation for their role in emotion processing (Mayberg et al. 1999), empathy (Mitchell et al. 2013), and affective disorders (Montag et al. 2009;Horn et al. 2010;Klumpers et al. 2010). There is a remarkable overlap between these clusters and peak decreases in 5-HT 1A receptor binding following treatment with SSRIs (Spindelegger et al. 2009) or electroconvulsive therapy (Lanzenberger et al. 2013), and increases in gray matter volume following electroconvulsive therapy (Sartorius et al. 2016).

Limitations
Next to the limitation that only four proteins were available for the current analysis, several other issues should be mentioned. Given the resolution of molecular imaging using PET, some areas may be affected by partial volume effects, that is blurring of signal at borders between adjacent regions. While partial volume effects are reduced by surface-based analysis (Greve et al. 2014), some regions may be still affected and result in the identification of clusters in areas of mixed signal. The delineation of Cluster 1 close to the corpus callosum is likely due to a partial volume effect rather than low binding in this area of the cingulate cortex. The localization of proteins on pre-or post-synaptic sites, excitatory or inhibitory neurons cannot be differentiated using PET, but is likely to vary across the cortex and affect the effects of pharmacological modulation. Another issue is the fact that all four proteins could not be quantified in the same subjects, to the end that covariances between these proteins are unknown and probabilistic information on the interindividual variation of the extent of clusters could not be generated. However, the protein maps in this analysis are generated from an adequately sized sample of healthy subjects and should constitute a fair representation of average protein distribution in the healthy adult population. Nevertheless, the study of the influence of demographic variables on clustering results cannot be performed with this sample size. The use of a different scanner for 5-HT 2A receptor imaging may reduce the conformity of resulting PET data, which may have been accommodated by the uniform quantification pipeline, smoothing and standardization of the data. The focus of the current work was parcellation of the cortex. Clustering of cerebellum was not performed based on its negligible concentration of three out of four proteins studied, for the quantification of which it serves as reference region. Given the lack of a validated quantification pipeline which allows for integration of cortical surface-based outcome measures with volumetric subcortical regions, these would have to be processed separately. Furthermore, an optimized normalization pipeline for registration of PET data obtained from subcortical regions using different tracers is lacking. Lastly, at the present time we could only provide support for the utility of our method with reference to previous research in pharmacological imaging. Direct application of the defined clusters for analysis of imaging data is still pending.

Concluding Remarks
Definition of brain regions using clustering of molecular imaging data results in a parcellation of the cerebral cortex which is of high explanatory value for the appraisal of the effects of psychotropic drugs. The application of parcels defined in this manner for the analysis of pharmacological imaging data a priori will increase sensitivity and power of these investigations using, e.g. functional magnetic resonance imaging, which is a pressing issue in the field. Therefore, the proposed method constitutes a highly promising approach towards integration of multimodal imaging data for research and development in neuropharmacology and psychiatry.

Supplementary Material
Supplementary material is available at Cerebral Cortex online.
Funding data on the cortical surface, a tool essential for processing of the study results presented here. This scientific project was performed with the support of the Medical Imaging Cluster of the Medical University of Vienna. Conflict of Interest: With relevance to this work there is no conflict of interest to declare.