Adaptive Baseline Fitting for 1H MR Spectroscopy Analysis

Purpose Accurate baseline modeling is essential for reliable MRS analysis and interpretation — particularly at short echo-times, where enhanced metabolite information coincides with elevated baseline interference. The degree of baseline smoothness is a key analysis parameter for metabolite estimation, and in this study a new method is presented to estimate its optimal value. Methods An adaptive baseline fitting algorithm (ABfit) is described, incorporating a spline basis into a frequency-domain analysis model, with a penalty parameter to enforce baseline smoothness. A series of candidate analyses are performed over a range of smoothness penalties, as part of a four stage algorithm, and the Akaike information criterion is used to estimate the appropriate penalty. ABfit is applied to a set of simulated spectra with differing baseline features and experimentally acquired 2D MRSI — both at a field strength of 3 Tesla. Results Simulated analyses demonstrate metabolite errors result from two main sources: bias from an inflexible baseline (underfitting) and increased variance from an overly flexible baseline (over-fitting). In the case of an ideal flat baseline ABfit is shown to correctly estimate a highly rigid baseline, and for more realistic spectra a reasonable compromise between bias and variance is found. Analysis of experimentally acquired data demonstrates good agreement with known correlations between metabolite ratios and the contributing volumes of gray and white matter tissue. Conclusion ABfit has been shown to perform accurate baseline estimation and is suitable for fully-automated routine MRS analysis.

quire these signals to be manually appended to the metabolite basis.
Control over the level of baseline flexibility (or smoothness) is a common and necessary requirement of each of the baseline modeling methods outlined above. In spline based approaches, a combination of the number of spline functions for a given frequency range and the smoothness penalty parameter control the baseline flexibility. For LCModel, the frequency spacing between the spline basis functions is dependent on data quality, and is set at maximum of 1.5 times the estimated full width at half maximum (FWHM) of the metabolite resonances or 0.1 ppm [8]. Similarly, for the FITT algorithm a fixed Lowess filter smoothing value is used and wavelet coefficients with scales less than twice the FWHM are excluded from the baseline model to ensure smoothness [12]. In the time-domain truncation approach baseline flexibility is primarily determined by the number of initial data points to be omitted from the fit evaluation.
For QUEST and TARQUIN the number of truncated data points, and therefore degree of baseline flexibility, is set at a default value that may be adjusted by the user.
Automated methods to determine the correct degree of baseline flexibility are important for obtaining accurate metabolite levels independently of the analyst. Furthermore, the manual adjustment of baseline flexibility for each individual spectrum is impractical for MRSI studies -where hundreds of spectra may be acquired in a single scan.
Whilst LCModel provides automated adjustment of baseline flexibility, a growing number of analysts choose to manually override the default analysis settings by adjusting the spline spacing parameter (DKNTMN). The first reported use of this manual adjustment was to improve the modeling of macromolecular resonances in rat brain at 9.4 T [14].
More recently, this parameter has been adjusted to encourage flatter baselines [15,16,17], suggesting the default LCModel baseline flexibility may not be optimal in some cases.
Finding the optimal degree of baseline flexibility is a crucial question in MRS analysis research, yet few studies have investigated this topic in detail. Using simulated data, Ratiney et al. demonstrated how the interference between metabolite and baseline signals was reduced by increasing the number of omitted data points, but this came at the cost of inflating errors due to noise [18]. More recently, Near et al. showed how the estimated baseline in LCModel can depend strongly on spectral SNR and metabolite FWHM, and that errors caused by baseline instability may dominate over errors from spectral noise in some cases [19]. The influence of baseline flexibility has also been explored using experimentally acquired data, with a recent study demonstrating a 15% difference in metabolite levels when comparing between a default and less flexible baseline model [20]. Baseline flexibility has also been shown to have a strong influence of the measurement of 2-hydroxyglutarate [21].
In this study, we introduce a new method to automatically determine the optimal degree of baseline flexibility for a frequency-domain spline based fitting algorithm. Firstly, background is given on the use of penalized splines for optimal data smoothing. A fully-automated fitting algorithm is presented, incorporating a novel method to automatically estimate the optimal level of baseline flexibility. Finally, the new method is validated with simulated and experimentally acquired MRSI data.

