Brain tissues have single-voxel signatures in multi-spectral MRI

Since the seminal works by Brodmann and contemporaries, it is well-known that different brain regions exhibit unique cytoarchitectonic and myeloarchitectonic features. Transferring the approach of classifying brain tissues - and other tissues - based on their intrinsic features to the realm of magnetic resonance (MR) is a longstanding endeavor. In the 1990s, atlas-based segmentation replaced earlier multi-spectral classification approaches because of the large overlap between the class distributions. Here, we explored the feasibility of performing global brain classification based on intrinsic MR features, and used several technological advances: Ultra-high field MRI, q-space trajectory diffusion imaging revealing voxel-intrinsic diffusion properties, chemical exchange saturation transfer and semi-solid magnetization transfer imaging as a marker of myelination and neurochemistry, and current neural network architectures to analyze the data. In particular, we used the raw image data as well to increase the number of input features. We found that a global brain classification of roughly 97 brain regions was feasible with gross classification accuracy of 60%; and that mapping from voxel-intrinsic MR data to the brain region to which the data belongs is possible. This indicates the presence of unique MR signals of different brain regions, similar to their cytoarchitectonic and myeloarchitectonic fingerprints.


Introduction
Since the seminal works by Brodmann (1909) and Vogt and Vogt (1919 ), it is well-known that different brain regions exhibit unique cytoarchitectonic and myeloarchitectonic features. Transferring the approach of classifying brain tissues -and other tissues -based on their intrinsic features to the realm of magnetic resonance (MR) is a longstanding endeavor. The idea of classifying tissues based on their T1 and T2 relaxation times can be traced back ( Chang et al., 1972 ;Damadian, 1971 ;Iijima et al., 1972 ;Odeblad and Lindstrom, 1955 ;Weisman et al., 1972 ) to the period before the advent of MR imaging (MRI) ( Lauterbur, 1973 ). In fact, Lauterbur used it to motivate MRI ( Lauterbur, 1973 ), and rightly so; the high soft tissue contrast arising from the disperse T1 and T2 times found in the human body is a cornerstone of current radiology that often uses relaxation time weighted MR images.
For this reason, atlas-based approaches that use spatial information to reduce the number of potential tissue classes at a given position were introduced with great success and largely facilitated the classification problem ( Collins et al., 1994 ;Mazziotta et al., 1995 ;Thompson et al., 1997 ). We can do no better than cite from the landmark paper by Fischl et al. (2002 ) published in 2002 to describe this situation: "…it is apparent why no global classification scheme can successfully distinguish structures from each other based only on intensity information -there is far too much overlap between the class distributions (even cortical gray matter and white matter overlap by more than 12%). While adding additional MRI sequences with differing contrast properties or different imaging modalities entirely can help separate the class distributions, spatial information is still required to disambiguate the classification problem. ".
Nonetheless, the approach of classifying tissues based on intrinsic MR properties was pursued after 2002. The focus was on classification problems involving only a few classes and on using additional contrasts. A prominent example is the detection of prostate cancer, a task supported by the upcoming availability of diffusion-weighted imaging sequences ( Artan et al., 2009 ;Artan and Yetik, 2012 ;Chan et al., 2003 ;Chen et al., 2018a ;Freitag et al., 2016 ;Langer et al., 2009 ;Laun et al., 2011 ;Ozer et al., 2010 ).
In the present study, we explore the feasibility of performing global brain classification based on intrinsic MR features. To this end, we use several technological advances. First, we used a newest-generation 7 Tesla (7 T) scanner that yields a high signal-to-noise ratio (SNR) and increased contrast-to-noise ratio (CNR) for many MR contrasts ( Ladd et al., 2018 ). Second, we used a novel kind of diffusion-encoding sequence within the q-space trajectory imaging (QTI) framework ( Westin et al., 2016 ). This allows us to not only measure voxel-averaged diffusion metrics, but also to assess the variance of diffusion tensors within a voxel, which we assumed to be relevant in many cortical gray matter regions featuring more than one dominant fiber orientation ( Aggarwal et al., 2015 ;Calamante et al., 2018 ;Kleinnijenhuis et al., 2015 ;Kleinnijenhuis et al., 2013 ;Leuze et al., 2014 ;McNab et al., 2013 ). Third, we used a chemical exchange saturation transfer (CEST) sequence. We deemed the ensuing obtainable magnetization transfer (MT) contrast as a suitable marker of myelination, which can be used to distinguish and segment different cortical regions ( Cohen-Adad et al., 2012 ;Dinse et al., 2015 ;Glasser et al., 2016 ;Glasser and Van Essen, 2011 ;Haast et al., 2016 ;Mangeat et al., 2015 ). Moreover, ultra-high field CEST is rich in information on several chemically relevant contributions such as pH ( Zhou et al., 2003 ), glutamate ( Cai et al., 2012 ), phosphocreatine ( Chen et al., 2020 ), protein content ( Jones et al., 2013 ;Zaiss et al., 2018b ) and structure ( Goerke et al., 2015 ), as well as lipids ( Görke et al., 2017 ). Thus, we believed that the CEST-derived information on the presence of chemical groups might be related to reported differences in neurochemistry ( Zilles et al., 2002 ). Finally, we used a current feed-forward neural network architecture to analyze the data; in particular, we use the raw data as well to increase the number of input features similarly as in ( Deoni and Jones, 2006 ;Roberts et al., 2017 ). We show that global brain classification thus becomes feasible.

