Benchmarking GABA Quantification: A Ground Truth Data Set and Comparative Analysis of TARQUIN, LCModel, jMRUI and Gannet

Many tools exist for the quantification of GABA-edited magnetic resonance spectroscopy (MRS) data. Despite a recent consensus effort by the MRS community, literature comparing them is sparse but indicates a methodological bias. While invivo data sets can ascertain the level of agreement between tools, ground-truth is required to establish accuracy, and investigate the sources of discrepancy. We present a novel approach to benchmarking GABA quantification tools, using several series of phantom experiments with iterated GABA concentration. Each series presents a different set of background metabolites and environmental conditions allowing comparison of not only individual estimates, but the ability of tools to characterise changes in GABA across a range of potential confounds. The methodology of the phantom experiments is presented, as well as characterisation of the data. We also perform an initial comparative analysis of several common MRS quantification tools and in-house code, to illustrate utility of the dataset and the potential bias introduced by different quantification methods. The GABA-to-NAA ratios reported by each tool are compared to the ground-truth, and estimation accuracy is assessed by linear regression of this relationship. While the linearity of GABA-to-NAA gradients is generally captured by all of the tools, a large variation in the slope, offset and environmental stability of the gradient is observed. The primary driver of differences in linear combination modelling is the choice of basis function. However, tools employing a common basis and pre-processing still produce differences on the order of 4%. Less-strictly parametrised fitting approaches appear to improve the robustness of quantification, but accurate modelling of edit efficiency calculations is still necessary to avoid systematic offsets. In general, the level of variation suggests that comparisons across...


INTRODUCTION
Gamma-aminobutyric acid (GABA) is the primary inhibitory neurotransmitter in the central nervous system, maintaining the brain's excitation-inhibition balance. Its prominent role in neurotransmission and metabolism has led to extensive study and a plethora of applications for its detection. Disruption of GABAergic processes has been observed in Schizophrenic patients 1,2,3 and GABA receptor dysfunction has associations with epilepsy 4,5 . GABA processes are prominent in type I diabetes 6 and autism spectrum disorders 7 . GABA is also the subject of intense study by the psychology community, with GABA levels influencing impulsivity 8,9 , drug addiction 10 , anxiety disorders 11 and depression 12 . The most common method for non-invasive detection of GABA in-vivo is magnetic resonance spectroscopy (MRS). However, at clinically relevant field strengths the GABA signal is overlapped by more prevalent metabolites, necessitating other detection methods for reliable quantification.
While specialised techniques such as 2D J-spectroscopy 13 are used, the most common method is edited MRS 14,15,16,17 , in particular MEGA-PRESS 18,16,19,20 . MEGA-PRESS GABA-editing utilises two interleaved acquisitions -'edit on' and 'edit off' -which differ only in the frequency of applied editing pulses. The edit frequencies are chosen to take advantage of the J-evolution of the target metabolites, such that upon subtraction of the two acquisitions, only resonances within the edit bandwidth and those coupled to them remain. The widespread use of edited MRS has seen the development of a range of software packages aimed at processing and quantification of the spectra 21,22,23,24,25,26,27 . The tools process data by some combination of apodization, frequency and phase calibration, macromolecular baseline subtraction and residual water signal removal.
Each step may be realised by a diverse array of methods, providing many potential sources of variation for the quantified signal. This effect is compounded by further differences in the quantification procedure. The most common approach to quantification is the generation of basis spectra -via phantom scans or more commonly by quantum mechanical simulation -followed by the application of a fitting algorithm to decompose the MRS signal into weighted combinations of the basis functions. For edited MRS, with its simplified spectra, other fitting approaches -using Gaussian/Lorentzian modelling of resonances and numerical integration -become more feasible.
A recent large-scale effort has been made to establish consensus in the MRS community 28,29,30,31,32,33,34 . While recommendations were made around the quantification of MRS data -Choi et al. 29 reference basis set methods and Gaussian peak fitting for GABA quantification -there was no consensus on the best approach. While it has been previously argued that different methods, suitably applied, exhibit similar variation and clinical observations 35,36 , it is acknowledged in the literature that there is a procedural dependency in MRS analysis 37,38,39,40,41 , which provides a statistically significant variation in reported metabolite amplitudes or relative concentrations 42,43 . With the exact pre-processing, simulation and fitting parameters only partially reported in many research articles -an issue discussed by Lin et al. 30 -the methodological dependence of MRS quantification propagates, reducing the reproducibility of MRS studies. Despite this, there has been relatively little comparative study and benchmarking of the MRS analysis tools.
The goal of this study was to generate benchmarking data sets for the purpose of assessing the relative performance of GABA quantification tools and investigate the methodological dependence of GABA-edited data analysis. Previous work in this area has utilised either in-vivo data -for which there is usually no ground truth to evaluate the results -or simulated data sets, which inevitably cannot capture all sources of experimental variation. For this study, MRS benchmarking data sets were generated using several series of calibrated phantom experiments. While phantoms fail to capture all influences of the in-vivo environment, they do allow the inclusion of many experimental factors, while still enabling a ground truth comparison and ability to control the environmental characteristics expected to influence results. Four series of phantom data were experimentally acquired. In each case, the concentrations of various metabolites -creatine (Cr), N-acetyl-L-aspartic acid (NAA), glutamine (Gln) and glutamate (Glu) -are fixed at approximately in-vivo levels, while the concentration of GABA is iteratively increased. At each GABA concentration step, MEGA-PRESS spectra were acquired. The four series of experiments were designed to increase the complexity of the spectra, starting with pH and temperature calibrated solutions containing only NAA, Cr and GABA (E1), then adding Glu and Gln to the solutions to increase the complexity of the spectra (E3), and finally progressing to a series of gel phantoms which more closely mimic in-vivo spectra 44 (E4 1 ). To investigate the influence of the molecular environment, a series of non-neutral pH solutions of NAA, Cr and GABA (E2) were also acquired. The phantom data set has been designed with the benchmarking of new quantification tools in mind, for example in the use of convolutional neural networks for GABA quantification 45 . As such, the full data have been made available for others to benchmark their own analysis software, or optimise existing tools (DOI:10.21227/ak1d-3s20) 46 .
Furthermore, to demonstrate the utility of the dataset, MEGAPRESS spectra were quantified using a variety of state-of-the-art MRS analysis software to investigate the potential methodological bias that could be introduced. The GABA-to-NAA concentration ratio was estimated for each tool by applying appropriate signal attenuation corrections where needed. The reported ratios are plotted against the known ground-truth, and fit using linear regression. The GABA-to-NAA gradients are used to investigate the systematic offsets that can be introduced into reported concentration changes, and to inform future comparative benchmarking studies.

