GABA, glutamatergic dynamics and BOLD contrast assessed concurrently using functional MRS during a cognitive task

A recurring issue in functional neuroimaging is how to link task‐driven haemodynamic blood oxygen level dependent functional MRI (BOLD‐fMRI) responses to underlying neurochemistry at the synaptic level. Glutamate and γ‐aminobutyric acid (GABA), the major excitatory and inhibitory neurotransmitters respectively, are typically measured with MRS sequences separately from fMRI, in the absence of a task. The present study aims to resolve this disconnect, developing acquisition and processing techniques to simultaneously assess GABA, glutamate and glutamine (Glx) and BOLD in relation to a cognitive task, at 3 T. Healthy subjects (N = 81) performed a cognitive task (Eriksen flanker), which was presented visually in a task‐OFF, task‐ON block design, with individual event onset timing jittered with respect to the MRS readout. fMRS data were acquired from the medial anterior cingulate cortex during task performance, using an adapted MEGA‐PRESS implementation incorporating unsuppressed water‐reference signals at a regular interval. These allowed for continuous assessment of BOLD activation, through T2*‐related changes in water linewidth. BOLD‐fMRI data were additionally acquired. A novel linear model was used to extract modelled metabolite spectra associated with discrete functional stimuli, building on well established processing and quantification tools. Behavioural outcomes from the flanker task, and activation patterns from the BOLD‐fMRI sequence, were as expected from the literature. BOLD response assessed through fMRS showed a significant correlation with fMRI, specific to the fMRS‐targeted region of interest; fMRS‐assessed BOLD additionally correlated with lengthening of response time in the incongruent flanker condition. While no significant task‐related changes were observed for GABA+, a significant increase in measured Glx levels (~8.8%) was found between task‐OFF and task‐ON periods. These findings verify the efficacy of our protocol and analysis pipelines for the simultaneous assessment of metabolite dynamics and BOLD. As well as establishing a robust basis for further work using these techniques, we also identify a number of clear directions for further refinement in future studies.

found between task-OFF and task-ON periods.These findings verify the efficacy of our protocol and analysis pipelines for the simultaneous assessment of metabolite dynamics and BOLD.As well as establishing a robust basis for further work using these techniques, we also identify a number of clear directions for further refinement in future studies.