Participants
We recruited 38 participants: eight young men (mean age 24.9 ± 1.6 years) and 18 healthy older participants (mean age 58.6 ± 6.5 years, 4 women and 14 men), and 12 participants with idiopathic Parkinson's disease (PD) (mean age 58.4 ± 8.2 years, 2 women and 10 men) from the Erlangen Movement Disorders Clinic, for an ongoing biomarker study. Potential disease-related brain changes were deemed not relevant because the patients were in the earliest disease stage. All participants provided written informed consent to participate in the study, which was approved by the local ethics committee.

MRI acquisition and processing
MRI acquisition and processing scans were acquired using a 7T MRI scanner (Magnetom Terra, Siemens Healthineers AG, Erlangen) with a 32-channel receive head coil and 8-channel parallel-transmit coil (ptx).
To generate susceptibility maps of each echo time, phase maps were unwrapped using a Laplacian-based phase unwrapping algorithm . A brain mask was generated on the magnitudes using the high definition brain extraction tool created by ( Isensee et al., 2019 ). Subsequently, the background field was removed using the sophisticated harmonic artifact reduction for phase data (V-SHARP) with varying spherical kernel sizes ( Wu et al., 2012 ). Finally, susceptibility maps were calculated by using STAR-QSM ( Wei et al., 2015 ). An echotime averaged quantitative susceptibility map (QSM) was created using the squared echo-time and the squared signal magnitude as a weight ( Chen et al., 2018b ). Susceptibility map-weighted images (SMWI) were created by multiplying a susceptibility mask four times and subsequently multiplying it with the magnitude of the first echo time ( Gho et al., 2014 ).
The SMWI/QSM data were used only to define the gold standard segmentation, but not in the automated classification approach.

Coregistration
CEST and SMWI were coregistered with the structural MPRAGE images using the FSL linear registration tool (FLIRT) with a betweenmodality functions correlation ratio as the cost-function ( Jenkinson and Smith, 2001 ). QTI was coregistered with MPRAGE with a boundarybased registration approach  ) using FLIRT and the FSL automated segmentation tool (FAST) ( Zhang et al., 2001 ). A single blinded rater visually checked for adequate coregistration quality.