METHODS
Phantom preparation procedures, experimental setup and scan parameters, and details of the analysis are described.

Phantom preparation
The  Table 1 and representative images of the phantoms are shown in Figure A1.
The solution series were prepared by dissolving the required concentration of background metabolites in de-ionised water to create a base solution. 290 ml of the GABA-free base solution was transferred to a round bottom glass flask, filling it to the neck. To avoid dilution of the background metabolites, a small amount of this base solution was removed from the flask using a syringe and GABA was added to prepare a concentrated GABA solution (1 mg/ml). To vary the GABA concentration between scans, a small amount of the solution was removed from the flask with a syringe and replaced by the same volume of the concentrated GABA solution. The amount of solution replaced in each step varied depending on the desired GABA concentration change, but was generally between 0.5 ml and 2 ml. This procedure was repeated several times, iterating GABA to a relatively high concentration with the same scan protocols applied for each concentration step. This procedure limited experimental errors solely to the GABA concentration. The pH of the solution was monitored between scans and adjustments were made using a 36% hydrochloric acid (CAS-7647-01-0) solution and a 3.99% sodium hydroxide solution (CAS-1310-73-2) to maintain a pH of 7.2 ± 0.2. The amounts of HCl or NaOH added after the initial pH calibration of the base solution were marginal (≤ 1 ml for the 290 ml flask) and any dilution of the overall solution was deemed negligible. The E2 series was one exception to this procedure, where the pH was intentionally not adjusted, to investigate the effect this had on the various quantification procedures. The final pH of E2 was measured to be 3.0 ± 0.2.
For the gel series, the phantoms were made in advance of scanning. 800 ml of a base solution of Glu, Gln, NAA and Cr was prepared and divided into eight 100 ml portions. Different amounts of GABA were added to each solution and pH calibration was performed as above. Finally, 1 g of agar (CAS-9002-18-0) 2 was added as gelling agent. The mixtures were then heated to 90°C to 100°C while being stirred until the agar had fully dissolved. A small hole (< 2 mm diameter) was created in a spherical plastic mould using a plastic welding tool and the solutions were injected via a syringe. The arrangement was then allowed to cool over night. Once solidified, the gels were examined and the opening was sealed using a small amount of silicone sealant. Spherically-shaped phantoms were deemed the most suitable to reduce magnetic susceptibility-induced field inhomogeneity and minimise the spectral linewidth.
The accuracy of the scale used for mass measurements is 1 mg and the graduated cylinders used for the volume measurements have an accuracy of 1 ml. Assuming negligible errors in the manufacturer's stated molar masses, the uncertainty ∆c in the molar concentration c = m M V , where m is the mass in grams, M is the molar mass of the solute and V is the volume of the solvent, is ∆c c = ∆m m 2 + ∆V V 2 . Assuming (generously) an uncertainty of 2% for the mass measurements and 1% for the volume measurement, the uncertainty in the concentrations of base metabolites is estimated to be less than 3%. The uncertainty in the amount of solution added or removed from the flask in each step via 5 ml syringes is slightly higher and could be up to 10%, assuming an uncertainty of 0.1 ml per 1 ml to account for slight misreadings of the plunger position and tiny air pockets within the syringe.