| Penalized spline smoothing
Baseline signals have a characteristically smooth spectral appearance, and must be accurately modeled to avoid biasing metabolite estimates. In good quality 1 H MRS of brain tissue baseline signals have a low intensity, relative to the primary metabolite resonances, and are therefore challenging to estimate in the presence of noise. Estimating a smooth function from noisy data is known in statistics as "scatterplot smoothing", and a number of approaches have been developed [22]. In this section we briefly outline the method of penalized splines in the simpler context of scatterplot smoothing, before describing their use as part of an MRS fitting algorithm.
A spline is a piecewise function made up of one or more polynomial segments joined together at points known as "knots". A wide range of spline functions are possible, however Basis-splines, more commonly known as B-splines, are a popular choice for smoothing applications due to their favorable numerical properties [23]. The degree of a Bspline function determines its overall smoothness, and third degree (or cubic) B-splines are often used for smoothing.
A cubic B-spline basis is show in Figure 1A with an offset added in the y-axis for clarity. B-spline bases consist of regularly spaced overlapping spline functions, spanning the full range of x values.
A B-spline basis may be used to obtain a smooth estimate of a signal using simple linear regression. Figure 1B shows the result of a spline regression, where each spline function has been optimally weighted, such that the sum of the functions (spline fit line) is the least squares fit to the data. The desired smoothness of the fit is controlled by adjusting the spline knot spacing, which in turn changes the density of spline functions. In the case of Figure 1B more detail to be captured by the fit, however too many functions results in an increased sensitivity to noise and the smooth estimate begins to exhibit random fluctuations. Conversely, insufficient density of spline functions results in the spline fit being unable to model genuine smooth trends present in the data.
An alternative to adjusting the number of spline basis functions to achieve a desired level of smoothness is to introduce a penalty parameter into to spline regression model. The smooth estimateŷ, of our data y, is calculated as:ŷ = Bâ, where B is a B-spline basis in matrix form, andâ is a vector of corresponding spline weightings to be determined. In simple spline regressionâ is found by solving the normal equations to minimize the sum of the squared differences between y andŷ. In penalized spline regression, the minimization function Q B is adjusted to incorporate an additional term to enforce smoothness in the estimate: The λ parameter controls the degree of smoothness by penalizing solutions forâ that interact with the difference matrix D. The difference matrix may be constructed to penalize the first, second or higher orders of differences betweenâ values. Here, we exclusively use a second order difference matrix, which is particularly effective for modeling the smooth baseline features typically found in MRS data.
An example of the second order difference matrix is given in (2), which shows how increased differences between adjacent weighting factors inâ consequently increase the penalty term in Equation (1) and regressing the augmented y vector on the augmented B matrix to yieldâ.
The general approach for penalized spline regression is to over-specify the number of B-spline basis functions, and primarily control the smoothness through the adjustment of λ acting on a difference matrix. We refer to this approach as "P-splines", first introduced by Eilers and Marx [24]. Whilst the value of λ has a direct effect on smoothness, it is an unintuitive parameter to interpret, as the optimal value often varies by several orders of magnitude. In addition, λ has a complex dependence on the density of the spline functions, the number of data points and other factors unrelated to smoothness. A more intuitive measure of the smoothness of a P-spline model is known as the effective dimension (ED), proposed by Hastie and Tibshirani [25]. For a given value of λ, B-spline basis B and difference matrix D, we calculate ED as follows: where tr denotes the trace of a matrix. smaller value of ED corresponds to a smoother baseline estimate.
Using simulated data we investigate the relationship between λ, ED and the optimal level of smoothness. The top left panel in Figure 2A shows a simple sine function with added normally distributed noise, shown in red and black respectively. Candidate P-spline smoothers, with differing values of λ, are shown in the remaining 5 panels. Since the true shape of the underlying smooth function is known, the error of the P-spline estimate may be calculated as the sum of squared differences between the true function and the estimate. A plot of the error as a function of λ is shown in Figure 2B. The sum of squared differences between the P-spline estimate and noisy data (residual) and the ED are also shown as a function of λ in part B.
Inspection of the error plot reveals the optimal λ to be approximately 20 -corresponding to an ED value of 12, and this can be intuitively verified from part A. For larger values of λ, the estimate approaches a straight line fitfailing to capture the details in the sine cure and resulting in an increasing fit residual and fit error (part B). Models with insufficient freedom to adapt to genuine trends in the data result in biased estimates and this is known as "underfitting".
Conversely, too much freedom in the smoothing model (small λ) results in the estimate becoming overly sensitive to random fluctuations in the data resulting in "overfitting". In least-squares fitting, there is a temptation to equate the smallest residual with the optimal fit, however Figure 2B clearly shows an increase in the error for λ values of less than 20 -despite a steady reduction in the residual.
Determining the optimal value for the smoothness parameter is one of the primary challenges for P-spline smoothing. A careful balance between instability from overfitting, and bias from underfitting, must be made for the most accurate estimate, and searching for the smallest fit residual is a useful, but insufficient metric of quality. Numerous approaches have been proposed to find the optimal smoothing value [22], and in this study we use Akaike's information criterion (AIC) [26]: where ln denotes the natural logarithm and n is the number of data points. The AIC is typically used to compare models, with lower values indicating an improved balance between bias and instability. A commonly used alternative is known as the Bayesian information criterion (BIC), however the AIC is chosen here for its simpler form, which may be intuativly modified -as shown in the following section. A plot of the AIC as a function of λ is shown in the lower panel of Figure 2B, and the λ value with the lowest AIC shows good agreement with the true minimum error.