Segmentation
A segmentation of each brain into 102 anatomical regions (c.f. Fig. 1 ) was achieved by merging the segmentations obtained with three different approaches. First, to obtain a detailed gray matter segmentation including the hippocampus, automated cortical reconstruction and volumetric segmentation of subcortical structures were performed with the FreeSurfer image analysis suite ( https://surfer.nmr.mgh.harvard.edu , Fischl, 2012 ). The upfront SPM bias-field correction of the ptx-MPRAGE and the high resolution FreeSurfer option were used. Second, the white matter was segmented by adjacency to cortical parcellation using builtin FreeSurfer methods as described in Salat et al. (2009 ). Third, manual segmentation of four brain regions was performed for each participant in an axial orientation using MITK (Medical Imaging Interaction Toolkit, http://www.mitk.org ): the corona radiata and internal capsule were segmented in axial orientation of the MPRAGE data guided by pre-existing FreeSurfer segmentation. The substantia nigra and nucleus ruber were segmented in an axial orientation of the SMWI data irre- Each voxel results in a 341 dimensional vector, the spectrum (example from the precentral cortex in the middle), which is fed to a feed-forward neural network in order to determine the tissue class (bottom). Additionally, the saliency of the network for this voxel is shown. spective of pre-existing FreeSurfer segmentation. The left-right paired regions were merged; higher level regions such as the entire brain cortex were omitted.

Classification using all datasets except for one fixed test dataset
The classification approach is visualized in Fig. 2 . A dataset was created from all participants except for one young man (the test participant, he was selected because he was last in alphabetical order). For each participant individually, 10% of the voxels of each anatomical region were randomly chosen from the MPRAGE space. For each chosen voxel, the local 15 QTI parameters, 210 raw b -tensor values, local four Lorentzian CEST amplitude parameters and 112 raw Z-spectrum values were extracted from the images, coregistered and interpolated to MPRAGE space and resolution, and were saved to a 2D array with 6 million rows and 341 columns. We call each of the row vectors a spectrum. The corresponding anatomical region was stored in a one-hot-encoded 2-dimensional array with 6 million rows and 102 columns. One-hotencoded refers to an array with just zeros, except for a 1 at the class position. Some QTI voxels were excluded due to different CEST and QTI slab thickness (4.2 cm and 6.48 cm, respectively). The two datasets were then shuffled, split, and normalized to zero means and unit variance and no longer contained spatial or neighboring information. We defined a TensorFlow Keras ( Abadi et al., 2015 ) dense neural network (DNN) with a linear stack of four fully connected feed-forward hidden layers with 4096 neurons each, to predict one out of the 102 one-hot-encoded anatomical brain regions. The model had around 52 million trainable parameters and a 50% dropout after each hidden layer during training. The activation functions were rectified linear unit (ReLU) for hidden layers and softmax for the output layer ( Goodfellow et al., 2016 ), loss measure was categorical cross-entropy plus an L2-regularization term for each layer ( = 10 − 5 ), optimizer algorithm was adaptive moments as proposed in Kingma and Ba (2014 ) with learning rate restricted to 10 − 5 , batch size 128, number of epochs 30.
After training, the DNN was used to perform voxel-wise prediction for the test participant. The gross accuracy, defined as number of voxels classified correctly divided by total number of voxels, was calculated.
To investigate the operating principle of DNN classification, we computed first the normalized features of brain tissues averaged over all Fig. 3. Overview of image data (double column). Parameter maps from the test participant, as fed to the neural network (normalized to zero means and unit variance). The maps were derived from QTI except the last four, which were derived from CEST imaging. study participant, and second, the saliency vectors ( Simonyan et al., 2014 ) averaged over the respective regions of the test participant. Five segments were omitted because they were not present in the test dataset (denoted with the index 0 in Fig. 1 ). This is due to our limited slab thickness and individual differences in brain anatomy. Saliency of a neural network input vector highlights features with the biggest influence on the output, and is defined as the local gradient for changing the result of the last layer, output input . Saliency was computed using the keras-vis package ( https://github.com/raghakot/keras-vis by Kotikalapudi, Raghavendra and contributors).
Additionally, the training was performed with subsets of the data: QTI parameters; QTI parameters and diffusion-weighted images; CEST parameters; CEST parameters and z-spectra; QTI and CEST parameters. The respective cross accuracies were calculated.

Classification with cross-validation
Additionally, we performed a 38-fold cross-validation in order to assess group statistics. The procedure was identical to Section 2.5.1 except for that the training time was reduced (10 instead of 30 epochs). In each pass, the DNN weights were reinitialized, and one participant was excluded from the training dataset. After fitting, the DNN was used to perform single-voxel classification in the excluded participant and the gross accuracy was calculated. The Pearson correlation coefficient r was computed for the correlation of gross accuracy and participant age. Spearman's rank correlation coefficient was computed for the correlation of gross accuracy and disease status (1 = healthy, 2 = early PD). For comparison with conventional structural imaging contrasts, we performed an analogous 38-fold cross-validation using only single-voxel T1 MPRAGE and T2 data (from QTI with b = 0 s/mm 2 ).

Classification using only single participant data as training data
In order to explore the feasibility of performing the classification with a small dataset, the network was trained again using just one par-ticipant for the training. Then the gross classification accuracy was calculated for the test participant. Fig. 3 shows a selection of the MR contrasts used. Fig. 4 shows the segmentation results for the test subject obtained with our gold-standard segmentation described in Section 2.4 ( Fig. 4 a  top, Fig. 4 c) and with the neural network-based classification using only voxel-intrinsic information ( Fig. 4 a bottom and Fig. 4 d, using datasets for the training as described in Section 2.5.1 ). Additional MPRAGE is shown in Fig. 4 b. The accuracy from 38-fold cross validation is shown in Fig. 4 e as grayscale intensities. Note that the borders of the thalamus, internal capsule, and basal ganglia are even smoother than the targeted gold-standard segmentation; although cortical boundaries are more scattered (c.f. Fig. 4 a). Fig. 5 shows the class-wise probability maps including the midbrain nuclei. Accurate detection of small structures in the midbrain like substantia nigra and nucleus ruber were particularly surprising. The displayed segments were selected as representative examples from all parts of the brain. All 102 softmax probability images can be viewed in the folder DEMO/SM_NII following the online link in the Data Availability statement.

Image data and segmentation results
The confusion matrix of classification is shown in Fig. 6 for the crossvalidation evaluation using only T1-und T2-weighted data ( Fig. 6 a) and for using the complete spectrum ( Fig. 6 b). In Fig. 6 a, it is apparent that certain predicted labels are chosen for a wide range of true labels reflecting that T1-and T2-weighted data have too large an overlap of class distributions. Using the full spectrum, however, yields more robust results: The segmentation accuracy is generally higher in unique nuclei like the pallidum and unique white matter structures like the limbs of the inter- Fig. 4. Visualization of the segmentation (double column). a) A slice from the test participant is classified using single-voxel information only (bottom) and compared to the gold-standard segmentation in the image domain (top). b) MPRAGE images from the test participant. c) FreeSurfer segmentation from the test participant. d) Single-voxel classification results from the test participant. e) Accuracy from 38-fold cross validation for all participants as grayscale intensities. For color-coding of regions see Fig. 2 .
nal capsule compared to the numerous cortex and white matter regions (c.f. also Fig. 1 ). Fig. 6 b shows that there are confusions between cortical fields and white matter regions, respectively. Furthermore, there are systematic confusions between each cortical field ant its subjacent white matter. Within cortical fields, the medial orbitofrontal and the insular cortex were most accurately detected.
The gross accuracies for the test subject are stated in Table 1 .