| INTRODUCTION
Blood oxygen level dependent (BOLD) functional MRI (fMRI) techniques have been widely adopted for studying patterns of activation in neural systems and networks of the healthy and the diseased brain (see References 1 and 2 for overviews).While informative, a more comprehensive understanding of complex mental phenomena requires consideration across multiple complementary 'levels of explanation'. 3Aligning neuroimaging findings with neurotransmitter changes at cellular level remains a key goal.
MRS techniques allow non-invasive, quantitative assessment of neurotransmitter concentrations, including measurement of glutamate (Glu) and γ-aminobutyric acid (GABA), the primary excitatory and inhibitory neurotransmitters respectively. 4These data would typically be collected separately from (functional) imaging sequences, usually without task stimulus.This absence of direct concurrency makes it difficult to draw firm conclusions regarding the instantaneous, dynamic interrelation of the measured signals: neurochemistry and BOLD markers.Thus, developing new approaches for simultaneous acquisition of MRS and BOLD is a priority challenge when it comes to relating brain activity to mental events in the sense of cognitive performance.
There are several challenges associated with MRS techniques, particularly when considered in relation to task performance.First, the metabolite signals of interest are several orders of magnitude weaker than the water signal, which forms the basis of BOLD-fMRI measurements.For reliable measurement, the MRS signal is usually acquired from a large volume, over a large number of repetitions.This is in direct competition with demands for good anatomical specificity and high temporal resolution, as necessary for observing subtle, short-term functional dynamics localized to a particular brain region or network.
The measurement of GABA presents further challenges: significant spectral overlap with other neurotransmitters makes reliable separation difficult at commonly available field strengths, 5 with spectral editing techniques 6,7 such as MEGA-PRESS often used to address this.Such techniques places additional demands on the amount of signal that must be acquired to achieve an acceptable signal-to-noise ratio (SNR), due to the need for two distinct yet well matched sub-spectra and losses associated with imperfect editing, compounding already weak signal due to comparatively low biological concentrations of GABA.Accurate frequency alignment of individual transients is also critical to avoid artefacts associated with imperfect subtraction of the sub-spectra.
Further to the competing demands of high SNR, good temporal resolution and reasonable anatomical specificity, MRS measured in the presence of a cognitive task may be confounded by BOLD-related changes in signal relaxation and local shim quality.An increased proportion of oxyhaemoglobin is associated with a longer effective transverse relaxation time (T 2 *) and a correspondingly reduced spectral linewidth.This may affect the separability of certain metabolites or their unambiguous extraction from background signals, and hence represents a potential confounding factor in the obtained concentration estimates. 8Reliable consideration of MRS data acquired under functional conditions is contingent on accurate characterization and consideration of such changes.
While problematic for metabolite modelling, BOLD-related changes in spectral linewidth allow for quantification of the BOLD effect: the strong water peak in a water-unsuppressed MRS acquisition has been shown to reflect localized task-induced haemodynamic changes. 9,10švalka et al. 11 combined water-unsuppressed and water-suppressed PRESS acquisitions in an interleaved manner, to allow robust characterization of the BOLD-related linewidth changes alongside accurate measurement of more subtle changes in the metabolite spectrum.A similar approach has recently been used with semi-LASER at 7 T, finding a non-linear relation between visual contrast, BOLD and the Glu response, 12 but neither BOLD effects on linewidth nor changes in Glu during response inhibition. 13Metabolite cycling 14 and interleaved fMRI-functional MRS (fMRS) (semi-LASER, 3D echo-planar imaging (EPI)) 12,[15][16][17][18][19] sequences have also been demonstrated at higher fields, including reports of significant positive correlation between BOLD and Glu signal changes 15,16 and task-related variation of the excitatory-inhibitory balance. 18However, none of these techniques is optimal for the measurement of GABA.Therefore, we extend the interleaved water-reference approach to a GABA-editing context, with a locally adapted MEGA-PRESS implementation that incorporates unsuppressed water-reference signals at a regular interval within the edited acquisition scheme.
In a typical MRS analysis, final spectra for quantification are often obtained by averaging across a large number of repetitions.In functional contexts with event-related sampling, a general linear model (GLM) approach offers greater flexibility in dealing with confounding factors 20 ; such an approach has been previously applied to modelled concentration estimates 21 and in full 2D fitting methods. 22In the present study, we adopt a novel linear decomposition procedure to extract reconstructed sub-spectra corresponding with functional events or blocks of interest, prior to fitting.
To demonstrate the efficacy of these techniques for simultaneous measurement of GABA, glutamate and glutamine (Glx), and BOLD response, a well known cognitive task based on the Eriksen flanker 23 was used.fMRS data are obtained from the anterior cingulate cortex (ACC), selected as a representative node of the generalized extrinsic mode network 24,25 (EMN), wherein it has a postulated role in cognitive control and inhibition. 26We hypothesized a strong correlation between BOLD assessed with the fMRI and fMRS methods during task performance, and a robust increase in Glx in task-ON versus task-OFF periods consistent with previously reported values, of the order of 5-9%. 27We also anticipated reliable estimates of GABA+, possibly showing a negative association with task performance. 28

| Subject recruitment and demographics
Healthy subjects (N = 81) were recruited from the Bergen local community and at Haukeland University Hospital, through posters and an article in a local newspaper (Bergens Tidende).There were 44 males and 37 females, with mean age of 34.37 ± 11.2 years, range 19-62 years.All potential subjects were screened for history of major head injuries, implanted medical devices, substance abuse, and neurological and medical illnesses.
Written informed consent was obtained from all subjects prior to the study.Subjects were instructed to abstain from caffeine and nicotine in the 2 h prior to scanning.These data were collected as part of a larger study including additional imaging modalities and a clinical cohort, which will be the subject of future analyses.The study was approved by the Regional Committee for Medical Research Ethics in Western Norway (REK Vest 2016/800).

| MR scanning protocol
MR data were collected on a 3.0 T GE MR750 scanner (GE Healthcare, WI, USA), with an eight-channel head coil.The scanning protocol included a high-resolution T 1 structural acquisition: fast spoiled gradient (FSPGR) sequence, 188 sagittal slices, 256 Â 256 isometric 1 mm voxels, 12 flip angle and T E /T R approximately 2.95/6.8ms respectively.GABA-edited spectroscopy data were collected from a 22 Â 36 Â 23 mm 3 (18.2mL) voxel placed medially in the ACC, centred on an imaginary line through the forward part of the pons, parallel with the brain stem (see Figure 1 and Supporting Information Section B.2). Spectroscopy data were collected with a locally modified sequence detailed in Section 2.2, while subjects performed the functional task described in Section 2.3.Subsequently, BOLD-fMRI data were collected as subjects performed a similar functional task, with an EPI sequence, T E = 30 ms, T R = 2500 ms, 90 flip angle, 36 slices of 128 Â 128 voxels (1.72 Â 1.72 mm 2 ), 3.0 mm slice F I G U R E 1 Placement of the fMRS VOI for the healthy subjects, mapped to standard space.Contours indicate [5, 50, 95]-percentile coverage of the achieved placement across subjects.Dashed lines in A illustrate landmarks used for voxel positioning: medial ACC, centred on an imaginary line through the forward part of the pons (red), parallel with the brain stem (indicated in blue).thickness with 0.5 mm gap (3.5 mm slice spacing); 240 volumes for a total acquisition time of 600 s.Note that gradient-heavy fMRI acquisitions were performed after the MRS, to minimize the impact of thermal drift. 29A summary of key sequence and hardware parameters is presented in an MRSinMRS 30 checklist, Table S1.