| MRS baseline estimation using P-splines
In 1 H MRS analysis, the baseline signal must be estimated in the presence of numerous overlapping metabolite lipid and macromolecule signals. Fortunately, these non-baseline signals have a known molecular origin and are therefore well characterized and accurately simulated from established parameters [27]. We can update Equation 3 to incorporate these additional molecular components by arranging into columns of the basis matrix M, and appending to the P-spline basis B: As with Equation 3, we regress the basis matrix on y to yieldâ, however since metabolite signal amplitudes are always positive, analysis stability may be improved by enforcing a non-negative constraint on the subset ofâ values corresponding to the weightings on the basis set M (â M ≥ 0). The active-set method developed by Lawson and Hanson [28] is used herein to find the least-squares solution under the non-negative constraint -while still allowing theâ values corresponding to the spline basis (â B ) to remain unconstrained. Note that starting values are not required for any of the metabolite, macromolecular, lipid or spline basis function amplitudes when using this method.
A simulation study was performed to investigate the relationship between baseline smoothness, the accuracy of metabolite estimates and the AIC. Metabolite signals were simulated from known parameters [27] at levels consistent with normal brain tissue [29] -listed in Supporting Information Table S1. An experimentally derived macromolecule profile was also included in the simulation and basis matrix to yield a realistic spectrum [30]. This profile was comprised of multiple broad resonances combined with fixed relative intensities, and corresponded to only one element in the basis set. Simulated experimental conditions consisted of a field strength of 3 T, and semi-LASER localization (TE=28 ms). 1024 complex points were generated at a sampling frequency of 2000 Hz, 6 Hz Gaussian line-broadening was applied prior to zero-filling to 2048 points and Fourier transform to the frequency-domain. Ideal pulses were assumed and relaxation effects were not considered. 8