Origin of the segmentation accuracy and of confusions
The confusions matched with explicit microstructural and chemical traits of brain regions, which showed highly region specific patterns in brain nuclei and only subtly distinct patterns in the cortical fields and white matter areas ( Fig. 7 a; Fig. 7 uses the evaluation from Section 2.5.1 ). On average, each brain region has a unique signature of QTI and CEST parameters computed from the image raw data ( Fig. 7 a), which are just a subset of the input spectrum that also contained the raw image data.

Fig. 5.
Single-voxel softmax probabilities of nine tissues, results for the test participant (single column). Output of the softmax probability layer of the fitted neural network for each voxel of different brain slices of the validation participant as grayscale intensity with brain outline for orientation. a) ctxcuneus; b) ctx-insula; c) ctx-parsorbitalis; d) wm-cuneus; e) wm-insula; f) post. internal capsule; g) putamen; h) substantia nigra; i) nucleus ruber. For terminology please refer to Fig. 1 . Fig. 7 b shows the saliency vector for every region of the test dataset. For the precentral cortex ( Fig. 1 and region 48, black arrowheads in Fig. 7 b), the DNN relied predominantly on magnetization transfer (MT ( Jones et al., 2013 ), semisolid large macromolecular proteins and lipids) and very high offsets and the water pool resonance. Concordant saliency peaks at − 4 ppm in both z -spectra shows that the DNN has discovered offset domains of the intramolecular NOE as a feature (NOE, Jones et al., 2013 , soluble proteins and lipids, spectrum number 240 and 296). The peak at − 1.5 ppm in the 0.7 μT presaturated z -spectrum (spectrum number 247) suggested additional usage of the phosphocholine CEST signal ( Zaiss et al., 2018b ). The overall saliency vector ( Fig. 7 c) confirms that the DNN is generally reliant on computed CEST and QTI parameters; and high offset, NOE and water pool raw z -spectrum domains; and raw b -tensors with zero b -value (T2-weighted). For specific regions such as the corpus callosum (region 22-26, empty arrowheads in Fig. 7 b), it also depends on raw linear and planar b -tensors. It hardly factored in the amide offset domains and raw strong planar and spherical b -tensors. Nevertheless, it considered parameters calculated from those (e.g., μFA). Fig. 8 shows the gross accuracy obtained in the cross-validation training, which was 51% on average and did not reveal a strong age dependency ( r = − 0.26) nor relevant dependency on disease status (healthy vs early PD, = −0 . 29 ). Presumably due to the reduced epochs, accuracy in the test participant decreased to 54%. Comparison to T1 and T2 based classification revealed improved tissue separation compared to conventional structural imaging. Fig. 9 shows results for the single-participant training, the gross accuracy was 35%, indicating that more training data improves accuracy.

