A 7T interleaved fMRS and fMRI study on visual contrast dependency in the human brain

Functional magnetic resonance spectroscopy (fMRS) is a powerful non-invasive technique that measures dynamic changes in neurometabolites. Unlike fMRI, which relies on hemodynamic coupling, fMRS provides a more proximal measure of neural activity. fMRS allows time-resolved measurements of metabolite concentration while the subject performs a mental task or responds to a stimulus. This makes it possible to examine the temporal behavior of metabolite levels involved in neuroenergetics and neurotransmission in human subjects.

Previous studies have demonstrated that stimulus-induced increases in neural activity result in concentration changes in certain metabolites (for reviews, see Koush et al., 2022; Mullins, 2018; Stanley and Raz, 2018). For example, visual stimulation resulted in increases in lactate (Lac) and glutamate (Glu), but also in decreases in glucose (Glc) and aspartate (Asp) in the primary visual cortex (V1) (e.g. Mangia et al., 2007; Prichard et al., 1991; Schaller et al., 2013). Moreover, motor tasks, painful stimuli, and even higher-order cognitive paradigms have also been found to elicit changes in cortical Glu and γ-aminobutyric acid (GABA) levels (e.g. Floyer-Lea et al., 2006; Koush et al., 2019; Mullins et al., 2005; Volovyk & Tal, 2020). In the initial development of fMRS task paradigms, prolonged sensory stimuli lasting several minutes were typically employed. This was partially related to the need for multiple spectral averages to ensure sufficient signal-to-noise ratio (SNR) for the stimulation and rest blocks. Moreover, changes in neurometabolite levels were hypothesized to be the result of increased oxidative metabolism to a new steady-state level (Mangia et al., 2007). This supported the use of sustained stimulation as considerable time was needed for these changes to become discernible. Yet, with the development of higher magnetic fields, increased SNR allowed for better spectral resolution and improved temporal resolution (Ladd et al., 2018; Mlynárik et al., 2008; Tkác et al., 2009). As such, more recent studies have also adopted shorter block paradigms (e.g. Ip et al., 2017, 2019), more in line with fMRI stimulus paradigms, and event-related experiments (e.g. Apšvalka et al., 2015; Koolschijn et al., 2023; Stanley and Raz, 2018). This approach has the advantage of reducing the habituation of the neural response, as well as allowing for a more detailed investigation of the temporal dynamics of fMRS.

Despite these developments, the exact temporal response profile for fMRS remains to be determined, highlighting the importance of further research in this area. For example, different neurometabolites may exhibit distinct temporal behavior, and different stimuli may elicit a unique temporal response. Additionally, the relation between neurometabolite response and stimulus intensity has yet to be fully characterized. It is hypothesized that, similarly to fMRI (Boynton et al., 1996; Horner & Andrews, 2009), the fMRS response will be linearly correlated with neural activity over a certain range. However, Ip et al. (2019) showed a linear relationship with image contrast for Glu and for the BOLD signal changes, using a 60 s visual stimulation paradigm. Characterizing the temporal and intensity-dependent response profiles of fMRS is crucial for a better understanding of the neural correlates of the changes observed in neurometabolite levels in fMRS.

To characterize the contrast dependency and the role of experiment timings in fMRS, we conducted an interleaved fMRS and fMRI experiment using a full-field flashing checkerboard visual stimulus. The experiment consisted of two stimulation blocks (STIM), with visual contrast levels of 10% and 100%, respectively, that were alternated with REST blocks. To limit habituation to the stimulus, each long STIM block was subdivided into shorter ON-OFF periods (30 s ON, 20 s OFF). This approach enabled us to analyze the temporal response function of neurometabolites by examining two time scales: the longer full STIM block (ON and OFF period included) and the shorter ON STIM blocks only. We employed a recently developed dynamic fitting approach for the analysis of our fMRS data (Clarke et al., 2023) and evaluated the contrast dependence for both fMRS and fMRI data.