| MEGA-PRESS sequence adaption
To facilitate simultaneous acquisition of time-resolved MRS data and robust assessment of BOLD signal changes in response to the cognitive task, local extensions to the MEGA-PRESS sequence were adopted.These build upon the standard GE MEGA-PRESS implementation, adding per-T R trigger pulses (500 μs after the end of the excitation RF pulse) to allow for precise synchronization with an external stimulation paradigm.Moreover, the local implementation allows the chemically selective saturation (CHESS) water suppression pulses to be disabled at a user-specified interval-thereby periodically acquiring a water-unsuppressed reference (WREF) signal within the regular GABA-editing sequence.
The paradigm was implemented in E-Prime 2.0 SP1 (2.0.10.353:Psychology Software Tools, Pittsburgh, PA, USA, https://pstnet.com/),and presented from a separate PC adjacent to the MR system console.Stimuli were presented through a set of goggles mounted on the head coil (NordicNeuroLab (NNL), Bergen, Norway, http://nordicneurolab.com/,note declaration of interest) in light grey on a black background, with relatively small font and spacing to remain near the parafoveal field of view 31 ; see Figure S1d.A central fixation cross (Figure S1c) was presented between trials and during task-OFF blocks.The subjects used a set of response grips (NNL) to indicate the direction of the central 'target' arrow.
Instruction was provided in Norwegian or English depending on each subject's preference (see Figure S1a,b).
The paradigm was run as a standard block-event design, beginning with a 60 s task-OFF block followed by alternating 30 s task-ON blocks and 60 s task-OFF blocks; there were 11 task-ON blocks for the fMRS acquisition, six for the fMRI.This is summarized in Figure 2C.Within each task-ON block, one trial was presented per T R (average interstimulus interval (ISI) of 1500 ms), giving a total of 220 trials for the fMRS acquisition, and 120 trials for the shorter fMRI acquisition-of which a randomly selected 40% were incongruent.
Trial onset timing was jittered randomly such that stimuli were presented T S-A = 100-350 ms before the excitation pulse, spanning the range examined by Apšvalka et al. 11 and most of the evoked gamma-band power changes reported by Lally et al. 32 Stimuli for each trial were presented for 350 ms, with a nominal 800 ms response window (measured from the beginning of presentation).Trigger signals were received by E-Prime via an NNL SyncBox.Since stimulus onset preceded receipt of triggers, a feedback loop was used to iteratively tune compensatory factors-as such, the achieved intervals deviated slightly from nominal.For the fMRI acquisition, triggers were only issued at the beginning of each new task-ON or task-OFF block.
Individual response data were processed using in-house tools (see Section 2.6).Missing and ambiguous responses were flagged for exclusion before evaluation of basic task performance metrics: reaction time (RT) and percentage response accuracy (RA).This was done separately for each subject, for each functional scan (fMRS and subsequent fMRI) and for each flanker condition (congruent versus incongruent).Achieved ISI and proportion of congruent versus incongruent stimuli were also assessed to verify correct execution of the paradigm itself.

| fMRS data processing and quantification
Functional MEGA-PRESS data were initially loaded using the GannetLoad function from Gannet Version 3.1, 33 although bypassing that software's standard processing pipeline.After import, eddy-current correction 34,35 and global zero-order phase correction using a dual-Lorentzian model for Cr and choline (Cho) were performed.Robust spectral registration 36 was applied independently for edit-OFF and edit-ON transients, and for each phase cycling step, before alignment across the sub-spectra by spectral registration 37 -all using basic Gannet functionality.Individual, uncombined transients were taken for further processing.

| Metabolite sub-spectra: a linear model for spectral combination
The conventional approach for extraction of sub-spectra involves averaging across transients, which relies on assumptions that can be difficult to uphold with sparse event-related sampling and offers limited flexibility in dealing with confounding factors.As such, a linear model was used for this purpose.A linear model of the form Y = XB + U was constructed, wherein each acquired transient (Y n,* ) is expressed as a combination of a normal, baseline edit-OFF spectrum (X *,1 ), with periodic contributions associated with editing (the baseline DIFF spectrum, X *,2 ), and nuisance factors: periodic covariates associated with each of the two phase-cycle steps (X *,3:4 ), and additional covariates partitioning the data according to inferred motion (X *,4 + [1:n_partitions] ).In this context, probable subject motion is inferred from abrupt changes in the estimated water frequency, as detailed in Supporting Information Section B.3.
Functional events were binned according to the achieved interval between stimulus and acquisition (T S-A ).Five bins (bin 1-5 ) were defined, with edges at [100, 183, 267, 350] ms, open at either end, lower limits inclusive; the inner three bins evenly cover the nominated 100-350 ms T S- A range.This resulted in approximately 48 task-ON metabolite transients per bin (with some individual variation), and 320 task-OFF metabolite transients.These bins were incorporated into the design matrix as columns X *,e1:e5 .Note that columns (3:e5) capture variability common to the edit-OFF and DIFF (sub-)spectra; these columns were subsequently replicated and masked to capture differential changes associated with the respective variables; corresponding differential masked event columns are designated X *,(e 0 1:e 0 5) .A representative design matrix for a single subject is presented in Figure 3A.
Two further variations of this design matrix were used, to extract sub-spectra corresponding with different functional stimuli, and with different periods within the task-ON block.For the former, the event-related columns X *,(e1:e4) instead selected transients according to the stimulus and response: congruent accurate, congruent inaccurate, incongruent accurate, incongruent inaccurate.Approximately 88 and 59 task-ON metabolite transients were available for the congruent and incongruent components respectively, with subdivision according to accuracy dependent on individual performance.For temporal subdivision of the task-ON blocks, event-related columns X *,(e1:e5) selected successive 7.5 s portions of all task-ON blocks ([0-7.5, 7.5-15, 15-22.5, 22.5-30] s), and the first part of the following task-OFF block (30-37.5 s).
The system was solved stepwise.
1.An initial fit was performed on only the components of interest (edit-OFF, DIFF and masked event bins), using the MATLAB decomposition function for complete orthogonal decomposition to yield a minimum norm least squares solution.3. Removal of 'bad' transients was performed based on the residual error after fitting, where the Z-score of the residual error for a given transient exceeded 2.5.

Fitting
Steps 1 and 2 were repeated until the relative improvement in norm_residuals was below 10 À8 . 5. Baseline edit-OFF and DIFF spectra were extracted directly from the corresponding components X *,(1,2) , while edit-OFF and DIFF components binned according to functional events (bin 1-5 ) were evaluated arithmetically from components X *,(e1:e5) and X *,1 , and X *,(e 0 1:e 0 5) and X *,2 , respectively.6. Linewidth matching was performed with a Lorentzian filter, using Cho and Cr from the associated edit-OFF sub-spectra to assess linewidth.
Line broadening was performed to match the maximum measured linewidth across all extracted spectra for the individual subjects, capped at 130% of the median to avoid matching to strong outlier targets.
7. Singular value decomposition (SVD) was used to remove residual water signal.
8. Extracted spectra were fit using Gannet's GABAGlx model, which fits GABA+ with a Gaussian peak around 3.02 ppm and Glx as a pair of Gaussian peaks around 3.71 and 3.79 ppm-both measured from the DIFF spectrum. 38Water-scaled estimates with tissue-class correction 39 are reported.

| Quality control
Three rejection criteria were applied to the extracted DIFF spectra, in series.
• R1 rejects fits where the full-width at half-maximum (FWHM) linewidth of the fitted GABA+ or inverted NAA peak exceeded 30 or 12 Hz, respectively.
• R2 rejects fits where the SNR of the fitted NAA peak is extraordinarily low, below 20 (noting that quite low SNR may be expected in some of the extracted event-related spectra, per Section 2.4.1).
• R3 rejects extreme outliers: GABA+ or Glx estimates differing from the median by more than five times the median absolute deviation (MAD) across all spectra surviving R1 and R2.

| Unsuppressed water signal: T 2 * and BOLD assessment
Unsuppressed water-reference spectra (WREF) were fitted with a pseudo-Voigt function superimposed on a linear baseline, filtering for outliers and discontinuities as detailed in Supporting Information Section B.4.A linear model, similar to that described for the metabolite signal in Section 2.4.1, was constructed.Unit impulses at the recorded functional event onsets were convolved with a dual-gamma haemodynamic response function (HRF) model to create a model for the expected BOLD signal, (X *,1 ), with covariates for phase cycling step (X *,2:3 ), and inferred subject motion (X *,3 + [1:n_partitions] ).The model was fitted with linear least squares regression, with the BOLD model coefficient Y 1 (roughly equivalent to ΔFWHM water , the overall change in water linewidth) reflecting the strength of the individual's BOLD response as assessed from WREF signals during the fMRS acquisition.

| fMRI data processing
fMRI block analysis was performed using FEAT (FMRI Expert Analysis Tool Version 6.00, part of FSL).Processing included motion correction using MCFLIRT, 40 masking of non-brain data using BET, 41 spatial smoothing (5 mm FWHM Gaussian kernel) and linear registration to a high-resolution standard space structural template using FLIRT, 40,42 subsequently refined with non-linear registration using FNIRT. 43,44The MNI-152 (non-linear, sixth-generation) template 45 was used as a registration target.Time-series statistical analysis was performed using FILM with local autocorrelation correction 46 ; the resulting statistical images in subject space were thresholded non-parametrically 47 using clusters determined by Z > 3.1 and a corrected cluster significance threshold of p cluster = 0.05.Individual registration outcomes were subject to visual inspection, and higher-level (group-average) statistics evaluated after transformation into standard space; group-average response allowed visual inspection of the task response and selection of control volumes of interest (VOIs), but was not subject to further quantitative analysis.
Per-subject VOIs were defined using the individually prescribed fMRS voxel geometry in the ACC (VOI fMRS,ACC ).Based on the group-average fMRI response in standard space, additional VOIs with volume similar to the fMRS acquisition (around 18.2 mL) were defined in the temporal pole
Summary statistics are reported using outlier-resistant methods (median and MAD).Where appropriate, compatibility of variance was assessed using the Fligner-Killeen method, 56 and normality was assessed with the Shapiro-Wilk method. 57In cases where Shapiro-Wilk indicated non-normality, the Wilcoxon signed rank test 58 was used for hypothesis testing between related samples.Confidence intervals derived from parametric bootstrapping (10 000 permutations) are denoted CI boot,95% .In cases where adjustment for multiple comparisons was performed, Holm-Bonferroni 59,60 correction was applied within the sub-analysis; adjusted p values are denoted p holm , with a corrected significance threshold defined as p holm < 0.05; uncorrected p values are denoted p unc. .Outcomes for all performed tests are supplied in Supporting Information Sections C.3-C.5 (Tables S6-S14).Where correlation coefficients are reported, these are robust Spearman correlations using the skipped method 61,62 to exclude bivariate outliers.

| Behavioural outcomes
Basic behavioural outcomes from the flanker task are summarized in Figure 4 and Table S2.Strongly significant differences ( p holm < 0.001) were seen between congruent and incongruent trials for the RA, RT and RA/RT parameters, during both the fMRS and fMRI task-ON blocks.
Comparing the fMRS and subsequent fMRI task performance, increases in RA and RA/RT were observed in the latter, which were significant ( p holm < 0.01) for incongruent trials.

| BOLD assessment by fMRI and fMRS
Group-average spatial response for the fMRI task is presented in Figure 5 S3-5), with additional slices shown in Figure S2.
A significant correlation was seen between the strength of the BOLD response as assessed with the fMRS and fMRI methods (see Figure 6) within the individually prescribed VOI fMRS,ACC (r = 0.34, p holm = 0.01, CI 95% [0.132, 0.536]) and within the fixed ACC region (not accounting for individual placement) VOI ACC (r = 0.33, p holm = 0.01, CI 95% [0.11, 0.52]).No significant correlations were observed for control regions in the precuneus, occipital or temporal pole (see Table S6).