Wilson
Ideal spectra were distorted with normally distributed noise resulting in a spectral SNR of 54 -typical for 1 H MRS acquired from the human brain at 3 T. Signal strength was measured as the highest data point from the nominal NAA singlet peak at 2.01 ppm. A broad Gaussian resonance was added at 1.3 ppm with a linewidth of 100 Hz to simulate baseline distortion originating from scalp lipids. The degree of P-spline baseline flexibility is defined as the baseline ED (effective dimension) per ppm, which is more easily compared across analyses due to its independence from the number of points in the fit and the number of spline basis functions used. For example, an ED per ppm of 5 corresponds to an ED value of 19 when fitting the spectral region between 4 and 0.2 ppm (5 × (4 − 0.2)). Metabolite level estimates (â) were calculated using Equation 7 from real-valued data points in spectral region between 0.2 and 4 ppm, over a range of 10 levels of baseline flexibility. 32 spectra were analyzed at each level of baseline flexibility, with only the noise samples -drawn from the same distribution -randomly differing between each spectrum. The metabolite estimation error was calculated from the sum of the squared differences between the true amplitudes, listed in Table S1 (a), and estimated values (â) for each spectrum. Errors in the estimation of macromolecular or lipid signal amplitudes in the basis set were not considered. The mean and standard deviation of metabolite estimation errors was calculated across the 32 spectra at each level of baseline flexibility to evaluate accuracy and consistency.
Error values are provided in absolute units as further scaling was not performed.
Metabolite estimation errors are shown in Figure 3A, displaying a comparable shape to the simpler model in Figure 2B. Over the range of baseline flexibility studied, underfitting with an inflexible baseline results in much greater metabolite estimate errors compared to overfitting -which can be verified by inspecting the fit result plots in Figure   3 parts C-F. Good agreement is also seen between the metabolite estimate error and the AIC curve (part B), with a low AIC value corresponding to the most accurate level of baseline flexibility.
From the results presented in following sections, it was found that the AIC had a tendency to slightly overestimate the optimal level of baseline flexibility, and that a simple modification to Equation 6 improved accuracy: Setting m to a value of 5 was empirically found to be a good compromise between bias and variance for all simulated and in vivo analyses presented herein. The value was kept constant for all analyses presented, placing a greater penalty on overly flexible baseline models -resulting in a smoother baseline estimate relative to the standard AIC ( Figure 3B).

| Adaptive Baseline Fitting Algorithm
In this section we describe a fully-automated 1 H MRS analysis method based on the P-spline fitting approach presented in the Theory section. The emphasis of the design is to automatically adapt the baseline flexibility for improved accuracy, and we refer to the full algorithm as: Adaptive Baseline fitting or ABfit. The algorithm consists of 4 main steps, which will be described in order of execution: 1. Coarse frequency alignment.

| Approximate iterative fitting (step 2)
The goal of this stage of the algorithm is to find estimates of the three parameters with the largest influence on the fit residual: 1) the signal phase φ 0 ; 2) an approximate lineshape parameter d g and 3) an improved estimate of f o . The phase and frequency offset are applied to the acquired data as follows: The lineshape parameter applies Gaussian line-broadening to each signal in the basis set M, and the scaling factor β is introduced to define d g as the Full Width at Half Maximum (FWHM) measured in Hz: where denotes element-wise multiplication of the line-broadening function for each column of the basis matrix.
A simple, one parameter, lineshape model was chosen for this step to reduce the chance of poor solutions being found by the optimizer due to presence of local minima. Gaussian lineshape broadening was applied to all signals in the basis set M, including macromolecues. Metabolite signals in this paper were simulated with pure Lorentzian lineshape to model T2 relaxation -resulting in a Voigt lineshape following the Gaussian broadening applied in this step [31]. Combining the modified basis with the P-spline basis, and solving forâ in the least-squares sense with the same constraints as the previous section: Note Nelder-Mead simplex algorithm with bound constraints [32], minimizing: A λ value corresponding to 1 ED per ppm is used with a P-spline basis containing 15 spline functions per ppm, resulting in a total of 33 (15 × (4 − 1.8)) components in the basis B, with the same density of spline functions being used for all subsequent steps. A value of 1 ED per ppm was empirically found to work well with all data presented herein, however may need to be adjusted for spectra with fundamentally different levels of baseline complexity -such as 31 P MRS. Constraints are placed on the optimization algorithm to restrict the linebroading parameter to be positive and produce additional broadening of less than 15 Hz FWHM. The frequency offset is also constrained to be ±10 Hz different from the coarse frequency alignment (step 1).