Subjects were scanned in two sites: at the Leiden University Medical Center (LUMC) and the Spinoza Centre for Neuroimaging. A total of 20 datasets were obtained from 19 healthy subjects (4 men and 15 women aged 21 to 43; one subject was scanned at both locations). The study adhered to the guidelines from the Institutional Review Board of the AMC and the LUMC (the Netherlands). Written informed consent was obtained from all subjects prior to the study. All were healthy normal participants who had normal or corrected-to-normal vision and had no 7T MRI contraindications. No additional screening for neurological or psychiatric illness was performed.

Data were acquired on two similar 7T MR systems (Achieva, Philips, Best, the Netherlands). A head coil consisting of a quadrature birdcage transmit and 32-channel phased-array receive coils (Nova Medical, Inc., Wilmington, MA, USA) was used for all measurements. A T1-weighted (T1w) structural scan was acquired for each participant (TR/TE = 5/2 ms; FOV(AP,FH,RL) = 246 × 246 × 180 mm³; voxel size = 0.85 × 0.85 × 1 mm³; flip angle = 7°). fMRS and fMRI were interleaved with a combined dynamic scan time of 5 s (Henningsson et al., 2015). fMRS data were acquired using a sLASER sequence with FOCI refocusing pulses (Arteaga de Castro et al., 2013) and VAPOR water suppression (TR/TE = 3600/36 ms; bandwidth = 3 kHz; 1024 data-points; volume-of-interest (VOI) = 14 × 31 × 14 mm³). Two unsuppressed water spectra were obtained for eddy current correction. The VOI was placed in the primary visual cortex (V1), based on a short checkerboard localizer fMRI sequence and identification of the calcarine sulcus on the T1w image (Fig. 1B). Two outer volume suppression bands were carefully placed to suppress signals from outside the VOI, particularly the extracerebral fat signal. fMRI data were acquired using a 3D-EPI sequence (TR/TE/FA = 31/20 ms/10°; FOV(AP,FH,RL) = 240 × 136 × 240 mm³; voxel size = 1.875 × 1.875 × 2 mm³). We used dynamically alternating linear shims and a shared set of static second order terms to optimize shim settings for both the fMRS and fMRI acquisition (Boer et al., 2020). An in-house developed dielectric pad or metasurface was placed directly below the occipital cortex to maximize the transmit magnetic field (B₁⁺) homogeneity and efficiency in the occipital regions as previously described (Webb et al., 2022).

The task consisted of two 4 min STIM blocks, each with different black/white contrast: 10% and 100%. The paradigm started and ended with a 2.5 min REST block and the STIM blocks were separated by a 5 min REST block. Each STIM block was subdivided into sub-blocks of 30 s ON and 20 s OFF. The total task duration was 18.5 min with a total of 224 individual acquisitions for both fMRS and fMRI. The stimulus consisted of a full-field radial checkerboard flickering at 8 Hz at either 10% or 100% visual contrast. A gray static image of average luminosity equal to the checkerboard with a central white fixation dot was presented during rest periods. Stimuli were generated using PsychoPy v3 (Peirce et al., 2019), viewed through a mirror mounted on the RF head coil, and displayed on a 32" BOLD screen from Cambridge Research Systems (Amsterdam) or back-projected onto a screen (Leiden). Subjects were instructed to focus on the fixation point during the entire experiment. For a schematic representation of the stimulus paradigm, see Figure 1A.

Processing steps were performed with in-house Matlab scripts (R2019b, The MathWorks, Inc., USA). Phase correction and amplitude weighting parameters were calculated for each of the coils based on a water reference scan acquired at the start of the time series. This was then applied to all transients of each coil before coil combination. Subsequently, data were corrected for eddy currents, and spectral registration was used to correct for frequency and phase drifts across individual acquisitions (Near et al., 2015). A basis set was simulated based on our MRS sequence in FSL-MRS (Clarke et al., 2021). The basis set was simulated using full pulse descriptions and 60 spatially resolved points in each dimension (Landheer & Juchem, 2019; Landheer et al., 2021). The data were exported from Matlab and converted to NIfTI-MRS format (Clarke et al., 2022). A second alignment step was performed using FSL-MRS's dynamic alignment tool. Briefly, this performs a full linear combination modeling fit of each time point, using the above basis set, with all parameters except global phase and frequency shift fixed across time. Subsequently, the fitted phase and frequency shifts are removed from the data.