| fMRS
From a total of 324 spectra derived from the per-subject model, rejection criterion R1 removed 25 fits (of which eight exhibited high NAA FWHM and 19 exhibited high GABA FWHM; some overlapped).Subsequently, one spectrum was removed by R2 (low NAA SNR) and eight extreme outliers were rejected by R3 (four for GABA+, five for Glx; some overlapped).Achieved spectral quality metrics (SNR and FWHM) and group-average metabolite estimates after fitting and quality control are presented in Table 1; resultant spectra and mean fit for task-OFF and task-ON conditions F I G U R E 5 Group mean Z statistics from the fMRI task, in MNI-152 standard space; red-yellow positive, blue-cyan negative Z-scores on the range. 7,25I G U R E 6 Relation of BOLD assessed by linewidth changes in the fMRS acquisition to that observed from the fMRI data regionally masked to the individual fMRS voxel (A) and a control region in the temporal pole (B), showing significant correlation in the former and no correlation in the latter.
T A B L E 1 Quality metrics and concentration estimates from the fMRS analysis, task-ON versus task-OFF.F I G U R E 7 Metabolite concentration estimates for task-OFF versus task-ON conditions (A), and from spectra recorded at varying times after event onset, T S-A (B).
are additionally presented in Figure 8 below.To verify the efficacy of linewidth matching, the standard deviation of Cho/Cr FWHM across extracted time bins (and task-OFF) for each subject was evaluated before and after matching: median SD before matching 0.165 ± 0.05 Hz, reduced to 0.049 ± 0.021 Hz after matching-a significant improvement, p holm < 0.001.
Assessing the discrete T S-A bins (Figure 7B, Table S8), while Glx appeared somewhat elevated in each time point relative to task-OFF, this was significant only for T S-A = 100-183 ms (14.4% increase, p holm = 0.02, CI boot,95% [1.92, 17.3]).Pairwise comparison between time bins revealed no significant differences (all p holm > 0.29; marginal p unc = 0.049 between T S-A = 100-183 and 267-350 ms, see Table S9), and no differences were seen for GABA+, either with respect to task-OFF or between time bins (see Table S10).For spectra modelled on temporal subdivision of the task-ON block, no significant differences were seen across timepoints within the block (see Table S11).Comparing metabolite estimates obtained for different task stimulus and response conditions (congruent versus incongruent, accurate and incongruent, inaccurate) also found no significant differences (see Table S12).