Scan Protocols and Data Acquisition
All MR scans were conducted at Swansea University's Clinical Imaging Facility using a MAGNETOM Skyra 3 T (Siemens Healthcare GmBH, Erlangen, Germany). The scanner room was temperature controlled to 20 ±0 6°C. Between each scan, the phantom was removed to update concentrations, so the influence of heating is deemed negligible. Signal acquisition was done using the built-in spine coils, specifically the four channel spine coil element 'SP2'. The spine coils provided the platform for the most reproducible set-up and SP2 exhibited the highest SNR of the coil elements available. The phantom was aligned with this element and raised to isocenter, using a cardboard phantom holder to ensure maximum field homogeneity for the scans. Double-spin-echo field maps were acquired for each phantom to assess the homogeneity, and manual shimming and frequency calibration was performed to optimise the spectral width. Further manual calibration was performed to adjust the TX reference voltage and water-suppression pulses.
A 20 × 20 × 20mm 3 voxel at the isocenter of the magnet was selected for single voxel spectroscopy. The GABA-edited spectra were acquired using the MEGA-PRESS pulse sequence WIP859G 3 with CHESS water-suppression 47  For each experiment, the raw data acquired, the single-average coil-combined spectra and combined average edit-on, edit-off and difference spectra produced by the vendor-supplied spectroscopy software have been made available 46 . For the purpose of this paper, the combined channel and average dicom spectra were used. While raw data is generally preferable in order to appraise frequency and phase drift correction, the averaged dicom data was most universally supported by all the software packages evaluated, allowing a more encompassing cross-section of the methods.
All spectra used in this study were inspected for reconstruction issues, prior to their inclusion.

Spectral quality
Prior to quantification, the 2 ppm NAA and water peaks in each spectrum are fitted by Lorentzians, from which the position, full-width-halfmaximum (FWHM) and maximum signal amplitude of the resonances are obtained. The noise floor was estimated by computation of the standard deviation, σ, of the signal in the region 8−9 ppm, observed to be free of metabolite signal, then the signal-to-noise ratio (SNR) calculated: S 2σ . Water suppression factors are estimated by comparing the height and area of the water peak in the water unsuppressed and water-suppressed spectra.