For each subject, all individual transients were averaged and the average spectrum was fitted using FSL-MRS (version 2.1.6) to determine spectral quality. Briefly, basis spectra are fitted to the complex-valued spectrum in the frequency domain, using our simulated basis set, that included 19 metabolites and a measured macromolecular baseline (which is publicly available: https://github.com/mrshub/mm-consensus-data-collection; Tkáč et al., 2009). The basis spectra are shifted and broadened with a Voigt lineshape model with parameters fitted to the data (one Gaussian and one Lorentzian line broadening parameter). All metabolite basis spectra are broadened with the same lineshape model and parameters except the macromolecules which had independent parameters. That is, in the fixed lineshape model, only four linewidth parameters were used for the model: a Lorentzian and Gaussian parameter for all metabolites, and a Lorentzian and Gaussian parameter for the MM basis spectrum. In the variable lineshape model, the same temporal model and parameters were used for all metabolites, and a separate set of parameters for the MM. In all cases, the Lorentzian parameter was constrained to be positive, as was the constant term of the Gaussian broadening. A complex second-order polynomial baseline is also concurrently fitted. Model fitting was performed using the truncated Newton algorithm as implemented in scipy (Virtanen et al., 2020). Metabolites were excluded from further analysis if the CRLB exceeded 20% in more than 50% of the subjects based on the average spectra. Additionally, we report the linewidth (i.e. full-width-half-maximum (FWHM)) of N-acetylaspartate (NAA) and the Cramer-Rao Lower Bounds (CRLB) of Glu as metrics for spectral quality.

All fMRS data were analyzed using dynamic fitting as implemented in FSL-MRS, which allows fitting of a dynamic signal model to all transients simultaneously (also referred to as spectral-temporal fitting or 2D fitting (Tal, 2023)). To model the temporal response of each metabolite, we use a general linear model (GLM). The GLM consists of a linear combination of regressors forming a design matrix (Fig. 1C, D). The GLM model is combined with linear combination modeling of the spectral response at each timepoint, with the model fitted to the whole data simultaneously in a least-squares approach.

The dynamic spectral fitting model of FSL-MRS was used to assess the time-dependence of the metabolite concentration on the design matrix, using our simulated FSL-MRS basis set. In the dynamic model, the phase, shift, baseline, linewidth, and concentration can be time-dependent variables. Here, we chose to fix the phase, shift, and baseline parameters across all timepoints, whereas the concentration and linewidth were constrained by the temporal model (design matrix) (Supplementary Fig. 1). We chose to vary the linewidth as well, as previous studies have demonstrated that stimulus-induced increases in the BOLD effect can induce line-narrowing in MR spectra as a result of changes in T2* (Mangia et al., 2006; Zhu & Chen, 2001). To assess the effect of varying the linewidth parameter, we evaluated two models: a fixed linewidth model and a variable linewidth model, which additionally applies the same design matrix to the Gaussian line-broadening parameter (sigma); the Lorentzian parameter (gamma) was fixed across time.

Design matrices (created using nilearn (Abraham et al., 2014)) were created, with both the 10% and 100% as separate explanatory variables (EV) in the design matrix, as well as a linear drift, a quadratic drift, and a baseline EV (Fig. 1C, D). We created separate temporal regressors for the concentration and the linewidth parameters. The linewidth was expected to reflect the hemodynamic response, so the STIM blocks were modeled following the ON-OFF blocks, convolved with the “glover” hemodynamic response function. For the concentration, the STIM conditions were modeled in two different ways: either considering the full STIM block as active condition (full-block model, including both ON and OFF parts of STIM blocks) or only modeling the sub-blocks within the STIM block as active condition (sub-block model, only ON parts of STIM blocks) (Fig. 1C, D). As we do not know the metabolite response function, we used no convolution for the concentration regressors. Spectral fitting parameters are as described above in ⁷ The fMRS acquisition and analysis parameters are reported in the MRSinMRS checklist (Lin et al., 2021; Supplementary Table 1). In order to display the individual data points of the fits, we extracted the results of the initial (single transient) fits of Glu for each subject. These results were then averaged across all subjects and smoothed using a moving average with a bin size of 32.

FMRI data processing was carried out using FEAT Version 6.00, part of FSL (FMRIB’s Software Library, www.fmrib.ox.ac.uk/fsl). Preprocessing included removal of the first two time points, motion correction, high-pass temporal filtering (0.002 Hz), and spatial smoothing using a Gaussian kernel of 5 mm FWHM. fMRI data were co-registered to the T1w scan using boundary-based registration, and then subsequently registered to MNI space. Similar to the fMRS analyses, a design matrix (“the sub-block model”) was created with both the 10% and 100% as separate EVs (but without the drift EVs), and fed into first-level analyses in FEAT.

Second-level GLMs (across subjects, random effect) were used to calculate group-level statistics using FSL’s flameo (Woolrich et al., 2004) for both the fMRS and fMRI data.

For the fMRS data, we evaluated the following statistical contrasts: a) 10% vs rest, b) 100% vs rest, c) mean activation vs rest (i.e. average over 10% and 100% vs rest), and d) 100% vs 10%. We calculated these contrasts for all metabolites that could be fitted with sufficient certainty. For some metabolites, only the combined neurometabolite level was considered, due to the high correlation between components: total NAA (tNAA, NAA+NAAG), total creatine (tCr, creatine+phosphocreatine), total choline (tCho, phosphocholine+glycerophosphocholine), and glucose+taurine (here referred to as Glc+). We ran these analyses both for the full-block model and for the sub-block model, and each using the fixed linewidth model and variable linewidth model. For visualization purposes, the data were normalized using the constant term in the design matrix (Fig. 1C, D last column). Please note that the linear drift terms are defined to be 0 at the central time point.