| Integrated outcomes
No significant correlations were seen between baseline metabolite concentration estimates (task-OFF) and BOLD response measured by either method; see Table S13.An exploratory analysis comparing behavioural measures (RT, RA, RA from incongruent trials only [RA incong ] and incompatibility slowing: RT incongruent À RT congruent ) with BOLD and with metabolite estimates for GABA+ and Glx in task-OFF and task-ON periods is detailed in Table S14; this revealed a significant correlation between incompatibility slowing and the strength of the fMRS-assessed BOLD response (r = 0.39, p holm = 0.049, CI 95% [0.151, 0.586]).A possible negative correlation between RA incong and task-OFF Glx levels (r = À0.379,p holm = 0.15, p unc = 0.006, CI 95% [À0.594,À0.112]) did not survive strict correction for multiple comparisons.

| Behavioural outcomes
The behavioural outcomes at group level showed strong incompatibility slowing effects, consistent with existing studies, 31,[63][64][65] and a corresponding decreased RA.Improved RA between the fMRS and fMRI tasks suggests learning effects; see Section 4.2.Absolute RTs and error rates are compatible with some prior studies, 63 although there is substantial variability in the literature.Therefore, we conclude that the task was performed effectively (at group level) and is expected to have elicited the desired cognitive load as evidenced in the BOLD outcomes.
F I G U R E 8 Mean extracted spectrum from the fMRS sequence for task-OFF, showing ±1 SD range and mean fit for both task-OFF and task-ON conditions; note substantial overlap for GABA+ around 3.02 ppm (□) and slight differences in baseline fit around 3.02 ppm (◇) and in Glx shape at 3.71 () and 3.79 ppm (▷).