Analysis methods
Analysis methods can be broadly separated into two categories: (i) basis set methods (TARQUIN 27 , jMRUI 22,26 , LCModel 24 ), which attempt to fit the data using a set of simulated or experimentally obtained basis spectra and (ii) peak fitting and integration tools (Gannet 21 , LWFIT -An in-house peak integration method used as a baseline in this study), which compute the area of peaks associated with various metabolites. As the absolute amplitudes are arbitrary and depend on the scanner design, transmitter calibration, coils used and the phantom itself. We focus on quantifying the relative contributions of metabolite signals, specifically the GABA-to-NAA ratio. The concentration ratio estimate, C, between a target metabolite, M , and reference, R, is given by: where S is the reported resonance amplitude, N is the number of contributing spins, R is the signal attenuation factor resulting from T1 and T2 relaxation, and E is the edit efficiency. For the purposes of our analysis, the difference between the relaxation properties of GABA and NAA are assumed negligible, as suggested by Choi et al. 29 . With this assumption, basis set methods report appropriately scaled metabolite amplitudes.
However, the peak fitting methods require corrections for N and E. Both Gannet and LWFIT quantify the 3 ppm GABA resonance and 2 ppm NAA singlet, resulting in N R N M = 3/2. By default, Gannet's edit efficiency correction factor is estimated as E R /E M = 0.5, and LWFIT's correction factor was calculated from simulation (Section 3.2). The concentration ratios reported by each tool were plotted vs. the known ground-truth ratio, and the relationship calculated via linear regression.
For the purposes of this study, the tools were employed as close to "out-of-the-box" as possible. Default settings were used with additional developer recommendations whenever these were readily available. While it is indeed possible to further optimise the parameters of each tool, our aim was to examine the potential methodological variation that could occur in GABA MRS studies. A brief description of each tool is provided below, and readers are directed to the referenced literature for further information.

Gannet
GABA-MRS analysis tool (Gannet) 21 is an open-source automated MATLAB-based analysis tool, specifically designed for automated processing of 3 T GABA MEGA-PRESS data. Gannet is able to process raw data, further reducing user dependence, performing channel combination, adding line broadening, frequency and phase corrections, outlier rejection and time averaging. Fitting is performed using non-linear least-squares algorithms (lsqcurvefit and nlinfi) with a five parameter Gaussian model to estimate the 3 ppm GABA signal in the difference spectrum, a six parameter Lorentzian model to estimate the 3 ppm Cr signal in the off spectrum and, if appropriate, a six parameter Gaussian-Lorentzian model to fit the unsuppressed water signal. Quantitative results are reported as integral ratios of GABA to Cr and a concentration relative to water, NAA or Glx (Glu and Gln combined). While Gannet attempts to remove user interaction from the analysis procedure, it allows user modification of the code and adjustment of the pre-initialisation script is encouraged depending on setup. For this paper Gannet version 3.0 was used and the pre-initialisation script was adjusted to account for room temperature phantom acquisition.

jMRUI
jMRUI 22,26 is proprietary software distributed under its own license terms and made freely available to registered users for non-commercial use.
It has a large suite of pre-processing, analysis and simulation options. For this study, two basis set methods were used to quantify the spectra.
The first, quantitation based on semi-parametric quantum estimation (QUEST) 48,49 , is a time-domain fitting tool, which includes a semi-parametric approach to handle spurious signals resulting from macro molecules and lipids. The second method, automated quantitation of short echo time MRS spectra (AQSES) 50 , uses a modified VARPRO variable projection NLLS algorithm 51 , which imposes prior knowledge in the form of upper and lower bounds on the non-linear parameters.
For this study, jMRUI version 6.0 beta was used. MEGA-PRESS basis sets were generated using NMR-SCOPE-B 52,53 by defining a single MEGA-PRESS pulse sequence with sequence parameters matching our experimental protocol. Two basis sets for GABA/Cr/NAA and GABA/Cr/-NAA/Glu/Gln were generated. Manual frequency and phase calibration and apodization were performed on the input spectra and fitting with both tools was applied using this common input. This allowed investigation of the influence of the fit algorithm itself, independent of pre-processing.

LCModel
LCModel 23,24 is a formerly commercial, Linux-based MRS analysis tool which recently moved to open-source. It can be used with its internal in-vitro basis set or any arbitrary basis set specified by the user. The baseline signal resulting from macromolecules is established using spline fits. Fitting is performed using a Gauss-Newton non-linear least-squares algorithm with Marquardt modification 54 with additional terms for T2, spectral shift and field inhomogeneities. For this paper, version 6.3-1L and the Purcell lab basis set 55 recommended by LCModel's creators were used.