For the fMRI data, we evaluated the following statistical contrasts: a) 10% vs rest, b) 100% vs rest, and c) 100% vs 10% for the change in BOLD signal (based on the sub-block design). Additionally, we extracted the BOLD time courses and BOLD contrast estimate parameters within the MRS VOI location using fslmeants and featquery, respectively.

Finally, we computed Pearson’s correlation coefficients between the interindividual parameter estimates of the 10% contrast, the 100% contrast, and the 100% > 10% contrast for the Glu (fMRS), sigma (fMRS), and the BOLD (fMRI) response to assess the association between these modalities. As we compared two models, that is, full-block and sub-block, we corrected the p-value for multiple comparisons using the Bonferroni correction, resulting in a corrected α of 0.025.

Changes in concentrations of alanine (Ala), gamma-aminobutyric acid (GABA), glutamine (Gln), and lactate (Lac) were not statistically evaluated as their CRLB exceeded 20% in more than 50% of the subjects. A full list of the reported metabolites can be found in Tables 1 and 2. The mean FWHM of NAA was 14.11 (SD = 1.90), with a range of 12.28 to 19.04. Additionally, the mean CRLB of Glu in the average spectra was 4.20 % (SD = 1.04 %), ranging from 3.07% to 7.46%.

Using the variable linewidth model, we found significant line narrowing for both the 10% and 100% contrasts (10%: z = -3.33, p < 0.001; 100%: z = -4.41, p < 0.001), but no difference between the contrasts (z = -1.32, p = 0.09). We here report the results from the variable linewidth model, whereas the results of the fixed linewidth model can be found in Supplementary Table 2. We found a significant effect of visual stimulation on Glu levels for both the 10% and 100% contrast (10%: z = 2.48, p = 0.006; 100%: z = 3.21, p < 0.001), as well as across both contrasts (z = 3.21, p < 0.001). We found no significant difference in Glu levels between the 10% and 100% contrasts (z = 1.04, p = 0.15). Decreased neurometabolite concentrations during visual stimulation were observed for Glc+ during the 100% contrast (Glc+: z = -2.22, p = 0.01) and across STIM blocks for aspartate (Asp) and Glc+ (mean activation Asp: z = -2.03, p = 0.02; Glc+: z = -2.68, p = 0.004), and tNAA during the 10% contrast (z = -2.53, p = 0.005). These results are displayed in Figure 2. All statistical results can be found in Table 1.