| Baseline smoothness estimation (step 3)
Following the determination of the frequency offset, spectral phase and approximate lineshape the next step of ABfit is to estimate the optimal value for λ. A set of candidate fits are automtically performed with differing level of baseline smoothness, whilst maintaining the three main fitting parameters constant at their optimized values -as determined in the previous step. The maximum candidate baseline flexibility is set to a λ value corresponding to 7 ED per ppm, and the minimum value set to an ED of 2 -equivalent to straight line fit. 20 candidate fits are performed with logarithmic intervals between the maximum and minimum values of ED per ppm, where λ is back calculated from Equations 4 and 5, and the mAIC is calculated for each fit (Equation 8). The optimal λ value, and therefore ED per ppm, is found by selecting the candidate fit with the lowest mAIC. In contrast to the previous step, a wider spectral range of 0.2 to 4 ppm is analyzed to include the majority of conventionally measured metabolites -corresponding to 57 spline functions.

| Detailed iterative fitting (step 4)
In the final stage of ABfit, the overall lineshape model is refined alongside f o , φ 0 and minor adjustments to the frequencies and linewidths of the known basis signals. An asymmetric lineshape model is generated in the frequency-domain by smoothly adjusting the Gaussian linewidth parameter as a function of frequency -as described by Stancik and Brauns [33]: Note, the frequency dependence on the linewidth function is eliminated when the asymmetry parameter a g = 0, leading to s = d g and pure Gaussian broadening. l FD is transformed to the time-domain and applied to each of the basis signals in M.
Minor individual changes in the frequency offset and linewidth are applied to each molecular signal in the basis to accommodate variations in the acquired data:  Simulated spectra and metabolite estimation errors were calculated as described in section MRS baseline estimation using P-splines unless stated otherwise.

| Validation with simulated data
An experimentally derived macromolecule profile, comprised of multiple broad resonances combined with fixed relative intensities, was included in the simulations to yield a realistic spectrum [30]. In the first simulation study a broad simulated peak with a Gaussian lineshape (FWHM 100 Hz) at 1.3 ppm was combined with metabolites and the macromolecular signal with a spectral SNR of 54 (Figure 4 and Supporting Information Figure S1). The second simulation study was identical to the first, however the broad simulated peak was removed to yield a perfectly flat baseline (Figures 5 and S2). The third simulation study was identical to the first, however the amplitude of the simulated broad resonance was increased by a factor of 2 (Figures 6 and S3). The data for the fourth simulation study was identical to the second, with a perfectly flat baseline and an experimentally derived macromolecular signal (Figures 7 and S4). However, the experimentally derived macromolecular signal was removed from the fitting basis set and replaced with the set of simulated lipid and macromolecular signals commonly used by default in the LCModel and TARQUIN algorithms (listed in Table S2). The fifth simulation study was identical to the fourth, with the exception of an additional broad simulated peak with a Gaussian lineshape (FWHM 100 Hz) at 1.3 ppm (Figures 8 and S5). The final simulation study was identical to the first, however 12 Hz linebroading was applied to each signal in the basis set (rather than 6 Hz) -resulting in a spectral SNR of 34 (Figures 9 and S6).

| Validation with experimentally acquired data
Whilst simulation studies are important to assess the true accuracy of a method, it is challenging to adequately model the true range of variation present in MRS data. Therefore, ABfit was tested on experimentally acquired MRSI data to ensure validity and robustness to common artifacts -such as baseline distortions from scalp lipids, shimming variations and minor shifts in metabolite frequency.
MR data were acquired from two healthy adults with a 3 T Siemens Magnetom Prisma (Siemens Healthcare, Erlangen, Germany) system using a 32-channel receiver head coil-array. T1-weighted MRI was acquired with a 3D- Following acquisition, the four corner voxels were excluded from the central 8 by 8 grid due to their close proximity to the diagonal saturation regions, and ABfit was performed on the remaining 60 voxels. ABfit was applied without any manual adjustments, and exactly the same algorithm was used to analyze the simulated and acquired data. The percentage of white matter, gray matter and CSF contribution to each voxel was measured from segmentation of the T1 MRI using the FAST method [36] as implemented in the FSL software package (v6.0.1). The gray matter fraction was calculated as the percentage volume of gray matter divided by the sum of gray and white matter volumes, and compared with metabolite ratios.
The acquisition of human data included in this study was conducted with the approval of an Institutional Ethics Board.