| BOLD functional outcomes
Group-average BOLD contrast evaluated from the fMRI data exhibited a characteristic 'task positive' structure consistent with a generalized EMN 24,25 during task-ON blocks (see Figure 5), including significant activation in the ACC region targeted by the fMRS voxel placement.Network nodes typically associated with the 'resting state' default mode network 66 showed notably increased activation in the task-OFF blocks.In this way, we conclude that the paradigm resulted in the anticipated patterns of brain activation and deactivation in response to the presence or absence of task stimuli.
Although a significant correlation was seen between the BOLD signals assessed by the fMRI and fMRS methods, this correlation remained moderate.This could in part be attributed to the inherent variability of each measure, and intra-session variability in the BOLD response, with learning effects (indicated by behavioural outcomes) and fatigue perhaps mediating the strength of the BOLD signal in the later acquisition.While a counter-balanced design may mitigate these factors, technical considerations relating to gradient heating and thermal drift 29 precluded counterbalancing in the current study.

| GABA
A recent meta-analysis covering fMRS studies of GABA and Glu/Glx 67 shows small effect sizes and large heterogeneity, both in experimental design and outcomes.Only one study in the meta-analysis investigated GABA changes during a cognitive (Stroop) task in the ACC, showing negative correlation between GABA change and BOLD signal change during task conditions. 28Outcomes across other tasks and locations [68][69][70][71][72][73][74] varied substantially, with no significant change found for GABA relative to baseline across studies.In this context, our lack of significant results when examining GABA+ changes in relation to task is unremarkable.
Similarly, while an earlier meta-analysis 75 shows a number of studies reporting negative correlation between baseline GABA+ levels and BOLD magnitude during task performance, [76][77][78][79][80] findings of positive associations (or non-associations) have also been reported, [81][82][83][84] suggesting a somewhat more nuanced relationship.Once again, our lack of significant results is unremarkable in this context.
The T R parameter adopted in the present study was somewhat lower than commonly used (often around 2000 ms; References 85 and 86).
This implies increased T 1 weighting, which may increase the contribution of macromolecule signals to the measured GABA+ due to their substantially shorter T 1 . 5,87While this trade-off allowed a higher number of functional events to be acquired, altered weighting and increased macromolecule contribution may have limited the sensitivity to any more subtle GABA changes.
Since selective refocussing of coupling to GABA spins around 3 ppm does not affect that peak uniformly (predominantly affecting the outer lobes), it is not clear that the DIFF spectrum modelling will be equally sensitive to transient, task-related changes affecting the edit-ON and the edit-OFF parts.Explicitly modelling this, and tailoring the experimental design accordingly (perhaps targeting more transients from the more responsive sub-spectrum) may yield stronger outcomes.
Finally, incorporating real-time frequency adjustment or breaking the long fMRS acquisition into shorter segments, periodically adjusting the centre and editing frequencies and perhaps re-shimming to account for any changes over time (due to thermal drift, 29 subject motion and so forth) may yield better editing outcomes for future studies.This factor would need to be incorporated into the model for spectral combination (Section 2.4.1), and would additionally benefit co-edited signals such as Glx.

