Abstract
Local measures of neurotransmitters provide crucial insights into neurobiological changes underlying altered functional connectivity in psychiatric disorders. However, noninvasive neuroimaging techniques such as magnetic resonance spectroscopy (MRS) may cover anatomically and functionally distinct areas, such as p32 and p24 of the pregenual anterior cingulate cortex (pgACC). Here, we aimed to overcome this low spatial specificity of MRS by predicting local glutamate and GABA based on functional characteristics and neuroanatomy in a sample of 88 human participants (35 females), using complementary machine learning approaches. Functional connectivity profiles of pgACC area p32 predicted pgACC glutamate better than chance (R2 = 0.324) and explained more variance compared with area p24 using both elastic net and partial least-squares regression. In contrast, GABA could not be robustly predicted. To summarize, machine learning helps exploit the high resolution of fMRI to improve the interpretation of local neurometabolism. Our augmented multimodal imaging analysis can deliver novel insights into neurobiology by using complementary information.
SIGNIFICANCE STATEMENT Magnetic resonance spectroscopy (MRS) measures local glutamate and GABA noninvasively. However, conventional MRS requires large voxels compared with fMRI, because of its inherently low signal-to-noise ratio. Consequently, a single MRS voxel may cover areas with distinct cytoarchitecture. In the largest multimodal 7 tesla machine learning study to date, we overcome this limitation by capitalizing on the spatial resolution of fMRI to predict local neurotransmitters in the PFC. Critically, we found that prefrontal glutamate could be robustly and exclusively predicted from the functional connectivity fingerprint of one of two anatomically and functionally defined areas that form the pregenual anterior cingulate cortex. Our approach provides greater spatial specificity on neurotransmitter levels, potentially improving the understanding of altered functional connectivity in mental disorders.
Introduction
Since the beginning of the 20th century, the human brain is understood to consist of regions with distinct microarchitecture (Brodmann, 1909; Vogt and Vogt, 1919). Anatomical features, such as cytoarchitecture, myeloarchitecture, and receptorarchitecture distinguish cortical areas and highly constrain the local processing capabilities of a region (Cloutman and Lambon Ralph, 2012; Eickhoff et al., 2015; Zilles and Palomero-Gallagher, 2017; Palomero-Gallagher and Zilles, 2019). Cytoarchitecture also shapes the functional repertoire of a region through specific input and output connections, which embed the region in complex distributed networks (Cloutman and Lambon Ralph, 2012). This mesoscopic functional repertoire can be assessed with modern functional magnetic resonance imaging (fMRI) techniques, but methodological and ethical constraints limit us in assessing its relationship to neurometabolism on the microscale in vivo. Nevertheless, such linking of functional connectivity (FC) to local processing (e.g., the local excitation/inhibition balance) could provide valuable insights into the pathophysiology of psychiatric disorders.
Though in vivo measurements of local metabolism are not possible at a microscopic level, the noninvasive method of proton magnetic resonance spectroscopy (1H-MRS) is commonly used to assess local neurotransmitter levels. Available sequences such as point resolved spectroscopy (PRESS; Bottomley, 1987), stimulated-echo acquisition mode (STEAM; Frahm et al., 1987), and edited sequences such as MEGA-PRESS (Mescher et al., 1998) allow for the detection of glutamate (Glu) and GABA in single voxels in the human brain. MRS measures of these metabolites may be used to better understand FC changes in psychiatric disorders (Horn et al., 2010; Moeller et al., 2016; Allen et al., 2019) or pharmacological challenges (Li et al., 2017), but such results may suffer from limited interpretability because of the low spatial specificity of conventional single-voxel spectroscopy. An MRS voxel may cover several known cytoarchitectonically distinct areas (Duncan et al., 2014). For example, MRS measures of glutamate and GABA follow differential receptor distributions along the anterior cingulate cortex (ACC; Dou et al., 2013; Li et al., 2017). MRS imaging (MRSI), which maps the spatial distributions of metabolites in the brain (Posse et al., 2013) could overcome this limitation. However, MRSI at sufficiently high spatial resolution has impedingly long acquisition times for use in multimodal (patient) studies and requires more expertise on the acquisition and data-processing side compared with single-voxel spectroscopy (Henning, 2018; Nassirpour et al., 2018). While single-voxel spectroscopy is more easily implemented in clinical and multimodal studies, it is not known whether the measured concentrations of glutamate and GABA are representative of the entire MRS acquisition voxel and whether the spatial resolution of conventional MRS could be improved by using more fine-grained weights informed by functional imaging.
The pregenual ACC (pgACC) is a region that has been frequently studied in the MRS literature. Thus far, these studies have not accounted for the considerable cytoarchitectonic and functional variation between areas that comprise this region. The pgACC is part of the default mode network (DMN) of the human brain and has been implicated in the pathophysiology of depression (Salvadore and Singh, 2013). Previous work suggests altered glutamatergic metabolism in the pgACC in patients with major depressive disorder (MDD; Walter et al., 2009; Horn et al., 2010; Colic et al., 2019) and pgACC levels of a marker of glutamatergic metabolism (glutamate + glutamine) were correlated with FC between pgACC and insula in patients but not in healthy control subjects (Horn et al., 2010). Based on its differences in cytoarchitecture and densities of multiple receptors, the pgACC has been divided into six distinct regions: p24a, p24b, pv24c, pd24cd, pd24cv, and p32 (left hemisphere, 2918 ± 447 mm3; right hemisphere, 3107 ± 447 mm3; Palomero-Gallagher et al., 2019). These areas were partly merged into a gyral component (areas p24a and p24b into area p24ab; left hemisphere, 680 ± 271 mm3; right hemisphere, 663 ± 228 mm3) and a sulcal component (areas pv24c, pd24cd, and pd24cv into area p24c; left hemisphere, 603 ± 195 mm3; right hemisphere, 711 ± 92 mm3) for computation of 3D probabilistic maps (Palomero-Gallagher et al., 2008, 2019). Studies in nonhuman primates have shown that homologous areas have distinct structural connectivity patterns (Pandya et al., 1981). In humans, meta-analytic connectivity modeling showed that these areas have largely distinct functional connectivity patterns, with activation in area p32 largely associated with tasks involving “theory of mind” and cognitive regulation of emotion, and areas p24ab and p24c with tasks involving the experience of one's bodily state and action inhibition, respectively (Palomero-Gallagher et al., 2019). In sum, disregarding this heterogeneity during MRS acquisition may hamper the interpretation of links between local metabolism and functional connectivity.
To overcome the problem of low spatial specificity of conventional MRS, we propose a novel, multimodal approach offering a more nuanced prediction of glutamate and GABA in a single MRS voxel. To this end, we used ROIs originating from a voxelwise FC-based (“functional”) parcellation of a pgACC MRS voxel and a cytoarchitectonic (“anatomic”) parcellation of the same region (Palomero-Gallagher et al., 2019). We assessed correspondence between the functional and anatomic parcellations and tested whether the prediction of pgACC glutamate or GABA was improved by parcellating the voxel (Fig. 1). Crucially, FC profiles of p32 but not p24 could robustly predict local glutamatergic metabolism. We explored why glutamate but not GABA levels were significantly predicted from area FC by examining differential associations between GABAergic and glutamatergic gene coexpression and FC. We then addressed the functional relevance of the differential prediction of pgACC glutamate levels by “decoding” partial least-squares regression (PLSR) β-weights. Overall, our results demonstrate that fMRI may improve the spatial specificity of local neurometabolites assessed with MRS.
Overview of the primary analysis pipeline. To improve spatial information on local Glu and GABAergic levels measured with MRS, an MRS ROI was parcellated based on seed-based FC to 132 atlas nodes. From this hierarchical clustering step, two clusters emerged (Extended Data Figure 1-1). Clusterwise FC profiles were used as input into elastic net and partial least-squares regression models to predict GABA+/tCr and Glu/tCr residualized for gray matter proportion in participants' MRS voxels.
Figure 1-1
Dendrogram. This dendrogram shows the results of hierarchical clustering of the MRS ROI based on functional connectivity to 132 target ROIs derived from the CONN atlas. Two main clusters emerged. Peach, Functional p32; yellow, functional p24. Download Figure 1-1, TIF file.
Materials and Methods
Participants
We included 143 healthy participants in this study. Data were pooled from three studies. All participants were screened for prior and current neurologic or psychiatric illness using the German version 5.0.0 of the M.I.N.I. (Mini-International Neuropsychiatric Interview; Ackenheil et al., 1999). All study procedures were approved by the ethical committee of the University of Magdeburg and conformed with the Declaration of Helsinki.
From the 143 measurements, we removed two datasets for which MRS but not fMRI data were available. We then excluded four datasets that did not fulfill MRS quality criteria (see Data acquisition and preprocessing) as well as two additional datasets with insufficient water suppression in the MRS data (identified through visual inspection of the baseline by L.C., M.L., and L.M.). Visual inspection of the MRI data by L.C., L.M., and M.L., with a final decision by M.L., led to the exclusion of 37 further measurements. These datasets had either insufficient coverage of the cerebellum, or striping or ghosting artifacts. Last, several participants participated in more than one MRS acquisition. In these cases, we discarded the measurement that did not meet initial quality criteria or, in case multiple measurements of the same participant were of sufficient quality, we selected the measurement with the best MRS and/or fMRI data quality in advance of the current analysis. This step excluded 10 additional measurements.
The final sample consisted of 88 participants with good-quality resting-state and MRS data (age, 28.81 ± 9.02 years; 35 females) for analysis of Glu/creatine + phosphocreatine (tCr). For GABA+/tCr, one female and four male participants could not be included in the analyses because of Cramer–Rao lower bounds (CRLBs) exceeding the cutoff of 20, leading to a reduced sample size of 83 participants (age, 28.6 ± 8.99 years; 34 females).
Data acquisition and preprocessing
MRS
Ultra-high-field data were acquired on a 7 T MAGNETOM scanner equipped with a 32-channel head coil (Siemens). Before MRS measurements, an MPRAGE T1-weighted scan was acquired. The echo time (TE) was 2.73 ms, repetition time (TR) was 2300 ms, and inversion time was 1050 ms. The in-plane field of view (FOV) was 256 mm. The flip angle was set to 5°. Images were acquired with a bandwidth of 150 Hz/pixel and 0.8 mm isotropic image resolution.
Participants' T1-weighted scans were used for accurate placement of the pgACC voxel, according to an established protocol of anatomic landmarks described in the study by Dou et al. (2013). Briefly, the pgACC voxel (10 × 20 × 15 mm3) was placed in the bilateral pgACC and centered on the sagittal midline to ensure maximal coverage of cingulate gray matter. Automatic shim routines were used to optimize B0 homogeneity. We applied a STEAM sequence with variable-rate selective excitation radio frequency pulse (Dou et al., 2013) with short TE/mixing time (20 ms/10 ms) and TR of 3000 ms. Metabolite spectra were acquired with 128 averages. A single water reference signal was acquired for eddy current correction. The bandwidth was 2800 Hz, and the acquisition time for one image was 731 ms.
MRS data were fitted in LCModel version 6.3.0 (Stephen Provencher; Provencher, 2001). The basis set used for fitting included Cr, Glu, myo-inositol, lactate, N-acetylaspartate (NAA), phosphocholine (PCh), taurine, aspartate, GABA, glutamine, glucose (Glc), alanine, NAA-glutamate (NAAG), phosphocreatine, scyllo-inositol, acetate, succinate, phosphorylethanolamine, glutathione, citrate, and glycerophosphocholine. The analysis window was set to range from 4.0 to 0.6 ppm. Eddy current correction was performed based on the water signal (LCModel parameter DOECC = T) and water suppression was performed (DOWS = T). The attenuation factors of water (ATTH2O) and metabolites (ATTMET) were assumed to be 0.67 and 0.69, respectively. The chemical shift of the singlet used for scaling (Cr) was set to 3.0 ppm. To account for uncertainty in the referencing between in vitro (the simulated basis set) and in vivo measurements, the default SD of shift (DESDSH) and the expectation of 1/T2 (DEEXT2) were set to 0.01 and 12.0, respectively. Additional changes in the SDSH from the new DESDSH were applied for NAA (0.004), NAAG (0.004), Glc (0.025), and PCh (0.025). Lactate, scyllo-inositol, and acetate were omitted (CHOMIT) from the basis set used for fitting.
In all analyses, we used glutamate and GABA as a ratio to total creatine (i.e., tCr). Metabolite values were considered of insufficient quality if the signal-to-noise ratio (SNR) was <20, if the line width [full-width at half-maximum (FWHM)] was >24 Hz or if the CRLB was <20%. In our final sample used in the analyses of glutamate (n = 88), mean ± SD SNR was 5.639 ± 1.907. For FWHM and CRLB, the mean ± SD values were 43.455 ± 6.710 Hz and 3.159 ± 0.933, respectively. For the analyses of GABA (n = 83), the mean ± SD SNR was 43.940 ± 6.421, FWHM was 5.613 ± 1.907 Hz, and CRLB was 8.928 ± 2.556. Mean raw and fitted data are shown in Figure 2. To rule out that effects can be attributed to the denominator tCr, we computed water-scaled (“absolute”) millimolar concentrations of Glu and GABA+ using the procedure outlined in the study by Giapitzakis et al. (2018).
Summary of MRS data quality. This figure shows the mean raw and fitted spectra across all participants (n = 88) in black and red, respectively. The shaded band around the raw data depicts the SD of the raw spectra. Mean residuals (n = 88) are shown at the top of the figure in black. The line plots below the mean raw and fitted spectra represent the mean fits for glutamate (n = 88) and GABA (n = 83).
We corrected metabolite values for voxel gray matter (GM) content. Partial volume corrections are necessary even when using concentrations relative to creatine, because although tCr, Glu, and GABA all have higher concentrations in GM compared with white matter (WM) and CSF, the metabolite distributions are not identical (Nassirpour et al., 2018). Additionally, tCr and Glu are not necessarily correlated with GM proportion to the same degree (Zhang and Shen, 2015). If this is not accounted for, small changes in the denominator may lead to disproportionally high variance in the ratios (Li et al., 2003). To perform this partial volume correction, we segmented participants' T1-weighted images using voxel-based morphometry (VBM) in the CAT12 toolbox for SPM (http://dbm.neuro.uni-jena.de/cat/index.html#VBM; Ashburner and Friston, 2005). MRS voxels were normalized to MNI space using the forward deformation field produced during segmentation and normalization of the structural scans. Participants' normalized gray matter tissue probability map produced by VBM were then masked with their normalized MRS mask, and the percentage of probable voxel GM content was calculated. We constructed linear regression models in which Glu/tCr or GABA+/tCr were predicted from voxel GM content. The residuals of these models represented the Glu/tCr of GABA+/tCr after the influence of the percentage of GM in the voxel has been removed. We performed all further analyses on these residuals.
rs-fMRI acquisition
For the acquisition of resting-state fMRI (rs-fMRI) data, participants were instructed to keep their eyes closed and think of nothing in particular. The acquisition time for this measurement was 13 min and 18 s, and TR/TE were 2800 and 22 ms, respectively. The image resolution was 2 mm isotropic, and the in-plane FOV was 212 mm. The flip angle was 80°. Sixty-two slices were acquired for a total of 280 volumes. In-plane parallel imaging was done with GRAPPA (generalized autocalibrating partially parallel acquisition) image reconstruction (Griswold et al., 2002) acceleration factor 3. The first 10 volumes of resting-state data were discarded to reach steady state.
ROI definition
For functional parcellation, we created an ROI based on participants' pgACC MRS masks. Participant-specific masks would introduce a bias in the functional connectivity profile and could therefore inflate associations with local neurometabolism. For this reason, we created a composite mask. We resampled normalized participants' MRS masks to the space and resolution of functional images (2 × 2 × 2 mm3). From the resulting masks, we created a composite group MRS ROI such that each voxel within the ROI was contained in the normalized MRS mask of at least two participants (i.e., threshold for inclusion of a voxel, >1). In addition, we used a recently published anatomic parcellation of the pgACC as a second, atlas-based ROI parcellation (Palomero-Gallagher et al., 2019). This parcellation consists of maximum probability maps (MPMs) of areas p24ab, p24c, and p32. The delineation of the areas reflects cytoarchitectonic differences and is based on 10 postmortem human samples.
rs-fMRI preprocessing
Preprocessing of resting-state data were done using the CONN toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012). Briefly, images were realigned and unwarped (motion correction) and slice-time corrected. Functional and structural images were then segmented using default tissue probability maps and normalized to MNI coordinates using direct normalization and resampled to 2 mm isotropic voxel size. To take the fullest advantage of the gained spatial resolution as a result of using ultra-high-field strength, we did not apply spatial smoothing. This approach also allowed us to limit the calculation of seed-based FC to the MRS ROI. Further denoising was performed using custom MATLAB scripts. In this step, time series from voxels with either a GM probability of ≥0.35, or from those voxels falling within the MRS ROI, were z-scored, despiked, and subjected to quadratic detrending, after which six motion parameters (estimated during the realignment step) and the mean WM signal were regressed out. We did not perform bandpass filtering of the time series, as most high-frequency fluctuations related to physiological noise were likely to have been removed by regressing out the mean WM signal (Kahnt et al., 2012) and low-frequency drifts were removed in the detrending step (Tanabe et al., 2002).
FC calculation
FC was calculated separately for the group MRS ROI and the anatomic ROI. As seed voxels, we selected only those fMRI voxels in ROIs. To be able to interpret the resulting cluster solution in anatomically informed ways and to reduce the feature space to k closer to our N value (N = 88), we selected the 132 CONN atlas nodes as target ROIs. These include Harvard–Oxford cortical and subcortical ROIs (Desikan et al., 2006) as well as cerebellar ROIs from the AAL (automated anatomical labeling) atlas (Tzourio-Mazoyer et al., 2002). Functional connectivity between seed time series and mean target time series was calculated as the Pearson correlation between the two. This resulted in a three-dimensional connectivity matrix with 1216 (MRS ROI) or 2160 (atlas-based ROI) rows (seeds), 132 columns (target ROIs), and 88 participants in the z-dimension.
Connectivity-based parcellation of the MRS ROI
We parcellated the pgACC MRS ROI into clusters of similar connectivity using resting-state functional connectivity-based parcellation. The aim of this method is to decrease within-cluster distance and to increase between-cluster distance (Eickhoff et al., 2015). First, the seed × ROI × participant matrix was Fisher z-transformed and averaged across participants (Kahnt et al., 2012). We then computed the correlation between the connectivity profiles of every seed. To finally parcellate contained voxels, we created a similarity matrix (Nseed × Nseed) with the Pearson correlation coefficient as the distance measure.
To assess the functional hierarchy within the MRS seed voxel, the resulting similarity matrix was then subjected to hierarchical clustering to cluster functional voxels according to their similarity in terms of whole-brain functional connectivity (Johnson, 1967). A major advantage of hierarchical clustering is that, unlike, for example, the popular k-means algorithm, a hierarchical approach does not require a predefined number of clusters. The dendrogram may be cut at any level, with k + 1 clusters always nested in k clusters (Cloutman and Lambon Ralph, 2012; Eickhoff et al., 2015). We used the “average distance” linkage algorithm for clustering. Based on inspection of the dendrogram and the FC matrices, a two-cluster solution was found to be optimal for the overarching goal of the study. Note that because of the hierarchical clustering algorithm, more fine-grained parcellations could be tested to further localize predictive voxel clusters within their cluster branch. Yet, this would come at the cost of multiple-comparison correction and makes specificity incrementally harder to demonstrate. Thus, we decided to focus solely on the two-cluster solution to demonstrate the utility of the method in principle.
On account of their functional similarity, we averaged the functional connectivity to target regions across all seed voxels for each cluster to reduce the number of features in subsequent analyses. This resulted in separate matrices for each cluster, representing the cluster-to-target FC for each participant. To assess effectivity of functional parcellation, we calculated the within-cluster distance (sum of squares) and compared it to the between-cluster distance (sum of squares). Both within- and between-cluster distances were compared with whole MRS ROI distance to mean FC.
Cluster validity
Studies using tractographic and functional connectivity-based parcellation showed substantial correspondence between their parcellation and probabilistic cytoarchitectonic maps (Blumensath et al., 2013; Wig et al., 2014; Gordon et al., 2016). Others have provided qualitative evidence (visual inspection) for overlap between cytoarchitecture and connectivity based-derived parcellation (Beckmann et al., 2009; Gordon et al., 2016; Balsters et al., 2018).
We compared our functional clusters to the anatomic parcellation of the pgACC. We summed the MPMs of clusters p24ab and p24c to create a mask of area p24. The MPMs of areas p24 and p32 were then resampled to the space of our functional clusters. MPMs were restricted to MRS ROIs. The overlap between the MPMs and the functional clusters was then calculated using Dice coefficients (DCs; Arslan et al., 2018).
Influence of individual participants' voxel placement
To investigate a possible mediating influence of participants' exact voxel placement on neurometabolite prediction, we compared the means of overlap of participants' individual MRS masks with functional p32 (meanDC = 0.357, SDDC = 0.125) and with functional p24 (meanDC = 0.482, SDDC = 0.112; t(87) = −5.108, p < 0.001). The Dice overlap between individual MRS voxels and functional p32 did not correlate with local Glu/tCr (t(86) = −0.952, p = 0.344) or GABA+/tCr (t(81) = 1.311, p = 0.193), confirming that individual participants' voxel placement did not influence our results.
Statistical analyses
Correlational analyses and t tests were performed in R (version 3.5.0) with the RStudio IDE (version 1.0.136). All other statistical analyses were performed in MATLAB 2019a. The α was set to 0.05, two-tailed unless noted otherwise. Because the methods—elastic net (EN) and partial least-squares regression—are not independent, but complementary, we did not correct for multiple comparisons.
Demographics
To assess the influence of possible confounders, for both Glu/tCr and GABA+/tCr we calculated Pearson's correlations with age and gray matter proportions. To test for a difference between Glu/tCr and GABA+/tCr between male and female participants, we performed a Welch's two-sample t test.
Characterization of functional ROI FC profiles
We aimed to characterize the functional connectivity profiles of the functional areas p32 and p24. To this end, we calculated participant-wise mean time series for each area, and calculated Pearson's correlations with participant-wise whole-brain time series. We performed a paired t test on both areas to test which brain areas were significantly more functionally connected to either functional area, compared with the other functional area. This analysis was conducted using an implementation of threshold free cluster enhancement (TFCE) in MATLAB (https://github.com/markallenthornton/MatlabTFCE). The threshold for significance was set at 0.05, TFCE corrected. For visualization purposes, we repeated the paired t test in MATLAB and plotted the unthresholded t values (p32 > p24) by Yeo network (see Fig. 4E).
Neurosynth decoding
To further illustrate the difference between functional p32 and p24 FC, we “decoded” average FC difference maps (p32 – p24) using data from the Neurosynth framework (Yarkoni et al., 2011). The Neurosynth framework comprises neuroimaging data and text extracted from 14,371 fMRI studies. The decoder toolkit implemented within this framework allows for decoding cognitive states from a given (activation) map (Rubin et al., 2017). Using this toolkit, we correlated the average FC difference map with all cognitive state maps available in the Neurosynth database (release 0.7). Each voxel in a cognitive state posterior probability map reflects the likelihood that a cognitive state term is used in a study if the voxel is activated (Rubin et al., 2017; Quintana et al., 2019). The top 40 terms correlated with the difference map are displayed in a word cloud (see Fig. 9F).
Partial least-squares regression
To predict Glu/tCr and GABA+/tCr from functional and anatomic p32 and p24 FC, we used PLSR (McIntosh and Lobaugh, 2004; Krishnan et al., 2011). PLSR is a method that is particularly suitable for high-dimensional regression problems, where the number of parameters is larger than the number of samples. PLSR projects the predictor variables into a latent space, while optimizing the prediction of the outcome. PLSR is similar to principal component regression, but while principal component regression only constructs components based on the captured variance of the predictors, PLSR aims to maximize covariance between factors extracted from predictors and the outcome variable.
PLSR was run in MATLAB using the plsregress function, which uses the SIMPLS algorithm (de Jong, 1993). New PLSR models were constructed for each area and each metabolite. One component was retained in each analysis. As outcome variables, we used the residuals from a regression model where voxel GM proportion predicted either Glu/tCr or GABA+/tCr. The resulting residuals were subsequently z-scored. Predictors were the p24 or p32 FC values for each participant. To statistically assess the obtained model fit (residual sum of squares; Abdi, 2010), we performed permutation tests with 1000 permutations of the outcome measure. To test whether one region explained more variance than the other area or the unparcellated MRS ROI, we first calculated the difference in R2 between two models (e.g., p24 and p32). Then, we created a null distribution by running PLSR models with 1000 permutations for both ROIs simultaneously, using the same permuted outcome vector for both ROIs. We assessed statistical significance by comparing the true difference in R2 to the permutation distribution.
Elastic net
EN models are a combination of least-absolute-shrinkage-and-selection-operator (LASSO) regression and ridge regression. EN models perform variable selection while simultaneously shrinking regression coefficients to prevent overfitting. We use EN as a complementary approach to PLSR. PLSR fits a model based on global information extracted from the feature space and outcome. It is thus able to pick up diffuse, global effects of functional connectivity on local metabolite concentrations. EN, in contrast, penalizes some regression coefficients (here FC to target ROIs) to zero, resulting in a sparse model. EN therefore more strongly enforces spatially specific effects.
Residualized GABA+/tCr and residualized Glu/tCr were predicted using EN models fit for both functional and anatomic p24 and p32 ROIs. α, the weighing term of LASSO and ridge regression in the EN, was set a priori to 0.5. The mean squared error of the model fit was estimated using 10-fold cross-validation. λ was set to the value with minimum cross-validation error. Robust β-weights used for predicting metabolite concentrations were derived from the median of 20 EN iterations. Descriptive statistics of model fit are given by the R2 of predicted metabolite values and actual metabolite values. Analogous to the PLSR models, we assessed the model fit (R2) using a permutation test where the order of the outcome vector was randomly permuted (N = 1000). Also analogous to the PLSR models, we compared the difference in explained variance between two areas using permutation tests with 1000 permutations of the outcome vector. Given the significant association between Glu/tCr and age in our sample, we ran EN and PLSR models to predict age rather than glutamatergic metabolism and compared model performance to models where age was regressed out of Glu/tCr and the predictors. To assess how spectral quality (CRLB and line width) relates to the prediction models, we computed the residuals (observed values – predicted values) of the model predicting residualized Glu/tCr, as a measure of the quality of the prediction for each participant. These model residuals were then correlated with the line widths and CRLBs, respectively. For the prediction of residualized Glu/tCr from all ROIs (functional and anatomic), lower residuals were associated with lower CRLBs (r values > 0.560, p values < 0.0001) and line widths (r values > 0.210, p values < 0.05), meaning that better MRS data quality is associated with better prediction of metabolites from functional connectivity. As predicted values for GABA+/tCr were consistently zero, these predictions were not associated with MRS quality measures.
Differential mRNA coexpression
To explore what could explain the differential association between glutamatergic metabolism and FC from the two functional areas, we investigated gene coexpression. mRNA expression data from six donor brains were obtained from the Allen Human Brain Atlas (https://human.brain-map.org). We selected a subset of genes that is associated with glutamatergic and GABAergic transmission, receptors, transport, and metabolism (Extended Data Fig. 7-1, full list of selected genes). When multiple probes were available for the same mRNA, we selected the probe with the highest differential stability (i.e., the probe with the lowest spatial variability between donors; Hawrylycz et al., 2015; Quintana et al., 2019).
We investigated correlations between functional connectivity and gene coexpression within networks. Previous research has demonstrated that regions that are functionally connected (i.e., are part of functional, distributed networks) show similar gene expression patterns distinct from those shared within other networks (Richiardi et al., 2015; Anderson et al., 2018; Huntenburg et al., 2018). The different resting-state functional connectivity from the two areas to the rest of the brain may be explained in part by glutamatergic or GABAergic gene coexpression patterns. Therefore, we further explored whether gene coexpression between the pgACC and the CONN atlas regions is differentially correlated with functional connectivity from functional p32 and p24.
To this end, we used the gene set described above and computed mRNA expression maps according to Quintana et al. (2019) for each donor. Donors' expression maps were z-scored and winsorized (threshold: absolute z-score = 3.5). Coexpression was calculated separately for each gene and functional area. mRNA coexpression was calculated as the z-transformed Pearson correlation between the area and each of the 132 target ROIs. Subsequently, we calculated the Pearson's correlation between coexpression and mean FC for each target ROI and averaged this across network. For each network and neurotransmitter, we tested whether coexpression–FC correlations differed from zero using one-sample t tests. For the networks showing nonzero correlations of gene coexpression with FC, differences in the correlation between mRNA coexpression and functional connectivity were compared using Kolmogorov–Smirnov tests for the associated neurotransmitter (GABA or Glu) each gene. For the ventral attention network (vAt) specifically, we tested whether coupling between canonical structural and individual functional aspects was lower for Glu compared with GABA. To this end, we calculated bootstrapped mean correlations of gene coexpression associations with metabolite levels and mean correlations of FC with metabolite levels (number of iterations, 50,000). We resampled the resulting bootstrapped vectors and calculated the distance between functional and canonical associations, for both glutamate and GABA separately. Statistical significance of the difference between glutamate and GABA was assessed with a one-sample t test.
Cognitive state correlates of p32 and p24 FC–glutamate association
To explore the potential functional relevance of the differential prediction of pgACC glutamate levels, we correlated the PLSR β-maps produced by models ran on functional p24 and p32, as well as the difference between the two β-maps with association Z maps, using the Neurosynth decoder (see Materials and Methods, subsection Neurosynth decoding).
Data and materials availability
Raw mRNA expression data are available from the Allen Human Brain Atlas (http://human.brain-map.org). Code for functional imaging and MRS preprocessing as well as statistical analyses and figures will be made available on request. Maps with trained loadings are shared on neurovault (https://neurovault.org/collections/CWXNJGOY/).
Results
Functional and anatomic parcellations
For the functional parcellation, we created a group ROI based on participants' pgACC MRS masks. Participant-specific masks would introduce a bias in the FC profile and could therefore inflate individual associations with neurometabolism. Therefore, we created a composite mask. In addition, we used a recent cytoarchitectonically informed parcellation of the pgACC as a second, atlas-based ROI parcellation based on 10 postmortem human samples (Palomero-Gallagher et al., 2019). This anatomic parcellation consists of MPMs of p24ab, p24c, and p32. We parcellated the MRS ROI into two clusters of similar connectivity using hierarchical clustering (for details, see Materials and Methods, subsection Connectivity-based parcellation of the MRS ROI). Seeds were fMRI voxels within the group MRS ROI; target ROIs were the 132 CONN atlas regions.
We compared the resulting functional clusters to the anatomic parcellation of the pgACC. For this purpose, anatomic maps of p24 and p32 were restricted to MRS ROIs. The overlap between the MPMs and the functional clusters was then calculated using DCs (Arslan et al., 2018; Fig. 3C). Cluster 1 overlapped primarily with anatomic area p32 (DC = 0.750), but less with area p24 (DC = 0.322). Cluster 2 overlapped with anatomic area p24 (DC = 0.706) but not with area p32 (DC = 0.079). We thus observed good concordance between the parcellation of this region based on local, mesoscale differences and a parcellation based on whole-brain, macroscopic functional differences. Hence, functional clusters 1 and 2 are referred to as functional p32 and p24, respectively, in the remainder of this text.
Functional and anatomic parcellations. A, Functional parcellation. Dark orange, Cluster 1; light orange, cluster 2. B, Anatomical parcellation. Green, p32; blue, p24. C, Overlap of connectivity-based parcellation results with maximum probability maps from Palomero-Gallagher et al. (2019). MPMs are masked with the MRS ROI for computation of Dice overlap coefficients, cl.: cluster. D, The z-scored and de-meaned FC matrix sorted by functional clusters. Dark orange, Cluster 1; light orange, cluster 2. E, Within-cluster distance decreases with functional parcellation (median: cluster 1, 2.25; cluster 2, 2.16; unparcellated, 2.83). F, Between-cluster distance increases with functional parcellation (median: cluster 1, 4.29; cluster 2, 4.27; unparcellated, 2.83).
Areas covered by the MRS voxel show differential associations with brain networks
To characterize the functional connectivity profiles of each functional area, we calculated area-to-whole-brain connectivity and performed paired t tests on the results. Figure 4A–D depicts only those voxels that showed significantly different connectivity (p < 0.05, TFCE). Compared with functional p24, functional p32 showed stronger connectivity to areas that are part of the DMN, including the precuneus and posterior cingulate cortex, inferior parietal lobe, ventromedial/medial prefrontal cortex, temporal pole, and lateral temporal cortex (Fig. 4A,B). It was also more strongly connected to the inferior frontal gyrus. Functional p24 had relatively stronger connections to areas that are associated with the vAt, including the striatum, anterior insula (AI), anterior mid-cingulate cortex, and amygdala (Fig. 4C,D). Moreover, according to the t value distributions of the seed regions, functional p32 FC is more broadly connected to most Yeo networks (Yeo et al., 2011) compared with functional p24 FC, which is more specifically associated with the attention networks [vAt and dorsal attention network (dAt); Fig. 4F].
Characterization of FC differences between functional p32 and p24. A–D, Results of paired t test on whole-brain FC from functional areas, p < 0.05, TFCE corrected. A, B, functional p32 > functional p24 (orange). C, D, Functional p24 > functional p32 (yellow). E, Word cloud of cognitive state terms associated with brain activity in the regions that are more strongly connected to p32 relative to p24. The top 40 terms most strongly associated with this spatial pattern are displayed. F, Density plots of t values from a paired t test on whole-brain functional connectivity from functional pgACC areas, p32 > p24, split by Yeo networks and the cerebellar network. Vis, Visual network; SoM, somatomotor network; Lim, limbic network; Cer, cerebellum.
To further investigate potential associations of average FC differences between functional p32 and p24 and cognition, we used the Neurosynth framework (Yarkoni et al., 2011), which comprises neuroimaging data from 14,371 fMRI studies (release 0.7). The decoder toolkit implemented within this framework allows for decoding cognitive states from a given (activation) map (Rubin et al., 2017). Compared with p24, p32 is more strongly connected to a set of regions that, when activated, are associated with cognitive states such as theory of mind, mentalizing, self-referential cognition, and social cognition, which are cognitive states in which the DMN is thought to be heavily involved (Spreng and Grady, 2010) that require strong connections to other networks (Barrett and Simmons, 2015; Teckentrup et al., 2019; Fig. 4E).
Glutamate is better predicted by p32 FC compared with p24 FC
To test whether parcellation of the pgACC MRS ROI into p32 and p24 improved the prediction of pgACC glutamate, we used two complementary machine learning approaches that enforce different degrees of sparsity. PLSR fits a model based on global information extracted from the feature space and outcome (Zeighami et al., 2019). It is thus able to pick up diffuse global effects of functional connectivity on local metabolite concentrations. We use EN as a complementary approach to PLSR. EN, in contrast to PLSR, penalizes some regression coefficients (here: FC to target ROIs) to zero, resulting in a sparse model (Pervaiz et al., 2020). EN therefore more strongly enforces spatially specific effects (see Materials and Methods).
First, we used PLSR and EN to test whether FC from areas p24 and p32 could predict pgACC Glu/tCr (residualized for gray matter proportion) better than expected by chance (Extended Data Fig. 5-1). We found that FC from functional p32 could be predicted using EN (R2 = 0.324, p < 0.001; Fig. 5B). Although the PLSR model indicated that the FC profile of functional p32 was associated with pgACC Glu/tCr (R2 = 0.181, p = 0.119; Fig. 5A), this effect did not reach statistical significance. Nevertheless, predicted Glu/tCr values of the EN and PLSR models were highly correlated (R2 = 0.543, p < 0.001), indicating that both methods picked up similar information in the connectivity profiles, and that feature selection in EN was beneficial. Analyses using anatomic p32 replicated the prediction of Glu by functional parcellation and both models were significantly better than chance (EN: R2 = 0.394, p < 0.001; PLSR: R2 = 0.263, p = 0.023; Fig. 5C,D). In contrast to functional and anatomic p32, p24 FC was less consistently predictive of pgACC glutamatergic metabolism. p24 FC was not predictive of Glu/tCr using PLSR (p values > 0.450) or of using EN (p values > 0.999; Fig. 5, Extended Data Fig. 5-1). Overall, pgACC glutamatergic metabolism was most reliably predicted from p32 FC.
Glutamate levels can be predicted from functional connectivity profiles of p32, but not p24. A, B, Results of models predicting Glu/tCr with FC profiles derived from the functional parcellation as predictors, using PLSR (A) and EN (B). C, D, Results of models predicting Glu/tCr using FC profiles derived from the anatomical parcellation as predictors, using PLSR (C) and EN (D). Violin plots denote the distribution of R2 in permutation tests (1000 permutations of metabolite levels). Bars denote the R2 value obtained using Glu/tCr (residualized for gray matter proportion in participants' MRS voxels) as the predicted variable. Detailed statistics are presented in Extended Data Figure 5-1. Results did not change when absolute concentrations of Glu were predicted (Extended Data Figure 5-2). When the influence of age was accounted for in the EN model, the distinction in predictive ability of the two areas increased (Extended Data Figure 5-3). *p < 0.05, ***p < 0.001.
Figure 5-1
*p < 0.05, ***p < 0.001. Download Figure 5-1, DOCX file.
Figure 5-2
*p < 0.05, **p < 0.01, ***p < 0.001. Download Figure 5-2, DOCX file.
Figure 5-3
Predictions of age and Glu/tCr were done using the FC from functional areas p24 and p32. The models in which the influence of age was regressed out used the residuals of Glu/tCr after regressing out both voxel gray matter content and age. *p < 0.05, ***p < 0.001. Download Figure 5-3, DOCX file.
FC from the unparcellated MRS ROI could not predict Glu/tCr better than chance using PLSR (R2 = 0.142, p = 0.295), but EN led to results comparable to those using p32 alone (R2 = 0.384, p ≤ 0.001). To test whether p32 FC could predict glutamatergic metabolism in the pgACC better than p24 FC or FC from the unparcellated MRS ROI, we compared the variance explained by two sets of predictors (e.g., p32 and p24 FC) using permutation tests. Overall, the explained variance of p32 FC was higher than p24 FC, demonstrating that p32 FC was more strongly associated with Glu/tCr compared with p24 (EN, functional: p < 0.001; PLSR, functional: p = 0.082; EN anatomic: p < 0.001; PLSR anatomic: p = 0.017). p32 FC by itself predicted Glu/tCr better or equally well compared with the unparcellated MRS ROI (EN, functional: p = 0.995; PLSR, functional: p = 0.093; EN anatomic: p = 0.110; PLSR anatomic: p = 0.001).
Given the significant correlation between unresidualized Glu/tCr and age in our sample (r(86) = −0.256; p = 0.016; 95% CI, −0.442, −0.049), we explored a possible effect of age on the prediction of Glu/tCr (Extended Data Fig. 5-3). Age by itself could be predicted from functional p32 FC (R2 = 0.441) and p24 FC (R2 = 0.310; p values < 0.001) using EN but not using PLSR (R2 = 0.210, p = 0.066; and R2 = 0.143, p = 0.289, respectively). When the influence of age was accounted for in the EN model by residualizing for age in both FC and metabolite levels, the distinction in predictive ability of the two areas increased. Functional p32 FC explained numerically more variance in Glu/tCr (R2 = 0.423) compared with the EN model that did not account for age (R2 = 0.403; Extended Data Fig. 5-3).
pgACC FC is not predictive of pgACC GABA+
To test whether pgACC GABAergic neurometabolism could be differentially predicted from the FC of parcellated voxels, we repeated the above analyses for GABA/tCr (residualized for voxel gray matter proportion; Extended Data Fig. 6-1). As the GABA signal at 3 ppm is likely to contain a contribution of macromolecule resonances, measured concentrations of GABA are reported as GABA+. A PLSR model built on the FC profile of functional p24 numerically explained most variance (R2 = 0.185). However, none of the PLSR or EN models built on FC profiles of anatomic or functional p24 and p32 could predict GABA+/tCr better than chance. FC from the unparcellated MRS voxel was also not predictive of pgACC GABAergic neurometabolism (p values > 0.05; Fig. 6). Although PLSR models yielded numerically better prediction, none of the functionally or anatomically informed models predicted GABA+ better than chance. We performed a post hoc analysis of the pgACC Glu/GABA+ ratio, a proxy of the local excitation/inhibition balance. None of the models revealed significant associations.
GABA levels cannot be predicted from functional connectivity profiles of the pregenual anterior cingulate cortex. Results of PLSR and EN models for GABA+/tCr. A, B, Results of models predicting GABA+/tCr with FC profiles derived from the functional parcellation as predictors, using PLSR (A) and EN (B). C, D, Results of models predicting GABA+/tCr with FC profiles derived from the anatomical parcellation as predictors, using PLSR (C) and EN (D). Violin plots denote the distribution of R2 values in permutation tests. Bars denote R2 obtained using GABA+/tCr as the predicted variable (residualized for gray matter proportion in participants' MRS voxels). Detailed statistics are presented in Extended Data Figure 6-1. Results did not change when absolute concentrations of GABA+/tCr were predicted (Extended Data Figure 5-2).
Figure 6-1
*p < 0.05, **p < 0.01. Download Figure 6-1, DOCX file.
Coexpression of GABAergic and glutamatergic genes
To investigate why glutamate may be more robustly predicted from area FC than GABA, we explored whether glutamatergic and GABAergic gene coexpression patterns are differentially associated with area FC. Previous work has shown that glutamate and GABA levels across the entire cingulate cortex follow glutamate and GABA receptor fingerprints (Dou et al., 2013). Neurotransmitter receptor density fingerprints shape the local excitatory/inhibitory balance, which influences baseline resting-state functional connectivity (van den Heuvel et al., 2019). van den Heuvel et al. (2016) showed an association between the ratio of excitatory and inhibitory gene expression and cortical resting-state FC. In addition, recent work has shown that regions within functional networks share gene expression patterns (“gene coexpression”) that are distinguishable from those shared within other networks (Richiardi et al., 2015; Anderson et al., 2018; Huntenburg et al., 2018).
The differential prediction of glutamate and GABAergic metabolism may be explained in part by different correlations between resting-state FC patterns of pgACC areas, and canonical glutamatergic and GABAergic gene coexpression patterns obtained from the Allen Human Brain Atlas (see Materials and Methods, subsection Differential mRNA coexpression; Extended Data Fig. 7-1, full list of selected genes). We first assessed whether coexpression and pgACC FC were correlated in our sample, as a certain amount of coexpression of GABAergic or glutamatergic genes may be needed for individual variation in metabolism–FC correlations to be meaningful. Intriguingly, two of the networks that were most strongly associated with p32 and p24 FC—DMN and vAt—also showed nonzero correlations between gene coexpression for both GABA and glutamate and pgACC FC (Extended Data Fig. 7-2). A third network, the frontoparietal network (FPN), also showed nonzero correlations between FC and gene coexpression between pgACC and the network targets. For these three networks, a significant portion of variance in group FC is explained by canonical, structural factors such as gene expression.
As a next step, we investigated differences between GABAergic and glutamatergic coexpression–FC correlations, focusing on the networks for which an association between FC and gene expression was apparent. Coexpression of GABAergic genes within the DMN and vAt was more strongly coupled to FC within those networks relative to glutamate (DMN: distance (D) = 0.360, p < 0.001; vAt: D = 0.283, p = 0.011; Fig. 7, Extended Data Fig. 7-2). This was not the case for the FPN (D = 0.147, p = 0.488). For glutamate-associated genes, this relatively reduced coupling between canonical structural and individual functional aspects could leave room for individual variation in Glu levels to influence FC. This could in part explain why compared with GABA+, pgACC Glu levels could be better predicted from individuals' FC profiles.
Compared with GABA, there is reduced coupling of glutamate-associated gene coexpression with pgACC FC to vAt and DMN. Density plots show the distribution of glutamate-associated and GABA-associated gene coexpression correlations with p32 and p24 FC within the same functional network. Gene coexpression is calculated as the Pearson correlation between pgACC p24 and p32 mean mRNA expression with mean mRNA expression in target ROIs. Coexpression–FC correlations are then calculated for each gene and averaged across Yeo networks. The full list of genes used is presented in Extended Data Figure 7-1. For detailed statistics, see Extended Data Figure 7-2 and Extended Data Figure 7-3. *p < 0.05, ***p < 0.001. Vis, Visual network; SoM, somatomotor network; Lim, limbic network; Cer, cerebellum.
Figure 7-1
AHBA, Allan Human Brain Atlas. Download Figure 7-1, DOCX file.
Figure 7-2
Results of one-sample t tests testing whether mRNA coexpression correlations with FC are different from zero, for each network and neurometabolite. Vis, Visual network, SoM, somatomotor network, Lim, limbic network; Cer, cerebellar network. **p < 0.05, ***p < 0.001. Download Figure 7-2, DOCX file.
Figure 7-3
Two-sample Kolmogorov–Smirnov tests comparing the distributions of z-transformed GABAergic and glutamatergic mRNA coexpression and FC correlations. Vis, Visual network, SoM, somatomotor network; Lim, limbic network; Cer, cerebellar network. *p < 0.05, **p < 0.01, ***p < 0.001. Download Figure 7-3, DOCX file.
By comparing FC associations with metabolites versus gene coexpression, we can illustrate this possibility using our data. For each target ROI, the PLSR β-weight represents the strength of the relationship between glutamatergic or GABAergic metabolism and FC to an area. Analogous to gene coexpression (Fig. 7), the distribution of PLSR weights across the vAt was significantly different from zero (Fig. 8A, Extended Data Fig. 8-1). For this network, the discrepancy of bootstrapped mean correlations between FC and gene coexpression versus ones between FC and local metabolism was smaller for glutamate (Fig. 8B), such that metabolite–FC correlations were relatively stronger, and coexpression–FC correlations were relatively lower (pboot = 0.008).
Associations of FC profiles with neurometabolites are constrained by canonical gene expression. A, PLSR β values by Yeo network (pooled across functional p24 and p32) show that FC to vAt is associated with GABAergic and especially with glutamatergic metabolism (Extended Data Figure 8-1, detailed statistics). B, Canonical gene coexpression of GABAergic genes within pgACC and vAt is more strongly correlated to FC than pgACC GABA+ levels. For glutamate, this discrepancy is less strong, suggesting that pgACC glutamate levels have more room to influence FC to the vAt. Density plots in B display distributions of bootstrapped r values (50,000 iterations).
Figure 8-1
Results of one-sample t tests testing whether PLSR β-weights (pooled across p24 and p32) are different from zero, for each network and neurometabolite. *p < 0.05, **p < 0.01, ***p < 0.001. Download Figure 8-1, DOCX file.
Cognitive state correlates of p32 and p24 FC–glutamate associations
Next, we investigated the potential functional relevance of the differential prediction of pgACC glutamate levels. PLSR β-weight maps for functional p32 and p24 show different spatial patterns (Fig. 9A,B, top, Extended Data Fig. 9-1, analogous results for anatomic areas). To explore differences between the two pgACC areas, we applied the Neurosynth decoder tool to their respective β-weight maps. Regions for which FC to p32 is predictive of glutamate, when activated, correspond to cognitive states related to the ventral attention network and cognitive control (Fig. 9A, bottom). In contrast, regions important in the prediction of glutamate from p24 FC, when activated, reflect cognitive states related to movement and somatosensory experiences (Fig. 9B, bottom). The difference between the two β-maps suggests that relative to p24, p32 FC to the precuneus and posterior cingulate (i.e., DMN regions), anterior insula (vAt) as well as subcortical regions is linearly associated with glutamate levels (Fig. 9C). In general, increased glutamate in the whole pgACC appears to strengthen p32 FC across a wider functional spectrum compared with p24, whose FC coupling is restricted to the somatosensory and motor regions (Fig. 9D).
Associations of glutamate loadings with Neurosynth terms show that p32 FC extends more to other brain regions. A–C, Top, PLSR β-maps show the association between pgACC glutamate levels and pgACC area-to-target FC. A–C, Bottom, Decoded PLSR β-weight maps show the 40 cognitive state terms most strongly associated with the PLSR β-weight maps. D, Spider plot of correlations with the top five terms (excluding duplicates) associated with regional brain activity in the spatial pattern represented by PLSR β-weights of p32, p24, and p32–p24 PLSR, respectively. AI, anterior insula; RI, Response inhibition; SMA, supplementary motor; S2, secondary somatosensory; PCC, posterior cingulate. Analogous results for anatomic areas are presented in Extended Data Figure 9-1.
Figure 9-1
Associations of glutamate loadings with Neurosynth terms show that p32 FC extends more to other brain regions. A–C, Top, PLSR β-maps show the association between pgACC glutamate levels and pgACC area-to-target FC. A–C, Bottom, Decoded PLSR β-weight maps show the 40 cognitive states most strongly associated with the PLSR β-weight maps. D, Spider plot of correlations with the top five terms (excluding duplicates) associated with PLSR β-weight maps of p32, p24, and p32-p24 PLSR, respectively. SMA, Supplementary motor; RI, response inhibition; S2, secondary somatosensory; PCC, posterior cingulate. Download Figure 9-1, TIF file.
Discussion
Psychiatric disorders such as depression and anxiety are characterized by changes in neurometabolites (Pollack et al., 2008; Horn et al., 2010; Colic et al., 2019) and functional connectivity (Mulders et al., 2015). One of the key challenges in biological psychiatry is to understand how changes in neurotransmission lead to changes in brain function so that pharmacological interventions could be used to “normalize” brain function and improve behavioral symptoms (Wang and Krystal, 2014; Allen et al., 2019; van den Heuvel et al., 2019). However, technical limitations of current neuroimaging techniques such as large MRS voxels encompassing heterogeneous areas make the link of metabolites to function of brain networks nontrivial (a “many-to-one mapping problem”; Paulus and Thompson, 2019). Here, we demonstrate that parcellating an MRS region based on functional connectivity or cytoarchitecture improves the prediction of local neurometabolism via global connectivity profiles. Restating the question as a problem of classification rather than localization—can resting-state FC from an area predict local neurometabolism or not?—allowed us to reduce the number of comparisons and address the many-to-one mapping problem. Specifically, we found that area p32 FC predicts glutamate better than chance regardless of the parcellation scheme, whereas we did not find converging evidence for the prediction of GABA+ from FC profiles. Moreover, prediction of glutamate from p32 FC explained as much or more variance than FC from the unparcellated MRS ROI, while providing additional spatial information. Collectively, our results show that multimodal imaging may help to overcome the fundamental limitations of a single method, as fMRI can improve the spatial specificity of local glutamatergic metabolism assessed with conventional MRS.
Hierarchical clustering of the pgACC MRS voxel recovered clusters in line with anatomic areas p24 and p32 (Palomero-Gallagher et al., 2019), with distinct functional connections during rest. The area most predictive of pgACC glutamate, p32, is well connected to most functional networks but, compared with p24, it has particularly strong FC to regions that are part of the DMN. Area p24 showed relatively stronger connectivity to parts of the ventral attention network or salience network. These findings are in line with previous work in humans (Beckmann et al., 2009; Palomero-Gallagher et al., 2019) and nonhuman primates (Pandya et al., 1981; Vogt and Pandya, 1987; Carmichael and Price, 1995). Moreover, similar distinctions between p24 and p32 functional domains were found using direct stimulation of the cortex using stereo-electroencephalography (Caruana et al., 2018). We demonstrated that functional connectivity-based parcellation of an MRS ROI can reflect both cytoarchitectonic areas and well established connectivity differences.
How should the differential prediction of glutamate FC from the two pgACC areas be interpreted? One possibility is that there is an association between glutamate and p32 FC, because of the stronger connectivity of the area to the DMN. This intrinsic connectivity network shows the highest metabolic activity at rest (Raichle, 2001). Intriguingly, our findings suggest another possible mechanism. pgACC glutamate increases FC to networks associated with a relatively broad functional range, but particularly to regions of the vAt. Glutamate concentrations in the pgACC may therefore be an important factor contributing to the ability of the region to switch between exteroception (vAt or salience network) and interoception (DMN). The interaction between these networks is often dysregulated in psychiatric disorders (Menon, 2011; Manoliu et al., 2014; Kaiser et al., 2015; Teckentrup et al., 2019). In a sample of 22 healthy participants, the sum of glutamate and glutamine could not predict functional connectivity between pgACC and the left anterior insula, part of the vAt (Horn et al., 2010). In patients with MDD, the authors found a reduced anticorrelation of the pgACC and AI. Although our sample is healthy, it is possible that with our larger sample, we have captured effects similar to those in the study by Horn et al. (2010) and a larger spectrum of pgACC Glu/tCr concentrations. In MDD, the glutamatergic cycling between neurons and glia is thought to be reduced, supported by findings of reduced astrocyte density and markers in the anterior cingulate cortex (Bernstein et al., 2015; Rajkowska, 2003) and altered glutamine/glutamate ratios—a proxy of glutamatergic cycling—in the pgACC (Colic et al., 2019). Our approach is particularly well suited to probe how changes in the glutamine/glutamate ratio relate to functional connectivity between pgACC and connected regions in MDD. p32 in particular may be an especially relevant target for future research on metabolic and functional changes in psychiatric disorders.
We demonstrated the potential for further applications of this method by prediction of age. Age, unlike gender or voxel gray matter proportion, was highly correlated with glutamate levels. This is not surprising, as previous work has shown age-dependent decreases in glutamate levels in a variety of brain regions, including the ACC (Schubert et al., 2004; Hädel et al., 2013; but see Brandt et al., 2016), motor cortex (Kaiser et al., 2015), and striatum (Zahr et al., 2008). Changes in GABA levels with age are less consistent—a meta-analysis (Haga et al., 2009) did not find evidence for age-related changes in GABA. In addition to glutamate-related changes, a wealth of research covers FC changes that occur during aging. Particularly, intrinsic connectivity in the DMN alters with age (Damoiseaux et al., 2008; Wu et al., 2011; Tomasi and Volkow, 2012; Ferreira and Busatto, 2013). Both p24 FC and p32 FC were predictive of participants' age using EN. The successful prediction of age from both areas is encouraging and suggests that this method could be used for differential prediction of other clinical characteristics. It also demonstrates that while p24 FC is not predictive of glutamate, it is predictive of age. What is more, after regressing out the influence of age from FC and residualized Glu/tCr, the variance explained in metabolism by p24 was reduced, whereas variance explained by p32 increased compared with the model in which age was not accounted for.
In contrast to the considerable predictive power for glutamate, we did not find converging evidence for the prediction of local GABAergic metabolism via global connectivity profiles. Previous work demonstrated a positive correlation between GABA levels in the pgACC and negative BOLD responses on emotional stimulus presentation (Northoff et al., 2007), suggesting a link between GABA and the potential to downregulate DMN activity with increasing cognitive load. As our measurements were acquired at rest, it may be that an association between GABA+ and FC becomes apparent after stronger recruitment of task-positive networks. Another explanation relates to the complex association between GABAergic metabolism and the BOLD response. Depending on the brain region and network, an increase in inhibitory activity may or may not lead to increased BOLD response (Bartels et al., 2008; Logothetis, 2008). Our results showed that GABA-associated genes were more tightly linked to pgACC FC compared with glutamate-associated genes. GABA coexpression–FC relationships may thus be less variable across individuals. GABA modulates glutamatergic excitation by acting on pyramidal neurons in cortical microcircuits, and local GABA+ concentrations may therefore represent mostly local processes (Buzsáki et al., 2007; Logothetis, 2008; Isaacson and Scanziani, 2011). We explored whether the local excitation/inhibition balance (Glu/GABA+) could be predicted from FC. None of the models revealed significant associations. This is in line with previous work that showed that the excitation/inhibition balance in a DMN node was predictive of intrinsic connectivity of the DMN, but not that of other networks (Kapogiannis et al., 2013). To summarize, our results suggest that pgACC measures of GABA+ are unlikely to be associated with patterns of long-range functional connectivity at rest, calling for alternative techniques in the future.
Strikingly, while functional and anatomic parcellation performed similarly, there were vast differences in performance between EN and PLSR models. For the region most predictive of glutamate, p32, predicted values from EN and PLSR were highly correlated, suggesting they pick up similar information. Nevertheless, EN models provide a more conclusive answer on whether a metabolite can be predicted or not. In case there is no relationship between outcome and predictors for EN, the built-in 10-fold cross-validation will lead to models with all predictors reduced to zero, because there is no λ yielding a better than chance out-of-fold prediction. In such situations, PLSR models can still result in high explained variance, because these are not equipped with a way to penalize predictors that do not remain predictive in held-out folds. EN models have previously been used to predict behavior (Kashyap et al., 2019) and disease (Teipel et al., 2017) from neuroimaging data. This method outperforms multiple regression (Teipel et al., 2017; Jollans et al., 2019) and also frequently outperforms other machine learning techniques in cases where the number of participants is similar or smaller than the number of datapoints (Jollans et al., 2019). Overall, based on our results, combining MRS measures of local neurometabolites with resting-state FC might be promising to identify candidate regions or networks.
The results presented here must be considered in light of several limitations. First, the mean Dice overlap between MRS voxels and the functional areas was significantly greater for p24. Nevertheless, functional area p24 was less strongly associated with glutamatergic metabolism compared with functional area p32. In addition, there was no significant relationship between participants' MRS voxel overlap with functional clusters and their glutamate or GABA+ levels. Therefore, it is unlikely that the larger MRS voxel overlap with cluster p24 had a significant influence on our results. Second, as in most MRS studies, there are potential confounds in the quantification of neurometabolites. While GABA is more challenging to reliably quantify compared with glutamate, it is unlikely that data quality played a role in our findings. At high field strengths like 7 tesla, increased signal dispersion allows for the separation of GABA peaks from larger, overlapping resonances. Moreover, after our stringent quality control, only five participants had to be excluded from the GABA analyses. Nevertheless, it cannot be excluded that signals from macromolecules that overlap with the GABA peak at 3 ppm have contributed to our negative result for GABA+. Future work could test this by incorporating a measured macromolecular baseline in the basis set used for fitting MRS spectra. With regard to glutamate, MRS measures are not limited to glutamate as a neurotransmitter. Glutamate fulfills other, metabolic, roles in the cell, including protein synthesis and energy metabolism, which cannot be separated from a neurotransmitter (Rae, 2014). It also appears that vesicular glutamate is not detectable by MRS (Kauppinen and Williams, 1991). The generalizability of the models needs to be assessed in an independent dataset. For this purpose, β-maps in MNI space are available for download (https://neurovault.org/collections/CWXNJGOY/). Last, the gene expression dataset was obtained from an independent sample (Human Brain Atlas) consisting of six donor brains. At present, it is unclear whether gene expression and MRS measures are well aligned and whether gene coexpression generalizes to a wider population of healthy living adults. Notwithstanding, our study shows that associations of FC with canonical gene expression can provide crucial insights into neurobiology.
To summarize, we demonstrated that combining complementary information from different neuroimaging modalities (MRS and fMRI) can provide incremental spatial information on the relationship between function and neurometabolism by capitalizing on the higher resolution of fMRI. Our results show that p32 is more predictive of pgACC glutamate compared with p24 and suggest that, although p32 as a DMN node is strongly connected to most networks, pgACC glutamate concentrations are particularly associated with p32 FC to the ventral attention or salience network. Unlike glutamate, GABA+ could not be reliably predicted from pgACC FC, as canonical GABAergic coexpression may be more influential. As smaller voxel sizes reduce SNR, this approach could be used as an alternative to extract more localized information about key neurometabolites and can be particularly informative when the MRS ROI cannot be restricted to one functionally distinct area. Importantly, this method could be applied to other multimodal datasets, including EEG-fMRI or PET-fMRI to improve the spatial resolution of inferences. Crucially, our novel combination of techniques can be readily used in existing datasets to uncover more spatially specific relationships between functional connectivity underlying neurometabolism in health and disease. Thus, a broader application of interpretable machine learning methods may lead to a better understanding of neurobiological mechanisms of common psychiatric disorders.
Footnotes
The authors declare no competing financial interests.
This research was supported by the German Research Foundation (DFG; Grants SFB779/A06 and DFG Wa2673/4-1 to M.W.) the Center for Behavioral Brain Sciences (Grant CBBS NN05 to M.W.); and the University of Tübingen, Faculty of Medicine (Fortune Grant 2453-0-0 to N.B.K. and V.T. received salary support). We thank Renate Blober-Lüer and Dr. Claus Tempelmann (Department of Neurology, Otto von Guericke University) for MR data acquisition. We also thank Larissa Katz for her help with the analysis of MRS data.
- Correspondence should be addressed to Nils B. Kroemer at nils.kroemer{at}uni-tuebingen.de or Martin Walter at martin.walter{at}med.uni-jena.de
