All studies were performed following the approval of Service de la consommation et des affaires vétérinaires du canton de Vaud (Switzerland) and according to the federal guidelines of the Animal Care and approved by the local authority (EXPANIM-SCAV). Male Sprague-Dawley rats (n = 15, 350 ± 40 g; Charles River, L'Arbresle, France) under isoflurane anesthesia (2–3%) vaporized in 30% O₂ in air were intubated, and mechanically ventilated. Two femoral arteries and one femoral vein were catheterized for blood gas sampling and blood pressure measurements as well as α-chloralose (an initial intravenous dose of 80 mg/kg was administered followed by a continuous intravenous infusion of 27 mg/kg/h at a rate of 2 ml/h) and pancuronium administrations. Respiration rate was monitored through a pillow (SA Instruments, Stony Brook, NY, USA) placed underneath each rat. Temperature was measured using a rectal sensor and regulated via control of the temperature of water flowing through tubing covering the body of each rat and linked to a temperature-regulated bain-marie. Less than 300 μl of arterial blood were sampled every 30 min and blood parameters directly measured using an AVL blood gas analyzer (Dotmed, USA). Mean Arterial blood pressure (MABP) was measured continuously using a transducer attached to the femoral artery catheter. Body temperature and blood parameters were maintained at physiological levels (T = 37.5°C ± 0.5°C; pH = 7.4 ± 0.05, pCO2 = 39.7 ± 7 mmHg and MABP = 148.9 ± 11 mmHg) throughout each experiment. An intravenous femoral injection of Pancuronium Bromide (Sigma, Switzerland) of 0.7 ml per hour was performed to minimize tremors. Rats were positioned in a dedicated stereotactic holder equipped with ear and bite bars which was tilted in the magnet (30–45°) for a better positioning of voxels for fMRS over the barrel cortex. fMRS measurements were conducted sequentially in thalamus and S1BF followed by BOLD fMRI in 15 rats. The resting-state fMRI analysis is presented as supplementary data in the present study.

Electrodes were percutaneously inserted in the left infraorbital nerve. Electrical stimulation of the left trigeminal nerve (TGN) was performed using an external stimulator (WPI, Stevenage, UK) as described in Just et al. (2010). The paradigm of stimulation (1 minOFF–1 minON…) for both fMRI and fMRS was repeated for 32 min with pulse duration of 0.5 ms, stimulation frequency of 1Hz and stimulation current amplitude of 2 mA (Sonnay et al., 2015).

Experiments were described schematically in Figure 1. Experiments were performed on an actively shielded 9.4T/31 cm bore magnet (Agilent, USA) with 12 cm gradients and a surface coil. Shims were adjusted using FAST(EST)MAP (Gruetter and Tkác, 2000).

Localized proton spectroscopy was performed using Spin Echo Full Intensity Acquired Localized Sequence (SPECIAL, TE = 2.8 ms, TR = 4 s) (Mlynárik et al., 2006) in 15 rats. The voxel of interest for the thalamus was chosen by reference to the Paxinos and Watson atlas (Paxinos and Watson, 1998) so that it encompasses the VPM and POm structures (n = 12). Its size was 3 × 3 × 4 mm³ and it was shimmed down to a linewidth of 12 ± 2Hz. In the barrel cortex, the VOI size was 1.5 × 3 × 5 mm³ and was shimmed down to a linewidth of 10 ± 2Hz (n = 15). For each VOI, ¹H spectra were acquired during 32-min of rest followed by 32 min of TGN stimulation corresponding to 480 scans (30 × 16) per period. The water signal was suppressed using the VAPOR module containing a series of seven 25 ms asymmetric variable power RF pulses with optimized relaxation delays. To improve the signal localization, three modules of outer volume saturation (OVS) were interleaved with the water suppression pulses. The raw ¹HMRS spectra corrected for frequency drift and summed were used for LCModel analysis with a basis set of 21 simulated metabolites (Macromolecules: Mac; Scyllo-inositol: Scyllo; Ala: Alanine; Ascorbate: Asc; Aspartate: Asp; β-hydroxybutyrate: bHb; Glycerophosphocholine: GPC, Phosphocholine: PCho; Creatine: Cr; Phosphocreatine; γ-aminobutyric acid: GABA; Glucose: Glc; Glutamine: Gln; Glutamate: Glu; Glutathione: GSH; myo-inositol: Ins; Lactate: Lac; N-acetylaspartate: NAA; N-Acetylaspartylglutamic acid:NAAG; Phosphatidylethanolamines: PE; Taurine: Tau). Absolute metabolite concentrations were obtained using unsuppressed water signal (8 averages) as an internal reference. The Cramer–Rao lower bounds were used as a reliability measure of the metabolite concentration estimates. Higher concentration metabolites (Glu, NAA, Ins, Tau, PCr, Cr) with Cramer-Rao lower bounds (CRLB) under 10% were considered to be reliably quantified whereas for Glc, GABA, Lac, Asp CRLB under 30% were considered acceptable and were kept for further analysis.

Neurochemical profiles in thalamus and S1BF of the same rat were acquired serially for 12 rats starting with the thalamic nuclei and S1BF at rest followed by successive 32-min TGN stimulation periods for S1BF and thalamus. A 15-min period of time was left before acquiring in S1BF at rest, which was used to tilt the holder and perform further adjustments (morphological T2 imaging and Shim) as well as to return the holder to initial position between acquisitions in S1BF and thalamus. A one-way ANOVA test with Bonferroni correction was used to compare metabolite concentrations at rest and during stimulation. The significance level was set at 0.05. All the results are presented as Mean ± SEM. A two-way ANOVA test with Bonferroni correction was performed to compare Glu, GABA, Lac, Glc, and Gln between thalamus and S1BF brain structures during ON (stimulation) and OFF (rest) measurement times.

BOLD responses in the barrel cortex and thalamus were measured post-fMRS measurements using single shot gradient echo EPI in 15 rats (7 extra rats underwent BOLD fMRI only, a total of 22 rats underwent BOLD–fMRI). First and second order shims were adjusted using FAST(EST)MAP (Gruetter and Tkác, 2000) resulting in water linewidths of 12 ± 3 Hz in a 216 μl (6 × 8 × 4.5 mm³) volume. After echo re-alignment using a reference scan, BOLD responses were assessed using single shot gradient echo EPI (TR/TE = 2500–2000/25 ms; FOV = 20 × 20 mm; matrix = 64 × 64; slice thickness = 1 mm; 8 slices, Bandwidth = 325 KHz, 900 volumes).

Images were analyzed with SPM8 (Matlab; The Mathworks; Natick, USA; Statistical Parametric Mapping, www.fil.ion.ucl.ac.uk/spm/) using the general Linear model (GLM) analyzing each voxel independently and creating a parametric map of statistical significance (Friston et al., 1995). Gradient echo (GRE-EPI) time series were realigned, motion corrected, slice time corrected, normalized to each rat's anatomical images and spatially smoothed with a 3D Gaussian kernel (0.6 × 0.6 × 1 mm³). Within each analysis, the mean global intensities were mean scaled to an arbitrary value calculated within SPM8. The design model tested was a comparison between “off” and “on” conditions within each TGN stimulation paradigm. The paradigm was convolved with SPM's haemodynamic response function defined as a gamma-variate function and high pass-filtered. Residuals for the realigned rat movement were taken into account by submitting the realignment parameters (translations and rotations) as regressors. T-maps were calculated on a voxel by voxel basis. Thresholding criteria of 5 adjacent voxels, each with a T-score > 3.0 were used to identify regions of interest. Only clusters comprising at least 5 voxels were considered significant (p < 0.0001, corrected for multiple comparisons.

With STIMULATE (University of Minnesota, Minneapolis, USA) (Strupp, 1996), regions of interest (ROIs) over the activated primary somatosensory barrel field cortex (S1BF) were drawn with respect to the Paxinos and Watson's Atlas (Paxinos and Watson, 1998). ROIs were delineated from the thresholded t-map of each rat and had the same size as VOIs for 1H-MRS. A representative average time-course was recorded for each animal. When needed baseline correction was performed. T-maps were overlaid on single shot gradient echo EPI images.

Correlations were performed between BOLD responses in thalamus and S1BF, and [Glu] and [GABA] measured at rest and during stimulation in thalamus and S1BF. In addition BOLD responses were also correlated to changes in [Glu] and in [Lac] in S1BF only. Only metabolites with Cramer-Rao lower bounds (CRLB) under 30% were used for statistical analysis. To examine the association between BOLD responses and metabolites and between metabolites, the Shapiro-wilk test in SPSS 22 was performed for all the variables demonstrating p-values above 0.05 and therefore normality of each variable distribution.

The correlation analysis was performed within Origin (OriginLab version 9, Massachusetts, USA). Data were tabulated and scatter plots were obtained. Subsequently, a linear regression analysis was performed allowing to calculate the Pearson's r coefficient of correlation. Linear regression was peformed concomitantly with an ANOVA test with a threshold defined as 0.05. However, appropriate statistical analysis is required to perform multiple comparisons. With the conservative Bonferroni correction (Bretz et al., 2011) one would assume that all variables were independent which was not necessarily the case. Using a principal component analysis method within Origin, the number of effective comparisons (or eigenvalues) was determined as decribed by Cheverud et al. (1983). The adjusted threshold after Bonferroni correction was therefore p = α/4 with α = 0.05. Uncorrected p-values were denoted as P and were reported as correlating positively or negatively for P < 0.05. These values were comparable to values reported in the literature for a similar sample size (Table 3).

In order to examine the robustness of the correlation analysis performed in the present study, we further performed a Spearman- Rho correlation analysis within SPSS 22 since for small sample size data sets, a permulation analysis is automatically performed (Bretz et al., 2011). Finally, multiple linear regression analyses were performed in SPSS 22 for a sample size of 15 animals and for a sample size of 22 animals where the animal number was increased by adding 7 animals with CRLBs below 40% for Glu and GABA, since this type of regression analysis is valid for sample sizes of at least 20 subjects. Using a Kolmonorov-Smirnov test, the multivariate normality was checked again. To avoid over-fitting, a stepwise model analysis was used together with F statistics. The validity of the regression analysis was assessed and a Durbin-Watson test was included to check for autocorrelation and included a colinearity diagnostic. Homoscedasticity and normality of residuals were tested using standardized plots. To test the robustness of the multiple linear regression analysis Cook's and leverage tests were performed.

The neurochemical profiles of the thalamus and the barrel cortex before and during TGN stimulation were obtained successively for each rat. High spectra SNR levels were reproducibly measured in S1BF (91 ± 8) and in the thalamus (62 ± 6).GABA, Glutamate, Glutamine, Glucose, and Lactate levels were measured at rest and during stimulation periods both in the thalamus (n = 12) and in S1BF (n = 15) in VOI represented in Figure 2A. In order to increase SNR levels for thalamic MRS acquisitions, the thalamic voxel encompassed several thalamic nuclei of interest [VPM, Reticular Thalamic nucleus (RTN), and Posteromedian Thalamic Nucleus (POm) that can be recognized by reference to the labeled map. Examples of labeled spectra (mean spectra averaged over rat population and over time)] in both structures acquired during stimulation and rest periods are shown (Figure 2B). Increased and decreased levels of Glu and Lac levels during TGN stimulations can be visualized in each structure respectively. Figures 2C,D illustrate positive correlations between BOLD changes and the percent change in amplitude of NAA and PCr+ Cr peaks respectively. As a consequence of BOLD effect, the linewidths of these peaks should decrease while the peak height should increase. No relationship was found in the thalamus. The neurochemical profiles (± standard error of the mean) are depicted (Figure 3) for each structure and each condition. In the present work, Glutamate and GABA concentrations were compared as well as Glucose and Glutamine in 15 rats (1 rat demonstrated lipid contamination during stimulation and was discarded from further analysis). In the barrel cortex, Lac, Glu, GABA levels were significantly higher during stimulation (p = 0.0015, p = 0.015, and p = 0.021 respectively, one-way ANOVA and Bonferroni correction). In the thalamus, Gln was significantly higher during TGN stimulation (p < 0.009, one-way ANOVA and Bonferroni correction). Although not significantly due to high variability across the rat population, glucose levels decreased in S1BF and thalamus. Regional differences were demonstrated between thalamus and S1BF metabolites at rest and during stimulation using a two-way ANOVA test with Bonferroni correction. Significant differences between brain structures were found for Glu (p < 0.0001), Lac (p = 0.02), and GABA (p = 0.005) during rest and stimulation periods. Changes in [Glu] and [GABA] in both S1BF and thalamus were significantly related to changes related to activation (p = 0.02).

Although rats underwent sequential prolonged MRS measurements in S1BF and thalamus successively during rest and stimulation periods, BOLD responses were still detected in S1BF and in the thalamus respectively. Unexpectedly, a negative correlation (Pearson's r = −0.68, p = 0.0024) was found as shown in Figure 4A, between BOLD responses in S1BF and thalamus. In order to investigate further this relationship, correlations between BOLD responses, [Glu] and [GABA] were performed within each structure and between structures. BOLD changes in the barrel cortex were negatively correlated to resting GABA levels (Pearson's coefficient = −0.72, p = 0.003) (Figure 4B). Resting Glu levels correlated positively with BOLD changes (Pearson's r = 0.72; P < 0.03, Figure 4C) and changes in glutamate (ΔGlu) were negatively correlated to BOLD (Pearson's coefficient = −0.70; P = 0.024) (Figure 4D). A trend toward positive correlation between BOLD changes and Lac changes (ΔLac) was also found but not significant (Pearson's r = 0.47; P < 0.08) (Figure 4E).

In the thalamus, resting or stimulated GABA levels were not correlated to thalamic BOLD responses. Resting S1BF GABA levels demonstrated a positive interaction with thalamic BOLD (Pearson's r = 0.67 and P = 0.016, Figure 4F). Resting thalamic Glu levels were not correlated with thalamic BOLD changes (Pearson's r = 0.57; P < 0.09) and with BOLD changes in the barrel cortex (Pearson's r = −0.55, P = 0.09). A negative correlation between S1BF BOLD and [Glu] during stimulation was also suggested (Pearson's r = 0.62 and P = 0.055, Figure 4G).

Since the metabolic status of a structure may be predictive of its functional status including connectivity (Kapogiannis et al., 2013; Castro-Alamancos and Gulati, 2014), the relationships between [Glu] and [GABA] within and between structures were also investigated. As the energetic status of the whole brain is tightly controlled at rest, changes due to neuronal activity also involve changes between related structures. As reported in Table 1, during TGN stimulation, there was no correlation between [GABA] and [Glu] whereas at rest, GABA and Glu levels in the barrel cortex were positively correlated (Pearson's r = 0.74; P = 0.014, Figure 5A) and resting and stimulated [Glu] in S1BF were positively correlated (Figure 5B). On the other hand, resting thalamic GABA levels were negatively correlated to barrel cortex GABA levels measured during stimulation (Pearson's r = −0.69; P = 0.024) (Figure 5D) but not at rest (Figure 5C). However, Glu levels in S1BF at rest strongly influenced stimulated thalamic GABA levels (Figure 5D; Table 1). In S1BF, resting and activated Glu levels were correlated (Pearson's r = 0.67; P < 0.05, Figure 5E).

Table 2 presents the results of the Spearman-Rho correlation analysis. Results were in agreement with previous Pearson's correlation (p-values under 0.01 and 0.05).

Stepwise non-parameteric multiple linear regression analyis in a sample size of 15 animals revealed that resting Glutamate and GABA levels added significant explanatory power to the regression model. In addition, resting GABA and glutamate levels in the barrel cortex (BC) explained between 52.3 and 89% of the variance. The Durbin-watson test (D = 1.631) allowed to assume that there was no significant first order linear autocorrelation in the multiplelinear regression data. The F-Test is the test of significance of the multiple linear regression. The F-test of was highly significant (F = 19.063, p = 0.012 for Glu_rest and p = 0.001 for GABA_rest,) thus it can be assumed that there is a linear relationship between the variables in our model. Since we have multiple independent variables in the analysis the Beta weights compared the relative importance of each independent variable in standardized terms. We found that a regression model encompassing BOLD and resting Glutamate and GABA levels had the highest impact on the regression (β = 0.955 and 0.383). Multicollinearity in our multiple linear regression model was higher than MC = 0.811 showing that there was no suspiscion of multicollinearity. Q-Q plots indicated that in our multiple linear regression analysis there was no tendency in the error terms. Finally, Cook's and centered leverage mean values were 0.223 ± 0.4 and 0.273 ± 0.2 showing that exclusion of data points wouldn't have changed the regression statistics substantially while the centered leverage value demonstrated that there was no significant influence of specific datapoints. The same analysis was performed with a sample size of 22 animals showing that resting glutamate levels only added significant explanatory power to the regression model explaining 80% of the variance confirmed by highest Beta weights (β = 0.896). Again, there were no significant influence of autocorrelations and multicollinearity (D = 2.3, MC = 0.655). F-tests were again highly significant (F = 44.837, p < 0.0001). Finally due to the slight increase in sample size, Cook#s and centered leverage values were further decreased 0.077 ± 0.258 and 0.048 ± 0.045.

Our study confirmed the cortical BOLD responses dependence on resting GABA levels as well as changes in Glu and Lac in response to stimulation but revealed that the direction of the correlations was specific to the system under analysis. In agreement with previous results obtained in the visual human cortex (Northoff et al., 2007; Donahue et al., 2010; Muthukumaraswamy et al., 2012; Bednařík et al., 2015), positive S1BF BOLD signals were negatively correlated to baseline S1BF GABA levels. In the human visual cortex, concomitant positive, and negative correlations with Cerebral Blood Flow (CBF) and vascular space occupancy (VASO) signals respectively, were interpreted as a greater influx of blood to the visual cortex region in participants with higher GABA levels in order to fuel the greater amount of excitatory activity required for overcoming the higher resulting inhibition. In S1BF, low [Glu]_rest predicted lower BOLD responses whereas lower Δ[Glu] predicted higher BOLD percent changes during TGN stimulation. Furthermore, Δ[Lac] suggested a positive correlation with BOLD fMRI responses (Pearson's r = 0.47 P = 0.07). These findings were in line with the study by Bednařík et al. (2015) where higher Glu changes (Δ[Glu]) as well as higher Lac changes (Δ[Lac]) predicted higher BOLD percent changes in the human visual cortex, interpreted as a response to increased energy demands of neuronal activation. Although, this study confirmed the existence of relationships between BOLD responses and both GABA at rest and Glu changes, the negative relationship with the latter was unexpected and difficult to interpret: If higher ΔGlu levels not compensating for the inherently higher GABA levels in S1BF of some rats and therefore leading to lower BOLD responses could represent a reasonable interpretation, null ΔGlu (at intercept) resulting in highest BOLD responses appeared contradictory to previous results. To date, the majority of fMRS studies reported weak mean Δ[Glu] ranging between 2 and 4% due to cortical stimulation (Mangia et al., 2007; Lin et al., 2012; Schaller et al., 2014). In these studies, fMRS was conducted after carefully positioning the VOI over the core activated area in BOLD maps for which fMRI acquisitions were optimized. In the present study, BOLD responses were measured after fMRS acquisitions. The position of the VOI in S1BF relied on previous work (Just et al., 2013) possibly missing the S1BF activated core and thus shifting the amplitudes of vascular and metabolic responses. However, positive changes in amplitudes of NAA and tCr peaks as a consequence of BOLD changes were detected as in Just et al. (2013) suggesting that the choice of VOI was adequate. Moreover, averaging the voxels yielding the highest statistically significant S1BF BOLD responses also revealed a negative correlation (not shown) with Δ[Glu] thus discarding effects related to thresholding. Weak changes in [Gln] in S1BF also corresponded to rats with highest BOLD responses and could not serve as an explanation. Bednařík et al. (2015) also found a non-zero Δ[Glu] intercept for BOLD = 0%. Although they point at different sensitivities of fMRI and fMRS signals and agree that larger populations need to be examined, their study and the present one point at other effects. In particular, the effects of other neuromodulators such as norepinephrine, noradrenaline, dopamine etc…also need to be taken into account in future studies and may explain these discrepancies (Fontanez and Porter, 2006; Qiao et al., 2007; Castro-Alamancos and Gulati, 2014). Moreover, the literature reported either positive (Enzi et al., 2012) or negative correlations between BOLD and baseline Glu concentrations, which depended on the health status of patients (Falkenberg et al., 2014) and were mostly inter-regional (Duncan et al., 2014). Therefore, the direction of correlation may be specific to the network under investigation and its physiological status under rest and stimulation. Present results also showed that BOLD changes were strongly dependent on resting state metabolic conditions in S1BF.

During TGN stimulation, GABA and Glu levels were uncorrelated within each structure and between structures. On the other hand, resting [Glu] and [GABA] in S1BF were correlated suggesting a strong interplay between these neurotransmitters in S1BF to maintain the EIB at stable levels as would be expected without specific brain activity. Such a relationship was not found in the thalamus and attributed to a lack of specificity since the region of interest taken into account encompassed several thalamic nuclei in order to increase SNR levels for single voxel MR ¹H spectroscopy. Notably, the VOI emcompassed the reticular thalamic nucleus described as the only source of thalamic GABAergic inhibition. In addition, it receives inputs from the cerebral cortex but does not project to the cerebral cortex. The inclusion of this nucleus may explain differences between thalamic and cortical neurochemical profiles but also uncorrelated thalamic BOLD responses and resting thalamic GABA levels. Many studies suggested that RTN plays a significant role in thalamic gating, facilitating relevant stimuli and inhibiting others by increased inhibition of the information transfer via other relay nuclei. The shape and size of RTN make this important nucleus very difficult to avoid. Other techniques such as proton MR spectroscopic imaging (MRSI) (Seuwen et al., 2015) may be more adequate for future investigations while fMRS techniques coupled to optogenetic stimulation may help unraveling the exact contribution of this structure (Jurgens et al., 2012).

A significant negative correlation between BOLD responses measured in both S1BF and thalamus was found. In the literature, weak or absent thalamic BOLD responses with regards to forepaw or barrel cortex activations were reported in the past (Esaki et al., 2002; Keilholz et al., 2004; Zhao et al., 2008; Mishra et al., 2011; Devonshire et al., 2012). Thalamic BOLD responses are usually expected to be weaker than cortical ones due to generally higher GABA levels linked to higher densities of GABAergic neurons (Mishra et al., 2011; Devonshire et al., 2012; Tiwari et al., 2013). Here, accordingly, thalamic GABA levels were higher than in S1BF at rest and during stimulation but there was no significant increase of mean GABA levels in the thalamus upon TGN stimulation. Moreover, thalamic BOLD responses and GABA levels were uncorrelated. While these findings may be attributed to the inclusion of RTN, rats with higher resting S1BF GABA levels had high amplitude BOLD responses in the thalamus correlated to low amplitude BOLD responses in S1BF. Moreover, [Glu]_rest in S1BF were negatively correlated to [GABA]_stim in the thalamus. Therefore, rats with low Glu levels and high resting GABA levels in S1BF demonstrated lower BOLD responses in S1BF. Consequently, the release of high GABA concentrations in the thalamus could be interpreted as a way to overcome the higher excitatory activity supported by significantly increased thalamic Gln levels during TGN stimulation and thus balancing the decreased Glu levels and inducing strong BOLD thalamic activation. In the thalamus, [Gln] increased significantly during TGN stimulation while [Glu] decreased without significance. Interestingly, identical changes were reported by Xu et al. (2005) in the forepaw cortex of rats and were attributed to an augmented release of Glu and a rapid conversion into Gln through the Glu-Gln cycle between glutamatergic neurons and astrocytes. (Hirata and Castro-Alamancos, 2010) suggested that neuromodulators can control the synaptic activity of thalamocortical and corticothalamic cells and therefore control the state of their respective targets without being released there directly, thereby ensuring that both structures are under specific conditions essential for a proper connectivity. This is furthermore supported by studies suggesting the control of the synchrony of neuronal activity by regional EIB (Kapogiannis et al., 2013) and its correlation with functional connectivity. Therefore, the present study confirmed the interest in investigating relationships between neurotransmitters within the thalamus and the barrel cortex and between these two regions as potential predictors of functional connectivity, which could be of interest for studying pathologies such as absence epilepsy or schizophrenia. Duncan et al. (2014) in their review pointed at the fact that most studies looking at the roles of GABA and glutamate in brain function remain at the correlational stage while most of them include small numbers of subjects thus limiting the validation of findings. In the present work, functional connectivity was conducted on 8 rats only and showed high variability across this population (Figures S1, S2). At this stage, findings remain only indicative and should be further validated with more robust resting-state connectivity analysis. Figure 6 illustrates a preliminary schematical overview of the Glutamate and GABA modulatory effects and functional connectivity on cortico-thalamo cortical loop in a healthy rat.

In the text, correlations between either BOLD responses and metabolites or between metabolites were reported as being significant for uncorrected p-values. Upon application of corrected p-values as stated in Materials and Methods most correlations did not pass significance level which was attributed to the small number of animals included in the study. As a consequence of this, the correlations reported here point at interesting effects that remain to be validated. Nevertheless, to date, most of the correlation studies involving BOLD-metabolite relationships were not corrected for multiple comparisons. Table 3 summarizes various results from these studies with Pearson's correlation coefficients, p values and subject numbers. The column “statistics” mentions the statical test applied to the correlation. Unfortunately, these studies were mainly performed in humans. The number of subjects taken into account was in general similar to the number of rats used in the present study while significant negative correlation of brain activity measures with resting GABA were consistently found across studies including this one but not correlations with Glu levels. Although significance levels were not always adjusted, all these studies demonstrate that BOLD variability can be explained by GABA levels at least in cortical areas and the negative relationship between BOLD-resting GABA is a hallmark of neuronal activity. Here, the robustness of the correlation study was further confirmed by the Spearman-Rho correlation analysis and the multiple linear regression analysis. Although we delibaretely added Glutamate and GABA levels estimates with lower reliability (CRLB <=40%) from seven animals, we found that F statistics were significant while calculation of distance values (Cook's and leverage) increased the robustness of the analysis compared to a lower sample size of 15 animals. In addition, these analyses emphasized the roles of S1BF resting glutamate and GABA levels. The multiple linear regression analysis further demonstrated the need for higher sample sizes to validate metabolic correlations in rats but also pointed at the need for improved methods for associating MRS values as the interpretation of multiple linear regression models with an increased number of variables becomes extremely complicated.

One of the main pitfalls of the present study was its sequential design that unfortunately did not allow simultaneous assessment of BOLD, thalamic and S1BF metabolites. Using magnetic MRSI techniques such as FIDLOVs (Seuwen et al., 2015) which would at the same time contribute to preserve SNR levels and time may help while there is definitely a need for increased targeting at smaller structures with higher SNR levels for an increased specificity of studies as mentioned earlier for RTN or in pathologies (Pan et al., 2015). Finally, BOLD responses may be measured directly from fMRS acquisitions using changes in linewidths of NAA or total Creatine peaks as a result of activation which would also decrease scanning time. Again these measurements require high SNR levels per subject and require validation against standard BOLD acquisitions. In order to increase SNR levels for the temporal assessment of metabolite changes during TGN stimulation (not presented here), more fids were acquired increasing significantly the scanning time. Adaptation/habituation effects may therefore have been induced in the present study and will certainly need to be taken into account in future investigations.

In the present study, fMRS findings in the cortex of rats under α-chloralose anesthesia were similar to findings obtained in humans (Schaller et al., 2014; Bednařík et al., 2015). BOLD-GABA and BOLD-ΔLac correlations had similar directions. Nevertheless, the impact of anesthesia on the modulatory effects of Glu and GABA on thalamic and S1BF BOLD responses as well as resting-state functional connectivity cannot be neglected (Williams et al., 2010).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.