When modeling the ON-OFF blocks during the STIM block, we found significant line narrowing for both the 10% and 100% contrasts (10%: z = -3.65, p < 0.001; 100%: z = -4.75, p < 0.001), and significantly more narrowing during the 100% compared to the 10% contrast (z = -2.15, p = 0.02). We here report the results from the variable linewidth model, with the fixed linewidth model results in Supplementary Table 3. Visual stimulation significantly increased Glu levels for the 100% contrast (100%: z = 2.65, p = 0.004), and across both contrasts (z = 2.47, p = 0.006). No effect of contrast was found on Glu levels (z = 0.20, p = 0.42). Additionally, decreased tNAA levels for the 10% contrast (z = -2.6, p = 0.004) were observed. These results are displayed in Figure 3. All statistical results can be found in Table 2.

Both the 10% and 100% contrast induced significant BOLD signal changes within the MRS voxel (10%: t = 9.17; p < 0.01; 100%: t = 14.50; p < 0.01) (Fig. 4C), and we observed a higher increase in BOLD in the 100% contrast compared to the 10% contrast (100% > 10%: t = 8.70, p < 0.01). This can also be observed in the whole brain analyses (Fig. 4A).

The temporal signal envelope of the fMRS and fMRI data is visualized in Figure 4B and D. We plotted the contrast dependency of both the Glu and fMRI response in Figure 5 and show a non-linear relation with visual contrast. We did not find any significant correlations between the interindividual parameter estimates of the BOLD signal (fMRI) and Glu levels (fMRS), nor of the BOLD signal (fMRI) and the linewidth (fMRS). This outcome was consistent for both the full-block and sub-block analyses (Supplementary Fig. 2; Supplementary Table 4).

We observed significantly higher BOLD activation for the 100% contrast compared to the 10% contrast (Fig. 5C). This finding aligns with literature showing that the fMRI response, which is related to neural activity via a linear transform, exhibits a non-linear response to visual contrast (Boynton et al., 1996; Logothetis et al., 2001). This is likely due to the non-linearity of neural activity in V1, as suggested previously (Boynton et al., 1996). Interestingly, we observe a similar non-linear relationship between visual contrast and the neurometabolic response of glutamate, as illustrated in Figure 5A and B.

Our results partially align with one previous study investigating the neurometabolite response to visual contrast (Ip et al., 2019), in the sense that we also observe a monotonic relation between visual contrast and the fMRS response. However, while we find significant Glu increases at 10% visual contrast, they report no significant Glu increases for contrasts lower than 100% (i.e. 3%, 12.5%, and 50%). This could potentially be explained by differences in terms of stimuli timings and analysis methods. Nevertheless, we were also not able to discern significant differences between the 10% vs 100% contrast, suggesting that Glu changes in V1 in response to differences in visual contrast are subtle and challenging to detect, thus requiring a larger sample.

The neurometabolic demands of different visual stimuli have been investigated in previous studies. Bednařik et al. (2018) presented chromatic and achromatic stimuli that activate two separate clusters of neuronal populations, known as "blobs" and "interblobs," that show comparable BOLD responses, but have different capacities for glucose oxidation. Despite their hypothesis, they did not observe any difference in neurometabolic responses elicited by the different stimuli. Another study investigated isoluminant chromatic stimuli at different flickering frequencies that are perceived and unperceived, but elicit comparable BOLD responses (diNuzzo et al., 2022). They observed markedly different neurometabolic responses, such that Lac and Glu increased only when the flickering was perceived, but not when unperceived. Our study adds to the evidence that changes in neurometabolism are particularly sensitive to the type of information processing taking place in V1, rather than specifically targeted neuronal subpopulations. Therefore, fMRS can provide complementary information about cortical processing, in addition to BOLD fMRI.