| Simulation
The first simulation test contained a baseline distortion to mimic a spurious signal from scalp lipids, and the results of ABfit analyses are shown in Figures 4 and S1. The metabolite estimation error plot as a function of baseline flexibility shows the same characteristics as the simpler analysis model used in Figure 3A, however the errors are to the automatically determined baseline flexibility is shown in part C, where the true shape of a broad Gaussian peak centered at 1.3 ppm is apparent from the baseline estimate. A comparison between the true and estimated baseline components is shown in Figure S1 part C with the largest errors shown around 1.3 and 4 ppm. These spectral frequencies correspond to the lactate and alanine resonances which have been overestimated (part B) due to their overlap with the strong baseline artifact and noise.
The baseline distortion was removed for the second simulation study to test the ABfit approach for ideal spectra where the basis set alone is sufficient for accurate analysis. The error plot is shown in Figure 5, illustrating the absence of bias due to baseline underfitting. A compromise between bias and variance is unnecessary in the ideal case, as the correctly determined baseline flexibility has the lowest error and variability. Compared to the first simulation study, the metabolite errors are much smaller and may be attributed to noise ( Figure S1 vs S2 part B).
In the third simulation test the amplitude of the broad distortion at 1.3 ppm was doubled compared to the first.
A reasonable estimate of the optimal baseline flexibility is found using the ABfit method ( Figure 6A) with comparable levels of accuracy relative to the reduced amplitude baseline distortion ( Figure 4A). Figure S3 parts B and C show the largest metabolite errors for Lactate and Alanine due to their strong overlap with the baseline distortion at 1.3 ppm.
In the fourth simulation test a set of independent broad lipid and macromolecular signals in the basis were used to model the experimentally derived macromolecular profile. Figure 7A shows the lowest mean errors are obtained at higher levels of baseline flexibility (part D) -indicating a stronger bias at lower level of flexibility (parts B, C) compared to the previous simulations. This is likely explained by inadequate modeling of the broad macromolecular components around 3.8 ppm since these are not present in the commonly used individual macromolecular and lipid basis (Table   S2). Figure S4 part D shows a much larger discrepancy between the experimentally derived macromolecular profile used to simulate the data "true" and the estimate "est." compared to the previous simulation results. A combination of metabolites (part B) and baseline contributions (part C) are used to model the discrepancy in the macromolecular signal components -resulting in increased metabolite errors.  The fifth simulation test was the same as the fourth, except a broad baseline distortion was added at 1.3 ppm. In the final simulation test, the first test was repeated with broader metabolite signals to evaluate the efficacy of the automated baseline determination for poorly shimmed data. Errors are larger compared to the first test, due to the reduction in SNR and spectral resolution, however a reasonable estimate of the optimal baseline flexibility of approximately 4 ED per ppm is found ( Figure 9A). The primary source of metabolite error is shown in Figure S6 part B, where the lactate and alanine signals are confused with the baseline distortion to a greater degree compared to the first simulation test.