LWFIT
Our in-house code 56 , LWFIT, written in MATLAB, aligns the edit-off and difference spectra so that the main NAA peak is precisely at 2.01 ppm and then calculates the areas for the NAA (1.91 ppm to 2.11 ppm, difference), Cr (2.90 ppm to 3.10 ppm, edit-off) and GABA (2.90 ppm to 3.12 ppm, difference) peaks, as well the two main peaks associated with Glu and Gln, GLX1 (2.25 ppm to 2.45 ppm, difference) and GLX2 (3.65 ppm to 3.85 ppm, difference), using the real frequency-domain spectra. This is done by numerical integration using piecewise-linear functions over the indicated fixed ppm ranges, selected to minimise contamination from other signals. Lorentzian, Gaussian and spline fitting and several baseline fitting and filtering methods were explored but initial testing indicated that numerical integration with minimal pre-processing yielded the most accurate estimation results of the methods considered. While more aggressive noise filters aesthetically improve the quality of the spectra, their application tends to exacerbate underestimation of GABA, especially for weak signals, suggesting that filtering eliminates some of the GABA signal present in the data 57 . For this work, only zero-filling (N = 4096 points) and a frequency shift to align the NAA peak in the spectrum to 2.01 ppm are performed.

Simulations
To further facilitate interpretation of the experimental results, quantum mechanical simulations were also performed. Simulation is fundamental in MRS quantification. Basis set fitting methods utilise simulated spectra to model the peak structure to which the experimental data are fit. While peak fitting methods are less tightly parametrised during the fit, they still must compensate for edit efficiency -the proportion of peak amplitude lost as a part of the editing process -in order to accurately quantify the metabolites. The simulations were used to investigate the influence of ideal RF-pulse simulation and GABA models on quantification, and to facilitate edit efficiency correction for LWFIT.
The experimental setup was emulated, modulating GABA concentration while maintaining a fixed concentration of other selected metabolites.
MEGA-PRESS simulations were performed using FID-A 61 , a simulation and data processing package for MRS. MEGA-PRESS difference spectra were generated, mirroring the sequence parameters of the phantom study with a finite-bandwidth editing pulse alternating between 1.9 ppm respectively. Simulated series were tested with all three GABA models to investigate any potential bias of our simulated data. For all simulated series, GABA concentrations were varied between 1 mM and 12 mM at 1 mM intervals and the combined spectra were analysed using our in-house code.

RESULTS AND ANALYSIS
We initially present the experimentally acquired spectra along with the corresponding calibration and quality measures. We then present the simulation and quantification results.

Calibration and characterisation of spectra
The stability of the experimental setup over the course of an entire series of experiments was assessed by considering the variation of core parameters such as transmitter voltage and frequency, shim gradients and water suppression settings. Table 2 shows that a small adjustment made to the experimental procedure for E3 and E4 to avoid table movements between successive scans significantly reduced adjustments of transmitter frequency, reference voltage and linear gradient offsets Gx, Gy, Gz, although all results are within acceptable limits.
The difference spectra for all four experimental series are shown in Figure 1. The expected pseudo-triplet structure, clearly visible for the wellcalibrated solutions E1 and E3, is still discernible in the gel phantom spectra E4, despite the broader linewidth. However, this structure is not present in E2 where we observe a broad, almost Gaussian resonance. Despite this, the highlighted area does appear to increase linearly with GABA concentration. Table 3 gives the linewidths of the water and NAA peaks obtained from Lorentzian peak fits, indicative values for the water suppression factors (WSF) and the SNR of the NAA peak. While the linewidths increase for the gel phantoms as expected, the width of the NAA peak in the difference spectrum remains around 4 Hz, suggesting good field homogeneity and frequency stability over the duration of the experiments. SNR over 200 for NAA in the solution series and over 100 for the gel phantoms are also well within the acceptable quality limits.

(a) Spectra for series E1 (b) Spectra for series E2
(c) Spectra for series E3 (d) Spectra for gel phantom series E4