Our study revealed that visual stimulation elicited changes in multiple other neurometabolites, including Asp and Glc+, which have been previously reported in the context of neural activation (see Koush et al., 2022). Consistent with prior research, we observed a decrease in these metabolites following visual stimulation, although the effect was more consistent across visual contrasts for the full-block analysis. We found no significant differences in the change in Asp and Glc levels between the 10% and 100% contrasts. Reduced Glc levels are typically interpreted as the result of increased CMR_Glc in response to increased energy demands, while the opposite changes in Asp and Glu are thought to reflect an increase in the malate–aspartate shuttle (Mangia et al., 2007; McKenna et al., 2006). Furthermore, our study identified changes in neurometabolites that are less often reported, such as alterations in tNAA and tCr. A previous study using functional diffusion-weighted MRS ascribed changes in tCr and tNAA diffusion properties to modifications in neuronal microstructure during neural activation, as well as to an increase in the energy-dependent cytoplasmic streaming associated with enhanced metabolism during visual stimulation (Branzoli et al., 2013; Sappey Marinier et al., 1992). Alternatively, given that changes in these metabolites tend to be observed mostly without correction for BOLD-induced linewidth changes (Mangia et al., 2006), these effects may be due to the incomplete accounting for T2* effects by the variable linewidth model.

The characteristics of the BOLD fMRI response to visual stimulation have been extensively studied and are well documented (Jorge et al., 2018). Visual fMRI experiments commonly employ block durations of tens of seconds, eliciting a BOLD response following the hemodynamic response function. Indeed, our fMRI results showed a strong BOLD response that closely followed the experimental paradigm. However, in the case of fMRS, the existence of an analogous “neurometabolic response function” remains uncertain. To further investigate this temporal response, we examined whether a sustained response was present throughout the entire stimulus block (full-block analysis) or if there was a more rapid response akin to fMRI (sub-block analysis).

Our results indicate that the full-block model may capture the data slightly better as compared to the sub-block model, as evidenced by stronger effects for the neurometabolites (z-scores) and better fits of the model (lower AICs). This suggests that the paradigm may have elicited mostly sustained changes in neurometabolite levels. As such, modeling the OFF blocks in between the ON blocks may have increased the residuals of predicted-to-measured neurometabolite changes. This is supported by results from an animal fMRS study employing a comparable design (Ligneul et al., 2021). They also used stimulus blocks with an ON-OFF design, but with a longer OFF period compared to our study, and while they found evidence for sustained increase in Glu during the stimulation blocks, they did not detect faster neurometabolic responses following the ON-OFF design.

We observed several differences between the full-block and sub-block model, particularly for Asp, Glc+, and tCr. We speculate that these differences between the models for Asp and Glc+ may be due to low SNR of these peaks, and that dynamic fitting of these metabolites is sensitive to noise. For tCr (and other large singlet peaks), the concentration and linewidth terms can be correlated, meaning that some variance may be erroneously assigned, although this would be dependent on the extent that the concentration and linewidth terms follow the same temporal model. Alternatively, it may be that either the full-block or sub-block model is actually a better fit for fluctuations of these metabolites. Indeed, it is plausible that different neurometabolites may exhibit dissimilar temporal response profiles. Certain neurometabolites may demonstrate fluctuations at faster temporal scales, while others may manifest changes only following sustained stimulation. This would be dependent upon the respective neurometabolite's function in neuronal energetics and neurotransmission. Insights from mathematical modeling have provided evidence for both rapid as well as prolonged changes in neurometabolites. For example, sustained neurometabolite changes have been attributed to increased activity in metabolic pathways, that is, via neurometabolite synthesis (Mangia et al., 2009) whereas recent work has suggested that compartmental shifts of metabolites could explain more rapid changes in neurometabolites during fMRS experiments (Lea-Carnall et al., 2023). Future studies should focus on experiments teasing apart the contribution of different cellular compartments to the overall metabolite signal, to characterize the “neurometabolite response function.” In addition, future advancements in dynamic fitting models could conceivably integrate varied models to account for the idiosyncratic temporal profiles of different neurometabolites and hemodynamic responses.