| Glx dynamics
The meta-analysis by Pasanta et al. 67 shows moderate effect sizes for Glu and Glx across studies, while the meta-analysis by Mullins 27 reports a mean change in Glu of 6.97% CI 95% [5.23, 8.72]-strongly dependent on experimental design.For the few studies reporting basic cognitive tasks in the ACC, 28,88,89 typical changes were closer to 3%.The changes reported herein (8.85% increase, CI boot,95% [4.83, 13.9]) are comparatively strong; it is unlikely that a change of this magnitude, on the observed timescale, could be explained by metabolic processes (such as Glu synthesis) alone.
A more plausible explanation would be a compartmental shift.It has been proposed that a substantial portion of Glu 90 and Gln 91,92 within the neuron may exist in pools where metabolite movement and tumbling is restricted-putatively including the synaptic vesicles.This restricted movement leads to higher T 2 relaxation rates, and hence reduced visibility of the associated MRS signal, particularly at longer echo times such as used in the present study.During neural activity, Glu released from vesicles may move to a compartment where it is more visible, 27,93 such as the cytosol or synapse, leading to an increase in the apparent, measured signal.
As well as being the primary excitatory neurotransmitter, Glu plays a key role in normal energetic processes of neural cells, and in the synthesis of GABA. 94,95The proportion of Glu involved in each process is unclear, and current MRS techniques are generally not sufficiently selective to distinguish between the different compartments-although differing T 2 across compartments means the relative visibility may be modulated by T E , with higher T E values being more sensitive to compartmental shifts. 27,93Disentangling these factors would serve to validate the hypothesized compartmental shift, and improve interpretability of the outcomes.
While the concentration estimates showed a robust increase in measured Glx between task-OFF and task-ON spectra, we note (with reference to Figure 8) that the 3.71 ppm sub-peak appears paradoxically slightly reduced in amplitude, while the most notable increase is seen in the leftmost part of the 3.79 ppm sub-peak, with an apparent broadening of the Glx peak in that area.While subtle variations in the peak shape may be consistent with a compartmental shift, one might expect a compartment of lower T 2 (increased MRS visibility) to yield somewhat sharper peaks-this was not evident in the present data.Another possible explanation is that the apparent increase at the leftmost edge of the 3.75 ppm C2 peak may reflect an increased contribution of glutamine (Gln) to the observed signal, with the Gln C2 peak at 3.764 ppm being 0.017 ppm to the left of the Glu C2 peak at 3.747 ppm.
The concentrations of Glu and Gln are closely linked, with the two in constant exchange in the neuronal-glial Glu-Gln cycle. 95,96While independent quantification of Glu and Gln may be feasible from GABA-edited data of sufficient quality, 97 this approach has not been broadly adopted.
Given the comparatively low SNR of fMRS data, no attempt was made to separate Glu and Gln in the current data.Nonetheless, this may be feasible in further analysis using a basis-set fitting approach and strict acceptance criteria, or in future studies at higher field strengths and/or with dedicated sequences and judicious choice of subecho timings. 98e typical trajectory of the Glu response over time remains uncertain; while evidence from block-related designs suggests a gradual increase (around 16-20 s from the start of a block), with compartmental shifts probably relating to increased metabolism associated with neural activity, event-related studies suggest a separate process on a shorter timescale (returning to baseline within 3-4 s), related directly to neurotransmitter release from the vesicles. 11,20,27,32For the event timing examined in the present study (T S-A nominally 100-350 ms), while statistically robust changes were only observed in the earlier time section, we do note a slight tendency towards baseline in subsequent times, as well as greater variance in the earlier part-perhaps suggestive of substantial interindividual variability in the timing and onset of the initial response, which may normalize later in the response.Although the block timing of our study is well suited to elicit a strong BOLD response and to assess the putative slower glutamate response function (GRF), with long task-OFF periods to allow for a robust return to baseline, the event timings (ISI $ 1500 ms) were not optimal for unambiguous separation of adjacent events if we assume a return to baseline after 3-4 s for the faster GRF.Future studies primarily intending to model the GRF should ensure a longer interval between successive events, perhaps coupled with a longer T R for more optimal GABA measurement.
While the present analysis makes no particular assumptions regarding the shape of the GRF (beyond the restricted T S-A range), we note that recent studies have begun to investigate putative metabolite response functions 99 ; once an appropriate response function is determined, incorporation of this into the present model (in place of the coarse binning) will be trivially accomplished, and may well improve sensitivity to subtle variations.Furthermore, while the present model reconstructs discrete sub-spectra which are fitted independently using existing 1D modelling tools, we note recent developments towards simultaneous, 2D spectral-temporal fitting 100 of multiple spectra linked by an arbitrary model, such as that implemented in the FSL-MRS dynamic fitting module. 22,101With appropriate constraints, incorporating our linear model into such an approach may further improve fitting performance, particularly for lower-SNR cases.
While Glu changes were observed in relation to the presence or absence of task, Glu was not found to modulate significantly with increasing inhibition demands, which may be expected to be specific to incongruent flanker trials.This absence of modulation between lower-and higherdemand cases is in contrast to the findings of Ip et al., 16 wherein significant Glu changes were only observed in the highest-contrast cases (for visual stimulus, in the primary visual cortex V1), with sensitivity limits implicated in that study.Although further investigation would be required, these findings may suggest a possible non-linearity in the Glu response to task intensity-compatible with recent reports for visual contrast. 12