FIGURE 1
Difference spectra aligned to 2.01 ppm NAA peak. The 3 ppm GABA peak is shaded according to LWFIT integration limits, and vertically offset spectra correspond to increasing concentrations of GABA for the phantoms.    Figure 2 shows the GABA-to-NAA amplitude ratio obtained by LWFIT vs. the actual imposed GABA-to-NAA ratio for three series of simulated experiments. Simulations using ideal excitation and refocusing pulses produce a gradient of 1.01 -within 1% of the actual ratio -an intercept of −0.05 and an R 2 value of 1.00 (up to 3 significant digits). However, for spectra simulated with finite-duration refocusing pulses with realistic pulse envelopes, the GABA-to-NAA ratio gradient is only 84% of its expected value for the NAA-GABA-Cre series and this is further reduced to 80% when Glu and Gln are added. Figure 3 shows the 3 ppm GABA peak for the edit-on, edit-off and difference spectra as a function of position on a 2D grid. The difference spectra near the boundary of the selected voxel indicate that the editing efficiency is significantly reduced. In an experimental setting, the editing efficiency will be further reduced by imperfect slice profiles of the excitation pulses and other factors such as field inhomogeneity. Assuming these issues lead to a further reduction of the editing efficiency by a factor of 0.84 and 0.80, respectively, the simulations  Edit efficiency is a crucial factor in edited MRS quantification. Basis set fitting methods should have this incorporated into the basis by appropriate simulation of the basis functions. Peak fitting methods, on the other hand, require an ad-hoc factor to compensate for the resulting reduction in peak amplitude. To ensure the best GABA estimates, full simulations with correct pulse envelopes are required to fully capture the J-modulated peak structure for basis set methods, and editing efficiency factor for peak fitting.

Quantification
The quantification results obtained by the various fitting approaches are summarised in Table 4. The results are reported as gradient, a, and intercept, b, with 95% confidence intervals and R 2 value of the linear fit. The concentration ratios are also plotted against the known ground-truth ratios in   Linear fit results for GABA-to-NAA ratios in all four experimental series. Reported are the gradient a and intercept b of the linear fit with 95% confidence intervals (CI) and the coefficient of determination R 2 (to two significant digits). Perfect quantification corresponds to a = 1, b = 0 and R 2 = 1. The percentage variation of reported amplitudes for NAA, whose actual concentration is constant, is also given. The large % variation of NAA amplitudes reported by TARQUIN for E1 is due to an unexplained rescaling of metabolite amplitudes between scans 10 and 11. The % variation of NAA for scans 1 -10 is 6.09%. There is no simple explanation for the large variation reported by TARQUIN for E4.  Figure 4b shows results for the intentionally miscalibrated solution series, E2. Unlike the other series, E2 did not exhibit the characteristic pseudo-triplet GABA resonance, and instead presented a single broad peak at 3 ppm. While a pH of 3.0 is beyond what is commonly observed invivo, this series was intended as a means to study the effect of chemical environmental changes upon quantification, a relevant consideration when investigating certain pathologies. Surprisingly, both TARQUIN and jMRUI give significantly more accurate GABA-to-NAA gradients, with a close to 1 and higher R 2 than for the pH-calibrated series. The TARQUIN change appears to be a consequence of its previous tendency to over-estimate GABA, while jMRUI's GABA gradient increase is likely driven by a reduced efficacy in NAA quantification. Peak fitting tools were expected to be more robust to environmental pH changes. With less tightly constrained models, the majority of signal is still captured by the tool, leaving only a systematic offset resulting from edit efficiency correction. We found that LWFIT and Gannet maintain a high R 2 , with LWFIT showing a reduced GABA-to-NAA gradient, likely a result of integration limits cutting off part of the broad resonance. Gannet, on the other hand, exhibits a modest increase, however there is significant overlap with the E1 confidence interval, so this finding may not be significant. Gannet's Gaussian model of GABA seems a suitable choice for this series, and the robustness of Gannet's estimations support this.  Table 2. Figure 4d shows the results for the gel phantom series with Glu and Gln. The fit quality for the gel series is generally lower, as expected due to higher FWHM and lower SNR for tissue-mimicking gels. The mean R 2 of E4 GABA-to-NAA gradient estimates across all tools is 0.65 for E4,