The data in this study were analyzed using a novel dynamic fitting algorithm implemented in FSL-MRS (Clarke et al., 2023). While spectral-temporal or GLM analysis has been employed previously (Ligneul et al., 2021; Tal, 2023), most previous studies have utilized block averaging to demonstrate neurometabolic changes in fMRS experiments. Previous work suggests that compared to block averaging, dynamic fitting approaches may potentially increase the sensitivity (precision) in detecting stimulus-induced neurometabolite level changes, by reducing the number of model terms and increasing the effective SNR by analyzing all transients simultaneously. As such, it permits a better estimation of highly correlated model terms (Clarke et al., 2023; Tal, 2023). Nonetheless, the accuracy of the model and design matrix that we chose for the GLM is crucial for the accuracy and validity of the results obtained from these analyses. Short-block and event-related designs are likely to be more sensitive to the exact timing parameters of the model compared to long-block designs, as is also known from fMRI experimental design. Importantly, as mentioned above, we currently do not know the properties of the neurometabolite response function(s). Future studies are needed to experimentally validate the precise specification of the statistical model, which includes the selection and definition of design variables, the estimation of model parameters, and the assumptions underlying the model.

In light of previous investigations illustrating that local magnetic field alterations provoked by neural activity can induce linewidth alterations (Mangia et al., 2006; Zhu & Chen, 2001), we assessed the effect of incorporating a variable (Gaussian) linewidth parameter into our dynamic model. Our analysis revealed that permitting the linewidth to fluctuate significantly contributed to the model for both the full-block and sub-block approaches. Comparatively, the outcomes stemming from the fixed linewidth model (Supplementary Materials) demonstrated more significant results for the neurometabolites than those associated with a variable linewidth model. For example, significant differences in the tNAA and tCr resonances were found to be partially accounted for by adding the dynamic linewidth parameter. This corroborates block-averaging studies (e.g. Mangia et al., 2007), which underscore the necessity of correcting for linewidth to prevent false positives.

One limitation of this work is that we consistently presented the 10% contrast prior to the 100% contrast. This was chosen in order to prevent potential spillover effects from the 100% contrast. Another limitation is that, based on the quality of the spectra, we could not interpret the results of Ala, GABA, Gln, and Lac statistically, even though these metabolites are of interest due to their role in neurometabolism. In contrast to our expectations and previous findings (e.g. Bednarik et al., 2015), we did not find a negative association between the BOLD effects as estimated by fMRI and linewidth changes in fMRS. This may be due to the large interindividual variation in the sigma estimates compared to the BOLD estimates.

In conclusion, our study provides evidence of a non-linear relationship between visual contrast and both the BOLD response and the glutamate response in V1. In line with previous literature, we also identified concomitant changes in several other neurometabolites, including Asp, Glc, Tau, tNAA, and tCr. Moreover, our dynamic fitting approach allowed us to compare both sustained and more rapid neurometabolite responses. Although the data are suggestive of an overall more sustained response, future studies should explore the temporal response profiles of different neurometabolites and further refine the statistical models used for fMRS analysis. Notwithstanding its limitations, our study provides valuable insights into the complex relationship between visual contrast, neural activity, and neurometabolism, highlighting fMRS as a complementary technique to BOLD fMRI.

Data available on request due to privacy/ethical restrictions.

A.S.: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing—original draft, and Writing—review & editing. C.N.: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing—original draft, and Writing—review & editing. C.J.: Methodology and Writing—review & editing. W.Z.: Investigation, Methodology, and Writing—review & editing. S.J.: Formal analysis, Investigation, Methodology, and Writing—review & editing. W.T.C.: Formal analysis, Investigation, Methodology, and Writing—review & editing. I.R.: Conceptualization, Methodology, Writing—original draft, and Writing—review & editing.

A.S. is supported by a Dutch Research Council Veni grant (016.196.153). W.T.C. is funded by the Wellcome Trust [225924/Z/22/Z]. S.J. is funded by the Wellcome Trust [221933/Z/20/Z and 215573/Z/19/Z]. The Wellcome Centre for Integrative Neuroimaging is supported by core funding from the Wellcome Trust [203139/Z/16/Z and 203139/A/16/Z]. For the purpose of Open Access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

The authors declare no competing interests.

The authors thank Dr Vincent Boer for his help with the data acquisition software.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00031.