| CONCLUSIONS
In the current study, we present a MEGA-PRESS sequence adaption suitable for the concurrent measurement of time-resolved GABA+, Glx and BOLD from a single voxel, at a ubiquitous 3 T field strength.We additionally present a novel linear model for extracting spectra modelled from functional stimuli, building on well established processing and quantification tools.With these tools, we demonstrate a robust increase in measured Glx concentration in the ACC in response to a task stimulus, measured concurrently with regional BOLD response.The latter is shown to significantly correlate with the BOLD response as assessed by traditional fMRI methods.Findings for Glx are consistent with theoretical models and with our hypothesis, and both fMRS and fMRI findings align with existing literature.While identifying a number of readily achievable optimizations for future usage (particularly with regards to GABA sensitivity), we conclude that the acquisition and analysis methods documented herein are effective for the concurrent measurement of GABA, Glx and BOLD, in relation to a functional task.
These parameters ensured that the sequence remained balanced with respect to the number of edit-ON, edit-OFF and WREF signals acquired at each phase-cycle step, and with respect to any carryover impact of residual, unrelaxed signals (particularly from the water-unsuppressed acquisitions) on subsequent transients, with equal proportions of edit-OFF and edit-ON sub-spectra at each phase-cycle step being preceded by a WREF acquisition.F I G U R E 2 Timing of the fMRS acquisition.A, A MEGA-PRESS sequence modified such that every third transient is acquired without water suppression (denoted REF, shown in blue).B, Synchronization to task stimuli is achieved with additional trigger pulses, with the interval from stimulus onset to acquisition T S-A jittered between 100 and 350 ms.C, Flanker stimuli are grouped into 30 s task-ON blocks between 60 s task-OFF blocks.

2. 3 |
Functional paradigm-flanker task Subjects performed a cognitive task based on the Eriksen flanker task, 23 using visual arrow glyphs.In each trial a central 'target' arrow is presented, surrounded by four 'flankers'.The task is to indicate the direction of the central target arrow.If the target and the flankers match

2 .
Covariate components (nuisance factors) were modelled to the residuals from Step 1.

F
I G U R E 3 A, B, Sample design matrix (A) and partitioning according to motion inferred from the water frequency (B), for the first 5 min of data from a single representative subject.C, Extraction of spectra from this model (see Section 2.4.1).

(
VOI temporal around [À38.7, 10.6, À33.6] mm), the precuneus near the posterior cingulate (VOI precuneus around [0.5, À56, 20] mm), the occipital (VOI occipital around [2.83, À93, 12.3] mm) and the ACC (VOI ACC around [0.5, 3.5, 45.5] mm) in standard MNI-152 template space, subsequently transformed to the individual subject geometry.These VOIs are shown in Figure S2.Median Z-score across each VOI was extracted, without thresholding.Subjects where the median Z-score within the spectroscopy voxel was not positive were considered non-conformant: either the individual's task performance did not elicit the expected functional response, or the achieved positioning of the fMRS voxel did not demonstrably capture this response; see Figure 6A below.
, showing strong task-related activations in frontotemporoparietal regions and the ACC/SMA, with parietal/precuneus and ventromedial inferior frontal regions more activated during task-OFF periods.Cluster F I G U R E 4 Outcomes from the behavioural task, showing significant increase in RT (A) and reduction in RA (B, C) for incongruent conditions, and slightly improved accuracy in the second (fMRI) session.statistics, location of local maxima, and corresponding atlas structures are presented in Supporting Information Section C.2 (Tables