General discussion
In this work, we explored the feasibility of performing a global brain classification based on voxel-intrinsic multispectral MR features. Different from earlier works on global brain classifications, we extended the Fig. 6. Single-voxel classification accuracy across brain tissues (double column), results from all participants ( Section 2.5.2 ). Row normalized confusion matrices of single-voxel classification (left: T1w and T2w, right: CEST and QTI) and true labels as determined by image segmentation. Labels 28-59 corresponds to cortex areas, labels 60-91 to white matter regions, the two diagonal parallel lines correspond to respective cortex-white matter pairs. For detailed labels please refer to Fig. 1 . input data space by novel contrasts (QTI and CEST) and by including raw images, similarly to Roberts et al. (2017 ) and Deoni and Jones (2006 ). We found gross classification accuracy of 60% in the test participant without performing any advanced postprocessing of the classification results ( Geman and Geman, 1984 ).
With the single-voxel classification approach, the spatial coherence of prediction results in adjacent voxels serves as an inherent metric of reliability (such as in Fig. 4 ). However, current state-of-art atlasbased classifications approaches outperform the accuracy we observed; moreover, they do not rely on advanced, not widely available MR sequences such as the ones we used ( Ashburner and Friston, 2005 ;Avants et al., 2011 ;Fischl et al., 2002 ;Glasser et al., 2016 ;Iglesias and Sabuncu, 2015 ). Our results indicate that adding advanced contrasts can improve the classification performance of atlas-based approaches, naturally with the drawback of limited availability to the general community and potentially increased scan time. However, we see the merit of our findings to a lesser extent in the potential avenue towards better classification results. The more striking finding is that decent gross classification accuracy can be achieved by using voxel-intrinsic information only. Bearing in mind that different brain tissues exhibit rather unique cytoarchitectonic and myeloarchitectonic characters, it is more than tempting -and presumably justified to some degree -to consider our MR-based patterns as analoguous to these classical tissue fingerprints. Similarly as for the prostate cancer applications mentioned in the introduction, such a voxel by voxel approach might become most important in application to pathologies, where atlas-based approaches face the challenge that more tissue classes can be present at each point in space. Imaging of stroke, tumors, and degenerative diseases such as Alzheimer's disease might represent important application fields.
Naturally, the question arises of whether there are confounding factors that imprint spatial information to some degree into the single-voxel spectrums. Potential sources of such imprints can be B 0 or B 1 inhomogeneities, which increase at 7T compared to lower field strengths ( Ladd et al., 2018 ), eddy current induced image distortions that can cause severe artifacts at interfaces between different tissue types, or of surface coil flair ( Hayes and Axel, 1985 ). In the CEST preparation, we addressed the B 1 inhomogeneity by means of the MIMOSA approach ( Liebert et al., 2019 ;Orzada et al., 2010 ). In addition, B 0 and B 1 correc-tion was performed employing acquired field maps . We considered the B 0 and B 1 inhomogeneity to be less relevant for the diffusion data and did thus not perform a correction, which may be inaccurate in lower brain regions. Coil surface flare may not be corrected as easily at 7T as at lower field strengths because there is no body coil that provides homogenous field distribution information and thus we could not use the vendor-provided prescan normalize option. Interpolating from different original CEST and QTI resolutions to the high MPRAGE resolution adds environment information to the single-voxel spectra. A remaining effect thus cannot be ruled out. Its effect size is rather difficult to assess. One assessment strategy would be to evaluate the dependency of the gross classification accuracy on the experimental settings. Data could be acquired at lower field strengths which experiences reduced B 0 and B 1 field inhomogeneities to assess their effect, although the SNR and CNR would also change, which renders the interpretation difficult. Full-fledged ptx-sequences that would minimize the B 1 field inhomogeneity could be used ( Deniz, 2019 ;Wu et al., 2018 ) and local shimming setups could be used to minimize the B 0 inhomogeneities ( Aghaeifar et al., 2020 ). In the present study, these techniques were not available and thus we considered the extended investigation of these effects to be beyond our scope. At least, we address overfitting due to interpolation by splitting our dataset not randomly into training and testing subsets, but in a participant aware fashion. Reassuringly, Nagy et al. similarly found that single-voxel features derived from diffusion data can be used for classifying cortical brain regions ( Nagy et al., 2013 ). They used 3 T data, which suffers to a much lesser degree from B 1 and B 0 inhomogeneities.
Neural networks are a powerful tool for retrieving subtle signatures of a tissue, but elucidating their classification strategy is notoriously difficult Sahiner et al., 2019 ). Potentially for that reason, many of the tissue classification papers mentioned in the introduction do not perform any analysis in this regard. But, for example, Özkan et al. performed an evaluation and visualization of the decision boundaries, which is helpful with their 4-component input vector ( Özkan et al., 1993 ), but may be less suited for larger spectra. We tried a 2-fold strategy to understand the classification strategy of our neural network. First, mapping the dominant features of the brain regions showed that distinct regions exhibit different patterns, which are Fig. 7. Origin of the segmentation accuracy and of confusions (double column). a) Normalized microstructural and chemical features of brain tissues, averaged over all voxels from all study participants. b) Saliency vectors for brain tissues of the test subject. Empty arrowheads pointing to the corpus callosum, which features elevated b-tensor saliency. Filled arrowheads pointing to elevated MT saliency in the precentral cortex. c) Average saliency vector for all voxels from the test subject.
rather unique for most brain regions except for the cortical region, which were the most difficult to classify. Second, we evaluated the saliency ( Sahiner et al., 2019 ;Simonyan et al., 2014 ) and found that the network relied more on evaluated parameters (such as μFA and NOE) than on the raw data (diffusion-weighted images, CEST-weighted images), although the raw data essentially contains the information of these computed parameters. Potentially the network relies on them because they are more readily available, or because of a signal-to-noise effect; as the quantita-tive parameters contain the averaged information from many noisy raw data points. Only for certain brain regions, raw b -tensors were highly salient.
Many recent deep learning segmentation approaches were built on convolutional neural networks (CNN) ( Goodfellow et al., 2016 ). Mostly morphological image data were used (e.g. Henschel et al., 2020 ), but successful segmentations based on diffusion data have been reported ( Li et al., 2020 ;Wasserthal et al., 2018 ). Another common approach for tissue classification is radiomics ( Parekh and Jacobs, 2020 ;Rizzo et al., 2018 ). Unlike our approach, these approaches also use spatial information. A disadvantage of these approaches is the often large demand on the amount of available data ( Yip and Aerts, 2016 ). In contrast, the voxel-wise approach we used works a bit even if the training involves only a single participant, for which we achieved 35% gross validation accuracy. This is relevant considering that many methodological studies in MRI are performed with sample sizes of < 10 and may be explained by the fact that even measuring a single participant may yield millions of data points ( Hanspach et al. (2021 ). Indeed, successful single-voxel neural network training of CEST data with few datasets has been reported ( Chen et al., 2020 ;Glang et al., 2020 ;Zaiss et al., 2019 ).
One drawback of using additional advanced MR contrasts is the generally increased scan time needed to sample the data. Here, we did not use acceleration approaches such as MR fingerprinting (MRF) ( Ma et al., 2013 ), but it might be crucial for future applications if the aim is to further extend the input data space. Extensions to CEST MRF already exist that infer exchange rates and concentrations via deep learning ( Kim et al., 2020 ;Perlman et al., 2020 ).
The FreeSurfer segmentation is not perfect and may suffer from missegmentation at the border zones between different regions. This error is especially prominent in regions with high spatial variance approximating the image resolution like in the cortex. Indeed, we observed the lowest accuracy of our method in the segmentation of highly folded cortical areas, while the accuracy in homogeneous regions like thalamus and basal ganglia was very high. One possible extension would have been to test a range of different segmentation tools. We found that the SPM upfront bias field correction generally worked very well with our 7 T data and we were not forced to use particular 7 T adaptations for the FreeSurfer segmentation.
A general limitation is that segmentations without direct histological correlation just shows implicit, imprecise spatial distribution of microstructural features. Another interesting extension of our approach would therefore be to map not solely to brain segmentation, but also to tissue fingerprints obtained ex vivo. For example, one could map to neurochemistry fingerprints ( Zilles et al., 2002 ) or features derived from cytoarchitectonic histological data.

Strengths and limitations
A strength of the approach is that even data of a single subject yields decent classification results, which makes it easy to explore the benefits of adding further contrasts. Further strengths include that the spatial coherence displays reliability, the potential to straightforwardly extend it to other fields of application, and the potential to link to chemical and microstructural tissue properties.
A limitation of the present study is that the generalizability with respect to using data acquired with further scanners was not explored. We anticipate that in particular field strength related alterations would limit the generalizability because the CEST contrast changes ( Ladd et al., 2018 ); and because we observed also a field strength dependency of QTI parameters in previous works ( Endt et al., 2019 ;Martin, 2019 ). Another general limitation is that due to SAR limitations, the power of our CEST scans is not high enough to obtain signals from fast exchanging molecular species. As mentioned above, further limitations include the demand for advanced, partly not widely available MR sequences, the currently rather long scan time, the rather indirect connection to underlying microstructural and chemical parameters, and the yet to explore influence of potential confounding factors such as B 0 or B 1 inhomogeneity.

Conclusion
In conclusion, we show that mapping is possible from voxel-intrinsic MR data to the brain region to which the data belongs. This indicates that unique MR signals of different brain regions exist, similar to their cytoarchitectonic and myeloarchitectonic fingerprints.

Data availability
The preprocessed training dataset, the face and skull stripped data of the test participant, technical specification of features and segments, the trained neural network, and a short demonstration video can be accessed at https://nx8708.your-storageshare.de/s/TWoxPAjyFpPnyAx .

Code availability
Code for inspecting and analyzing the training and test dataset is available at https://github.com/a-german/brain _ voxel .

Author contributions
AG proposed the single-voxel classification and performed MRscanning and data analysis. AG generated the figures with support by FBL. AG and FBL wrote the manuscript with comments from all authors. JW, AG, AD and AM recruited participants. MS determined fitness of the patients. AD and MU organized and supervised all MRI exams. JHe, JHa, JM, TAK, AN, AL and MZ implemented and established the MPRAGE, QSM, QTI, and CEST MR sequences and postprocessing software. AM, AG, MZ, and FBL designed the MR-protocol. JW, AD, MU, MZ, and FBL supervised the work. The present work was performed by AG in fulfillment of the requirements for obtaining the degree "Dr. med. " at Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU).

Funding
This work was supported by the Interdisciplinary Center for Clinical Research Erlangen (IZKF-Erlangen, ELAN DR-17-12-14-1). AG received a MD thesis fellowship by the IZKF-Erlangen. JW is a member of the research training group 2162 "Neurodevelopment and Vulnerability of the Central Nervous System " funded by the DFG-270949263/GRK2162. The work was further supported by a Heisenberg professorship from the DFG to FBL (DFG LA 2804/6-1 and 2804/12-1).

Declaration of Competing Interest
AG was employed by Siemens as a scientific assistant.