DISCUSSION
The main aim of this work was to develop an iterative phantom-based method for benchmarking GABA quantification software. The phantom data exhibit a high SNR with narrow metabolite line-widths, well above quality thresholds outlined in consensus work. The four experimental series provide unique challenges to quantification software, allowing elucidation of several confounds which would not be possible using in-vivo data, alone. E1 represents the simplest case with only GABA, NAA and Cr, E2 introduces environmental changes to pH, E3 adds Glx to investigate its potential to interfere in GABA quantification, and the gels of E4 present broader and lower SNR spectra. LWFIT -intended to serve as a baseline sanity check -reproduces the expected gradient to a high fidelity for E1, E3 and E4. LWFIT does, however, underestimate the GABA-to-NAA concentration gradient for E2, with its pH offset. The broad 3 ppm GABA resonance present in E2 is likely insufficiently captured by the integration limits used by LWFIT, furthermore, edit efficiency calculations are likely invalidated somewhat by the apparent structural changes in the spectra.
While all tools generally estimated a reasonable linearity in the GABA-to-NAA gradients, there is a systematic offset in these estimates. For TARQUIN, jMRUI and LCModel, these changes are primarily driven by the choice of basis set. This can be seen by comparing the two jMRUI methods -QUEST and AQSES -which share a common basis. Their mean estimates of GABA-to-NAA gradients differed by just 4% across all the experimental series. The discrepancy increases to around 106% when comparing TARQUIN and LCModel. The basis issue is further complicated by differences in the default GABA models used by the various tools. LCModel includes the Govindaraju and Kaiser model, while jMRUI uses its own modified version of the Govindaraju and Kaiser models. GABA spectra simulated using the Govindaraju, Kaiser and Near model, shown in Figure A2 exhibit clear differences in the structure of the 3 ppm GABA peak. Generally, the Kaiser and Near models appear to approximate the experimental spectra better than the older Govindaraju model, in agreement with observations in the literature. However, despite this, many tools still maintain the older Govidaraju GABA model by default. It is therefore not possible to fundamentally ascertain the performance of these tools without a systematic study with common basis sets and optimised settings. An interesting follow-up study would be to generate a common optimised basis set to be used in several basis set fitting methods. The phantom series can be further enriched by adding degrees of noise, linebroadening, experimentally acquired macromolecular baselines and other co-edited resonances to investigate their influence on quantification, providing many avenues of further study.
With the exception of TARQUIN, the tools present a tendency to underestimate GABA-to-NAA ratios to varying degrees. Simulations suggest realistic refocusing pulses reduce the editing efficiency and can lead to a reduction in reported GABA-to-NAA ratios of up to 20%. Including additional losses in experiments due to line broadening as a result of imperfect shimming, and excitation pulses with realistic slice profiles could therefore explain the underestimation of GABA-to-NAA gradients by many of the tools discussed. These losses can be taken into account by appropriate simulation during the creation of basis functions, or as a part of the edit efficiency correction for peak fitting tools. Table 5 contains a meta analysis of the reported GABA-to-NAA gradient values, specifically the mean gradient estimate and standard deviation of this value across all experimental series. The mean provides an indication of each tool's tendency to under or over estimate the changes in GABAto-NAA. Whereas the standard deviation, σ, provides an indication of the robustness of each tool's estimations to changes in the background conditions. A notable observation is that the peak fitting methods -Gannet and LWFIT -have far lower σ than all basis set methods. This suggests that the less strictly parametrised approaches might be preferable for edited spectroscopy analysis, where the landscape is greatly simplified compared to conventional MRS. This is in agreement with the observations of Marzola et al. 43 , who also reported the robustness of these methods in-vitro but note the poor performance of integration methods in-vivo. However, it should be noted that when series E2 is excluded, LCModel exhibits a similar results to Gannet, with a mean and standard deviation of 0.5583 and 0.0939, respectively.
The fit error is one aspect that was not approached in this study. While basis set methods generally report the Cramer-Rao lower-bound (CRLB) as the preferred measure of fit quality for a particular basis signal, this functionality is not present in Gannet and not valid for LWFIT, limiting comparability. Furthermore, the use of CRLB as a quality filter has been shown to cause bias in mean concentration estimates 66 . Fit residuals could provide some insight into fit accuracy, however these will be modulated by baseline approximation procedures, which differ substantially between tools, making comparisons potentially unfair. While efforts were made to ensure reasonable fit quality in this study, future comparative work will need to establish a means of common error reporting, using the same, sensibly justified metric. This is a fundamental consideration of quantification benchmarking.

CONCLUSION
The phantom benchmark dataset presented here provided a novel, interesting way to compare GABA MRS analysis software. The dataset is intended to form part of a larger systematic benchmarking of analysis tools, a necessary future step for the MRS community. While indeed preliminary, the initial comparison of software seems to indicate that a less-strictly parametrised peak fitting methods might be the best option for edited MRS, where the simplified spectral landscape does not necessitate complete modelling of the peak. This is, however, still a necessary step, as modelling of edit efficiency is important in reducing systematic offsets in this approach. The methodological dependence of the quantification results observed also suggests that care must be taken when comparing results across studies, where analysis pipelines may differ and under-reporting of fitting methods muddy the waters for study replication.