| Experimental
Significant scalp lipid contamination was found in some of the voxels from the first MRSI scan, most likely resulting from subject movement between the T1 anatomical and MRSI acquisition. ABfit result plots are shown in Figure 10   High quality MRSI was acquired for the second scan with only minimal baseline distortion observed across all analyzed spectra. A strong correlation between metabolite levels and the underlying tissue contribution is visually apparent from Figure S7 part A, with an increased tNAA / tCr ratio in white matter compared to gray matter (tNAA = NAA + NAAG, tCr = PCr + Cr, tCho = GPC + PC, Glx = Glu + Gln). A linear regression of select metabolite ratios with the gray matter fraction is plotted in Figure S7 parts B, C and D, with strong correlations observed -in good agreement with high field observations [37,38]. The mean full-width at half maximum resolution across the voxels analyzed was 0.032 ppm (3.9 Hz), measured from the tNAA resonance. The mean SNR was 85, with the noise region defined as the real valued data points between -0.5 and -2.5 ppm.
Two example fits from one of the voxels in the second MRSI data set are shown Figure S8. In part A, the ABfit method is applied as described previously, and in part B the lineshape asymmetry parameter (a g ) is heavily constrained to enforce a symmetric lineshape model. A smaller fit residual in the tNAA spectral region is found for the asymmetric lineshape model, justifying the minor increase in modeling complexity associated with an additional fit parameter.

| DISCUSSION AND CONCLUSIONS
A new algorithm to automatically determine the optimal level of baseline flexibility has been developed and validated using simulated and acquired MRS data. LCModel is currently the most widely used approach for automatically estimating the optimal baseline flexibility, where the optimal penalty parameter (α B ) is chosen by gradually increasing its value until the boundary of the 50% confidence region for the fit is achieved -estimated by comparing successive fits to the first in the series [39,8]. In contrast to LCModel, the method presented here uses a modification to the AIC to determine the optimal penalty factor -as part of a four-step fitting procedure. An additional difference is the use of P-splines in ABfit compared to smoothing spline approach [40] used in LCModel.
Sima and Van Huffel proposed the use of the classical generalized cross-validation (GCV) criterion, combined with a golden-section search, to determine the optimal penalty parameter value [41]. High accuracy was demonstrated for simulated data, however the method was only tested with good starting values for the non-linear fitting parameters, which are not typically available for experimentally acquired MRS. In ABfit, these non-linear parameters are estimated using a simplified initial fit (step 2) and subsequently refined (step 4). More recently, Zhang and Shen showed that a measure of the baseline uncertainty is also a useful criterion to determine baseline smoothness for simulated and experimentally acquired data [42].
ABfit was shown to find a reasonable compromise between bias and variance for the majority of simulation tests.
However, in the fourth simulation test, where a set of independent approximate macromolecular signals were used to fit a realistic macromolecular model, the optimal flexibility was less easily determined. In this case, the AIC penalty for an overly rigid baseline was minor, since a low residual could still be achieved through the increased freedom afforded by using a set of independent macromolecular signals ( Figure 7, part B). However, in this case, a more rigid baseline introduces a greater level of metabolite estimation bias due to the mismatch between the simulated and modeled macromolcular profiles. This represents a significant challenge for automated baseline selection for data with lower SNR -a problem previously identified by Near et al [19] using the LCModel package.
Potential solutions to reduce metabolite estimation bias associated with low SNR include the use of more accurate macromolecular modeling [30] and opting to use a fixed level of baseline flexibility. Whist the default approach for ABfit is to automatically select the level of baseline flexibility, fitting options are implemented to specify a fixed degree. Alternatively, baseline flexibility may be systematically adjusted by changing the mAIC scaling parameter m, and this may be necessary for data with significantly different baseline characteristics such as 31 P MRS. Each of these approaches has advantages and weaknesses and should be justified depending on the study aims. For example, in functional-MRS a small change in metabolite levels is generally sought, and therefore a less flexible baseline may be preferred -since any metabolite estimation bias is eliminated when the change is normalized to a well determined signal. As a general rule, the baseline should be sufficiently flexible to eliminate any broad features present in the residual. However, poor spectral SNR will mask these features, and therefore caution is advised when comparing metabolite levels between spectra with greatly differing SNR -particularly when macromolecular signals are only partially modeled by the basis set.
In conclusion, new MRS analysis method with adaptive baseline modeling is presented and validated on simulated and experimentally acquired data. The approach is fully-automated and integrated into a free and open-source software package -providing a transparent and reproducible platform for future MRS studies [43].