Introduction

Magnetic resonance imaging and spectroscopy is a powerful tool for non-invasively assessing anatomic and metabolic changes that occur in brain diseases. In clinical spectroscopic studies and especially for magnetic resonance spectroscopic imaging (MRSI) data acquisition, short repetition times (TR) are often required to meet scan time constraints, but accurate metabolite longitudinal relaxation time values (T 1) are then needed to correct the metabolite concentrations for the T 1-weighted effect. The metabolite T 1 values are likely to be important for quantifying results that make comparisons between patients and normal controls. Moreover, the knowledge of 1H metabolite longitudinal relaxation times can by itself give insight into the properties of a given region of interest.

In many previous studies, the estimations of the metabolite T 1s were performed using single voxel acquisitions. Short echo time spectra coming from either progressive saturation [13] or inversion recovery experiments [46] were collected and T 1 values were usually derived from amono-exponential fit. The inversion recovery experiments typically used long repetition times (TRs are usually equal to 6 s) with varying inversion times [46], which is prohibitively long for MRSI experiments. These single voxel approaches assume a single T 1 over the whole voxel regardless of its tissue composition. To obtain white matter or gray matter T 1 values and good SNR, large (usually greater than or equal to 8 cc), and relatively heterogeneous single voxels were typically acquired. In most cases, gray matter (GM)-T 1 results were obtained from voxels containing 60–70% of GM, while white matter were obtained from voxels containing 70–90% of WM. At the same time, different MRSI studies [710] using linear regression demonstrated how metabolite concentrations (and thus metabolite signal intensities) can be different according to their tissue origin. Thus a common concern about the single voxel studies is whether the content of GM and WM in the examined voxel has an influence on the metabolite T 1 results. Moreover attempting to reduce the size of the single spectroscopic voxel to reduce the voxel tissue heterogeneity would result in increasing the number of averages and the scan time. In contrast, MRSI techniques offer the possibility to acquire simultaneously several spectra over a wide brain region and at a resolution allowing tissue analysis. The first goal of this paper was, therefore, to develop a MRSI post-processing method to estimatemetabolite T 1s while accounting for the voxels’ tissue content. To date, no published studies investigated the use of MRSI data to estimate metabolite T 1 relaxation times.

Then, as for any quantitative measurement based on model fitting, an assessment of the precision of the T 1 estimation is desirable. A benefit of MRSI is that it provides several spectra and thus several data points for the metabolite T 1 fit which can be resampled in a bootstrap manner to estimate standard error. Therefore, a second goal of this study was to develop an approach to obtain metabolite T 1 standard error by bootstrapping.

The last contribution of the paper was to apply the new techniques to measure metabolite T 1s in different regions and tissues of the brains of healthy subjects.

The proposed method estimates metabolite T 1 relaxation times by using 2D MRSI data at different repetition times. The progressive saturation method was chosen instead of inversion recovery method for scan time concern. While conventional techniques spend time in averaging single voxel acquisitions to obtain good SNR, we use this time to acquire multi-voxel data and investigate regional and tissue specific metabolite T 1 differences. The post-processing takes advantage of the combination of segmented MRI and spatially distributed spectroscopic data to investigate either WM versus GM metabolite T 1 values and/or regional differences in longitudinal relaxation times. The proposed method relies on three major concepts:

  1. 1.

    Increasing the SNR by averaging voxels according to their WM/GM content and their location (for example anterior vs. posterior) since the SNR of the metabolite signals is very low using the proposed acquisition parameters (number of excitations (NEX) and number of voxels in the slice) at 1.5 T.

  2. 2.

    Estimating metabolite T 1 for gray and white matter using a non-linear least squares algorithm. The underlying model function used in the fitting procedure associates WM/GM content of a voxel to the metabolite signal intensity.

  3. 3.

    Using a bootstrap technique to assess uncertainty on the metabolite T 1s and taking into account this confidence when calculating a group mean value.

This paper presents the techniques developed to utilize MRSI data for metabolite T 1 measurement. A validation of these techniques is then proposed through Monte Carlo data simulations, demonstrating the statistical performance of the proposed method. Finally the fitting procedure is applied to 2D conventional MRSI data acquired at 1.5 T from eight healthy subjects.

MR data acquisitions

Study subjects

A total of eight healthy volunteers (five females and three males, mean age 31.5±9.5 years) were examined to validate our method. Written and informed consentwas obtained from all participating subjects. The study was approved by the UCSF Committee on Human Research.

MR parameters

The healthy volunteers were scanned on a Signa 1.5 T clinical imager from GE Medical Systems (GE Healthcare Technologies, Waukesha, WI, USA) using a quadrature head coil. Short echo time (TE = 35 ms) 2D MRSI data sets (12×12, 1 cc resolution) were acquired using a PRESS volume selection at five different TRs (TE = 0.850, 1, 2, 4, 8s).The number of excitations (NEX) acquired were as follows: NEX=3 for TE = 0.850s, NEX=2 for TR = 1s, NEX=1 for TR = 2, 4 and 8s. Oblique Fast Spin Echo images were used to guide the positioning of the spectroscopic acquisition. Care was taken that the PRESS box avoided the ventricles and was centered in the anterior-posterior middle of the corpus callosum body. T 1-weighted 3D spoiled gradient recalled (SPGR) images were also acquired for segmentation of the anatomic images. The setup and data acquisition time for the anatomic and spectroscopic imaging was approximately 55 min.

Methods

The first goal of the analysis was the formulation of an approach to estimate regional metabolite T 1s in cortical gray matter and white matter.

The model function

With the assumption made by several previous MRSI studies [79] of a linear relationship between the voxel WM/GM content and the metabolite peak signal intensity, we applied the following simplified model for the measured signal intensity S (for example S =signal amplitude of NAA) in the nth MRSI voxel

$${S_n}({\rm{TR}}) = p_n^{{\rm{WM}}}S_0^{{\rm{WM}}}f({\rm{TR, }}T_1^{{\rm{WM}}}) + p_n^{{\rm{GM}}}S_0^{{\rm{GM}}}f({\rm{TR, }}T_1^{{\rm{GM}}})$$
(1)

where n runs through all the voxels in a given brain region (for example anterior and posterior). p WM n and p GM n are the calculated fraction of WM and GM, respectively within the nth spectroscopic voxel, S WM0 and S GM0 are respectively the signal intensity of a fully relaxed (TR ≫ 10s) resonance pertaining to assumed pure WM and pure GM, f is a function that characterizes the T 1-weighting at a certain TR. We used the usual mono-exponential function [2], \(f({\rm{TR}},{T_1}) = 1 - {e^{ - {\rm{TR}}/{T_1}}}\). S denotes either a resonance intensity corresponding to only a part of a metabolite or a whole metabolite signal intensity depending on whether T 1 values are assumed to differ between different parts of the molecule. We then, for example, split the creatine signal into two parts, assuming that the CH3 singlet at 3.02 ppm has a different T 1 than the CH2 singlet at 3.91 ppm. Note that this model does not necessarily require the WM and GM T 1s to be distinct, but rather relaxes the constraint of having only one T 1 independent of the tissue type.

The T 1 estimation procedure

Metabolite T 1s are obtained after the following steps:

Step 1. Calculation of p WM n and p GM n

To compute the fraction of WM and GM, p WM n and p GM n of Eq. 1, for each spectroscopic voxel, the anatomic T 1-weighted images (SPGR) were segmented into GM, WM and cerebrospinal fluid (CSF) compartments using the Automated Segmentation Tool, FAST [11]. Then these segmentation masks were resampled using nearest neighbor interpolation into the coordinate system of the T 2-weighted image used for the prescription of the spectroscopic grid. The transformation used was calculated from the position and orientation of the T 1 and T 2 images, assuming there was no motion between the scans. Finally, for each voxel, the fraction of WM and GM within a voxel was computed from these segmentation masks by counting the number of WM or GM pixels within a spectroscopic voxel. Thus, the WM and GM proportion maps are at the spectroscopic resolution as shown in Fig. 1. The chemical shift displacement was minimized by exciting a larger region than desired (using a larger PRESS selected volume than usual) and then using very selective saturation (VSS) pulses to eliminate extra signals and to obtain the final, desired selected region [12].

Fig. 1
figure 1

a Segmentation of the white matter, b Percentage of white matter at the spectroscopic resolution, c Segmentation of the gray matter, d Percentage of gray matter at the spectroscopic resolution, e Example of original spectra coming from the four highlighted voxels TR = 1s, NEX = 2, SNRNAA ≈ 2

Step 2. Voxel averaging

At 1.5 T, the SNR (defined for example as the ratio of NAA time-domain singlet amplitude to the standard deviation of the time domain noise) of our acquisition was very low, typically between 1.5 and 2.8. Consequently the estimation of the metabolite intensity S is not accurate. For each TR, we proposed to generate N gen new signals from the N total original signals of the MRSI grid by averaging N avg chosen signals to increase the SNR. The N avg signals are randomly selected among the nearest neighbors of the voxel in terms of WM/GM content to track the original WM/GM distribution while reaching a greater SNR. For each averaged signal k, the new fractions p WM k and p GM k are also calculated by averaging the fractions corresponding to the N avg chosen signals, see Fig. 2. It is possible, with this technique to have more voxels than originally, N genN total as we can draw several sets of N avg voxels for one original voxel. Figure 2b, c show the SNR gain between the original in vivo spectra and the averaged spectra.

Fig. 2
figure 2

a Graph showing the WM-GM percentage distribution of the original Ntotal voxels (black) as well as the Ngen averaged voxels. b Example of four spectra coming from the original voxels (Navg = 0, SNRNAA ≈ 2). c Example of four spectra coming from the averaged voxels (Navg = 6, SNRNAA ≈ 4.4)

Note that to perform regional analysis, this averaging procedure is applied on a specific part of the spectroscopic grid (for example in the anterior part of the slice) by working only with the voxels belonging to the region of interest.

Step 3. Quantification of S

The signal intensity S of Eq. 1 has to be estimated for each metabolite of interest, at each TR value, and for n running through all the N gen averaged voxels. A quantification or a peak picking method can be used as is done for single voxel studies. We applied the time-domain algorithm QUEST [13] using a simulated basis set suitable to the acquisition. In a preprocessing step, the residual water was eliminated using the HLSVD method [14]. At 1.5 T, six metabolite patterns were estimated: the whole metabolite pattern for N-Acetyl compounds (NAAt = NAA + NAAG, mainly located at 2.02 ppm that we will call summarily NAA), choline compounds (3.21 ppm), glutamate (Glu, ∼2.3 ppm), and myo-Inositol (mI, ∼ 3.6 ppm), and the two patterns for creatine (Cr-CH3, 3.03 ppm, Cr-CH2, 3.91 ppm). The metabolites were simulated using the NMR-SCOPE[15] program of the jMRUI software package [16]. The background contamination coming from broad macromolecules, lipids and surrounding broad pattern of metabolites of low concentration were automatically modeled and taken into account in the semi-parametric approach utilized in QUEST [17]. The T 1s of the NAA singlet (T 1,NAA), Cr-CH3 (T 1,Cr-CH3), Cho (T 1,Cho) and mI (T 1,mI) were investigated. Note that the T 1 of Cr-CH2 was not estimated because of its poor signal quality due to its proximity to the water in the spectrum. To ensure the quality of the four-parameter fit, good quantification results were selected prior to step 4. Voxel results were selected following these two criteria:

  1. 1.

    estimated relative Cramér-Rao lower bound [18] of the metabolite amplitude, (rCRB) < 15%.

  2. 2.

    metabolite peaks that had smaller than 9Hz linewidth at half peak height.

After this selection, we have N final,TR points at each TR for the metabolite T 1 fit.

Step 4. Four-parameter fit

The four parameters S WMO , S GM0 , T WM1 , T GM1 of Eq. 1 were fitted using a non-linear least squares algorithm. We used the lsqnonlin method from the Optimization Toolbox of Matlab (The MathWorks, Natick, MA, USA). See Fig. 3. If there are N final,TR selected results and there is a number N TR of repetition times, one has \(\Sigma^{N^{TR}}\) N final,TR (usually 100–300) data points to fit the four parameters.

Fig. 3
figure 3

Map showing an in vivo example of the four parameter fit for NAA compounds. In this case \(\Sigma^{N^{TR}}\) N final,TR = 480 data points (5 TR, 96 points per TR) were available for the four parameters fit

Estimating metabolite T 1 uncertainties using bootstrap

Bootstrap is an empirical, non-parametric statistical technique based on data resampling. It is used to make statistical inference such as the variance estimate of some fitted parameters. Although well-known and used in other MR modalities such as fMRI [19] or diffusion tensor MRI [20], it has rarely been applied, to our knowledge, to MRS parameters [21]. This computer-based method relies on the drawing of some bootstrap samples [22]. In our case, a bootstrap sample consists, for each TR, in a random sample of size N final,TR, say S j/*(TR), 1 ≤ jN final,TR, where S j/*(TR) is drawn with replacement from the original data points S j (TR). To avoid some downward bias on the estimated standard error, we applied a bootknife approach [23] which is a resampling technique combining the features of jackknife and bootstrap. For each TR, one point is first randomly omitted from the original sample of size N final,TR (jacknife), then, from the remaining sample with size N final−1, a bootstrap sample of size N final,TR with replacement (bootstrap) is drawn. This is done at each TR and the four parameter fitting procedure is then applied to the bootknife samples to obtain a ^T WM1 and a ^T GM1 replications. This operation is done B times (typically we used B = 200), see Fig. 4, and a standard error (SE) can be estimated as follows:

$${\rm{\hat S}}{{\rm{E}}_{{T_1}}} = {\left( {{{\Sigma _{b = 1}^B{{[\hat T_1^b - {{\bar T}_1}]}^2}} \over {B - 1}}} \right)^{1/2}}$$
(2)

, where ¯T 1 is the mean value of the estimates over the B bootknife samples and ^T b1 is the bth estimation of T 1 (or replication) fitted from the bth bootknife sample set.

Fig. 4
figure 4

The bootstrap algorithm adapted from [22] for estimating standard error of metabolite T 1; The T 1 estimate is fitted from the original sample and can be either WM or GM T 1. The bootstrap replications ^T b1 , b=1 … B, are used to calculate its standard error estimate (SE). B is usually between 25 and 200

To take into account the confidence in the T 1 value estimation from a subject when calculating the mean T 1 value over a group of subjects, we proposed using the estimated standard errors in a weighted average calculated as follows:

Weighted sample mean:

$$T_1^* = {{\Sigma _{i = 1}^{{\rm{nbofscans}}}{w_i}T_1^i} \over {\Sigma _{i = 1}^{{\rm{nbofscans}}}{w_i}}}{\rm{with }}{{\rm{w}}_i} = {1 \over {\hat SE_i^2}}$$
(3)

, where nbofscans corresponds to the number of subjects scanned in the study and contributing to the group mean value.

The corresponding individual standard deviation of this mean value is estimated with:

$${\rm{stdev(}}T_1^i) = \sqrt {{{\Sigma _{i = 1}^{{\rm{nbofscans}}}{{(T_1^i - T_1^*)}^2}*{w_i}} \over {\Sigma _{i = 1}^{{\rm{nbofscans}}}{w_i} - {{\Sigma _{i = 1}^{{\rm{nbofscans}}}w_i^2} \over {\Sigma _{i = 1}^{{\rm{nbofscans}}}{w_i}}}}}{\rm{ }}} {\rm{with }}{w_i} = {1 \over {{\rm{\hat SE}}_i^2}}$$
(4)

, where T il is the estimated T 1 for a given metabolite and subject i, ^SE i is the bootstrap estimate of the standard error of the T 1 estimation for the subject i. Figure 5 summarizes the steps used in the proposed method for the case of noisy signals.

Fig. 5
figure 5

Scheme of the proposed procedure for a metabolite T 1 fitting where “SE” stands for bootstrap standard error estimates

Monte carlo simulations

Simulations were conducted to validate the method. The bias and standard deviation of metabolite T 1 estimates were determined using Monte Carlo studies for different N avg, with or without a macromolecular background contamination. We show the ability of the bootstrap technique to estimate the metabolite T 1 standard error. We used previously described MRSI data simulation programs [24]. In this simulation technique, aMRSI k-space is generated by performing the product of the k-space distribution of a given object and its corresponding spectroscopic signal. The effect of the PRESS box selection is also taken into account. For our simulations, one object for WM and one for GM were generated and the results summed. A short echo time signal containing NAA, choline, creatine, Glu, and mI was created for each TR and each WM or GM object. A background signal reproducing the effect of macromolecular contamination in the short echo time signal and mimicking, in the metabolite region, the macromolecule baseline signal found in [3] was also added to the signal. The metabolite relative amplitudes, T 1, SNR and the background components that were used in the Monte Carlo studies are summarized in Table 1. Without noise, for one voxel with a certain WM/GM content and one TR, the simulated metabolite signal amplitude was set exactly to Eq. 1. White Gaussian distributed noise signals were added to each noise-free signal of the simulated MRSI grid. These steps were repeated to obtain a total of N MC = 100 noisy MRSI sets of signals. Metabolite T 1 values were obtained using the proposed method for those 100 realizations. For each metabolite, we obtained a gold standard SE (standard error) on T 1 by calculating the standard deviation of all of the T 1 estimations. This gold standard SE was compared to the mean value of the standard errors calculated by the proposed bootstrap approach for each realization.

Table 1 Parameters used in the Monte Carlo simulation studies

Study 1: Effects of averaging

In the first Monte Carlo study, we evaluated the statistical performance of the method in terms of bias and standard deviation for SNR = 2 along with an increasing number of voxels in an average, N avg = 2, 4 and 6, and the macromolecular background signal (see Table 1) added to the simulated short echo time signals. For each TR, N gen = 120 voxels are created by the proposed averaging method from the original MRSI grid. We also tested our method in the case of no averaging. In this case, designated by N avg = 0, only the 48 voxels (N gen = N total) in the MRSI grid contribute to the fitting method.

Study 2: Macromolecular background effect

We tested the macromolecular background effect on the T 1 estimation in the second Monte Carlo study. We compared the statistical results obtained with N avg = 6 with and without a macromolecule signal in the simulation. The weighted average using the bootstrapped standard error estimation is tested in these two cases.

Results

Monte Carlo studies

Results of the Monte Carlo simulation are shown in Figs. 6, 7 and Table 2. SE, bootstrap estimates of the SE and bias are expressed as a percentage of the true T 1 value.

Fig. 6
figure 6

From Monte Carlo simulations, gold standard standard errors of metabolite T 1s compared against standard errors estimated by the bootstrap technique with varying numbers of voxels used for signal averaging N avg. The sign (*) indicates cases where the standard deviation of the bootstrap estimate is larger than its mean value or where the bootstrap estimate is biased by more than 100% from the actual value. Note the general decrease of the SE with N avg

Fig. 7
figure 7

The bias (b) on metabolite T 1 and the bias for a weighted average (b w) calculated with Eq. 3. b and b w are displayed as percentage of the true metabolite T 1 value

Table 2 From Monte Carlo simulations (N MC = 100), gold standard standard errors (SE), standard error (SEw) corresponding to the weighted average and calculated with Eq. 4, bias (b) and bias for a weighted average (b w) calculated with Eq. 3

Study 1: Effects of averaging

Figure 6 shows the mean and the standard deviation of the bootstrap estimates, as well as the gold standard of the SE for each metabolite, for a range of N avg, SNRNAA equals 2, and the set of 5 TRs. (*) indicates cases where the bootstrap estimate is biased by more than 100 %from the actual value.

For the four metabolites of interest (NAA, Cr-CH3, Cho, mI), higher SE were found in the GM than in WM consistent with the discrepancy between the number of GM voxels versus the number of WM voxels available in the masks we used in the simulation. The SE globally decreases with the number of voxels used in an average N avg. For N avg = 6, the gold standard SE is below 20% of the true T 1 for NAA, Cho and Cre in the WM.

For all of the considered metabolites in the WM and the GM and for N avg ≥ 4, the SE estimated by the bootstrap technique was within 50% of the actual SE value. Less biased bootstrap estimates were generally obtained for N avg = 2 and the bias between the bootstrap estimates and the actual SE value increased with N avg. For mI, the bootstrap approach successfully estimated the SE with a large standard deviation in the WM and only for N avg = 6 in the GM. For N avg = 0, the bootstrap estimates have reasonable values in WM for NAA, Cho and Cre and failed to estimate the SE in the other cases.

Figure 7 shows the bias (denoted by b) on T 1 values for the four metabolites of interest and for a weighted average calculated by Eq. 3 and denoted by b w. b and b w are displayed as percentages of the true metabolite T 1 value. For (WM/GM)-T 1,NAA (WM/GM)-T 1,Cr-CH3 and WM-T 1,Cho the bias is below ±10% of the true T 1 value. We note that, in the case of our simulation, increasing the N avg did not reduce the bias for WM-T 1,Cr-CH3. As seen in the next study, we think that this bias is more due to the interaction with the macromolecules than to a lack of SNR. Increasing N avg seems to reduce the bias for (WM/GM) T 1,NAA, GM-T 1,Cho, GM-T 1,Cre, and (WM/GM) T 1,mI. For WM and GM-T 1,mI, the weighted average makes the bias below 10% of the actual value for N avg ≥ 2.

The best bias and standard deviation trade-off is obtained for N avg = 6 at the expense of slightly biased bootstrap estimates.

Study 2: Macromolecular background effect

Table 2 shows the statistical results (regular SE, SEw calculated with Eq. 4, b and b w) of the T 1 estimation procedure for N avg = 6, SNR = 2, with and without a background contamination in the signals to process. All the metabolite T 1 estimations show the same trend with a larger standard deviation and a bigger bias in the case of macromolecular contamination as compared to the absence of a background signal. In the absence of a macromolecular background, the weighted average that takes into account the bootstrapped standard error successfully reduces the bias (b w as compared to b) on the metabolite T 1.

In the presence of the simulated background signal, the SE of the T 1 estimate increased by 10%(for WM-T 1,Cr-CH3) and more than doubled for the T 1,mI. The bias was also affected (but usually reduced) by using a weighted average. The GMT 1,Cr-CH3 and WM-T 1,Cho estimation were more downward biased in the presence of a background signal when using a weighted average than when using the standard mean calculation.

Metabolite T 1 estimation on in vivo data

For the eight subjects, the anterior-posterior center of the PRESS Box (12 × 12, 1 cc) was centered in the anterior-posterior middle of the corpus callosum body visualized in a sagittal plane. Voxels anterior to the center of the PRESS Box were analyzed as part of the anterior brain region (Ant.) and the rest of the voxels were evaluated as posterior voxels (Post.).

The bootstrap estimate of standard error on in vivo data

Figure 8 shows, for each subject of the study, a histogram of B = 200 bootstrap replications of ^T WM1,NAA calculated in the posterior part of the PRESS box. For each subject i(i = 1,…, 8), these replications are used to estimate standard error SE i . In this example, subject 3 presents a small standard error and thus will have a larger weighting in the proposed weighted average in Eq. 3 while the results from subject 2 or 5 will have smaller weighting due to their bigger bootstrap estimated standard error. This bootstrap standard error estimate gives good insight about the reliability of the fitted metabolite T 1. In the case of subject 8, the bootstrap histogram reveals a peak quite distinct from the average when fitting ^T WM1,NAA from the selected voxels. The dotted line indicates the estimated T 1 when all the selected voxels are considered while the bootstrap histogram is calculated with the bootknife approach that randomly removes one voxel at each calculation. Therefore, this case shows that some voxels were essential in determining the T 1.

Fig. 8
figure 8

Histograms of B = 200 bootstrap replications of WM-T 1,NAA, calculated from in vivo data (posterior region) of eight healthy subjects. A broken line is drawn at the parameter T 1 estimate

Anterior versus posterior metabolite T 1

The estimated in vivo T 1 relaxation times (mean±SD) for “pure” WM and “pure” GM at 1.5 T for N avg = 6 and computed from the anterior region (Ant.) or from the posterior region (Post.) are given in Table 3. These results are calculated with Eqs. 3 and 4 using the bootstrap standard error estimates and are given as weighted means ± the corresponding standard deviations. For statistical analysis, a two-tailed paired t-test was used. See Appendix.

Table 3 From eight healthy volunteers, estimated T 1-relaxation times (in seconds) of NAA, Cho, Cr-CH3, mI at 1.5 T, in pure WM and pure GM, using MRSI data, (N avg = 6)

The standard deviation in the GM was larger than in the WM, for most of the time, as in the Monte Carlo simulations. From these results, no significant differences were found in metabolite T 1s between GM and WM in the posterior region. In anterior region WM-T 1,NAA, \(WM-T_{1,Cr-CH_{3}}\), and WM-T 1,mI tend to be higher than GM-T 1s but this trend did not reach statistical significance (0.1 < P < 0.2). The T 1 of NAA in the anterior part of the WM was significantly longer than in the posterior part of the WM (P < 0.05).

Discussion

This work presents a novel method for estimating metabolite T 1s using MRSI data, enabling estimates within tissue types (GM and WM) and across different regions (anterior and posterior were demonstrated here). The smaller voxel size of the MRSI data versus previous single voxel studies partially addresses the issue of large partial volume artifacts between gray matter and white matter. Additionally, incorporating the information regarding tissue type composition obtained from higher resolution MR images and the multiple voxel data obtained with MRSI allows better correction of partial volume artifacts than possible with the previous single voxel data. This method has been validated and tested on simulations and applied in vivo. We also proposed assessing a confidence interval in the fitted T 1 results by introducing a bootstrap approach. We showed that a weighting average that takes into account the confidence assessment can reduce the bias on the estimated T 1 value for a metabolite with a small SNR. This method relies upon a group mean metabolite T 1 approach. The proposed algorithm yielded results that are in agreement with the literature and support the hypothesis of regional differences in T 1 in the brain.

Methodology

The important parameters in the proposed method are the SNR of the considered metabolite signal, the number of voxels used in an average, N avg, and the number of available, mostly GM or mostly WM voxels in the MRSI grid.

  1. 1.

    The SNR, as expected, appeared to clearly play a role both on the bias and the standard deviation of the T 1 estimation. In Figs 6 and 7 the results for NAA and Cr-CH3 which have the greatest SNR in our simulation, present good biases and standard deviations on T 1 estimation compared to the ones for Cho (especially in GM) and mI. In order to realize a robust four-parameter fit with low SNR, non-reliable voxel quantification results should be rejected from the analysis. The use of criteria, such as Cramér-Rao lower bounds [18] or linewidth thresholds, is necessary to determine the quality of the metabolite amplitude quantification and to perform voxel selection.

  2. 2.

    Increasing the number of voxels used in an average N avg, improves SNR and so typically reduces both bias and standard deviation. Note that, by averaging the voxels, the independence between the averaged voxels is reduced and the bootstrap technique can tend to underestimate the real standard error. Also note that this averaging is performed while taking into account the WM-GM distribution and the introduced dependence tracks the initial tissue content. It was also shown by the Monte Carlo simulation results that, for an original SNR of two, the use of six voxels in an average corresponds to a good trade-off between the SNR gain and the lost of voxel independence and leads to reliable metabolite T 1 estimations. Furthermore, this voxel averaging introduces some partial voluming with CSF to the generated voxels, especially for mostly GM voxels originating from the thin cortical ribbon at the midline. Nevertheless, as the percentage of GM and WM (and thus, CSF) in the voxel are explicitly taken into account during the fitting procedure, the regression presented in Eq. 1 enables a correction of this partial voluming effect.

  3. 3.

    The number of voxels available for a specific tissue type influences the standard deviation of the T 1 estimation. Consequently, the dispersion of the results is larger for the GM than for the WM values from both our simulation and the in vivo data.

The first Monte Carlo simulations showed that for an original SNRNAA of 2 in the MRSI data, an almost unbiased estimation of the T 1,NAA, T 1,Cho, T 1,Cr-CH3 is possible. A weighted average, taking into account the bootstrap estimates of the uncertainty on T 1 could also yield to an unbiased estimation of T 1,mI. This approach tends to reduce the dispersion due to bad data points and makes a group mean value more accurate.

We conclude from the second Monte Carlo study that the presence of a macromolecular background signal has an important effect on the dispersion and the bias of the results and should be considered as another effect, besides the set and number of TRs and the type and parameters of the sequence used [25], leading to the discrepancy between the different published T 1 values. The bootstrapped standard error estimation is also hampered by the macromolecular contribution. In the proposed simulation, all the voxels were contaminated in the same way and the background signal was arbitrary and particularly elevated under the mI making its T 1 estimation harder.

Finally, the proposed bootstrap procedure is a novel method to estimate metabolite T 1 standard error and is enabled in this study by the use of MRSI data. The bootstrap technique is a non-parametric method that does not require a complex implementation. This approach may be beneficial in future uncertainty and/or bias estimation studies of spectroscopic quantification.

In vivo results

The range of our results for the four metabolites is in good agreement with existing literature [14], when using the weighted average and N avg = 6. The method gives reliable results, for an original SNR around 2, when processing equivalent 6 cc voxels, which corresponds to an SNR around 4.9 for NAA. Voxel averaging combined with the quantification procedure, which constrained the additional damping factor allowed for each metabolite, was required to achieve the accuracy of the results.

We found significantly greater WM-T 1,NAA values in the anterior (frontal) part (1.38 ± 0.15, mean ± SD) than in the posterior part of the slice (1.28 ± 0.10), but we were not able to see this result for the GM or for the other metabolites. Brief et al. [1] also reported similar results between WM frontal (1.59 ± 0.10), and WM parietal (1.35) regions. In the literature, WM-T 1,NAA values can differ to some extent. Our anterior WM T 1,NAA is lower than the one reported by Brief et al. or Kreis et al. [3] (1.88 ± 0.09), but still higher than others, as for example the value reported by Ethofer [2] (1.19±0.09 in the fronto-parietal region). This discrepancy may be partly due to the way the macromolecular background signal was fitted.

The effect of regional variation in metabolite T 1 values may be necessary to take into account in the estimated metabolite concentration, depending on the ratio of TR/T 1 used, on the ratio of the regional T 1s, and on the accuracy of the metabolite amplitude estimation. In our study, a regional T 1 variation for NAA of 7.8% would result in a difference of only 5% in the NAA signal amplitude for a short TR of 1 s. Considering the biological variability and the accuracy achievable for the NAA signal amplitude estimation, this difference might be negligible for a concentration estimation point of view. In the case of healthy versus diseased brain, the metabolite T 1 difference may be larger than 7%. The measured NAA signal amplitudes could have important differences (greater than 5%), due solely to the T 1 variation and not to a tissue concentration difference. Conversely, the T 1 differences could mask the concentration changes due to the disease. Of course, when the repetition time exceeds three expected T 1(TR > 3T 1), the difference in regional T 1 can range from 0 to as much as 50%, and the difference (due to the T 1 weight) in metabolite signal amplitude will remain below 5% of the actual concentration value. Then the regional T 1 variation will have effectively no effect on the estimated concentrations. In practice, the use of a long TR increases scan time, especially for MRSI acquisition and is therefore avoided.

The T 1 value found for NAA in the WM of the posterior part of the slice (1.28 ± 0.10) is close to the values reported by Rutgers et al. [5,26] (1.30 ± 0.14) in the centrum semiovale. The T 1 relaxation times found for the other metabolites Cho, Cr and mI, are also essentially the same as other reported values [2,5,26].

Especially in gray matter, where the glutamate signal is higher, the macromolecular signal, the NAA and the glutamate signals and some contribution from metabolite present at low concentration such as GABA are unknowingly entangled at 2 ppm. Moreover, the amount of macromolecule signal compared to the metabolite signal differs at each TR, as macromolecules have a shorter T 1 than metabolites. As a consequence, the variability of the quantification results increases. We think that the higher standard error found for the GM-T 1,NAA is partly due to this variability and partly due to the few number of available gray matter voxels.

We also observed without reaching statistical significance that, as opposed to water T 1, the T 1 for NAA, Cr-CH3 and mI could be greater in the white matter than in the gray matter in the anterior part while no difference was seen in the posterior part. Although a difference of WM/GM voxel distribution in the anterior and posterior regions (see Fig. 1) might have influenced this result, this observation supports the assumption of different underlying mechanisms for water and metabolite relaxation times. While tissue composition and difference in anisotropy may be involved in water T 1 relaxation process, the intra-cellular metabolite T 1, may be more dependent, as suggested by Ethofer et al., on micro-structural characteristics and viscosity properties.

Conclusion

Brain metabolite T 1 measurements were calculated using a novel MRSI voxel averaging and bootstrapping approach. The proposed method takes advantage of the multi-voxel acquisition provided by MRSI and enables the investigation of regional variations in metabolite T 1 values. It also introduces a bootstrap technique for estimating a standard error on metabolite T 1s. Significant differences were found between anterior WM-T 1,NAA and posterior WM-T,1,NAA. This result emphasizes the need to take into account tissue and regional T 1 differences in MRSI metabolite quantification. Finally, the method only requires a multi-WM/GM voxel acquisition and is not restricted to short echo time 2D MRSI acquisition. The presence of a macromolecular background made the metabolite T 1 estimation less accurate and substantially increased the dispersion. The principle of the method can be applied and extended to other field strengths (to increase the SNR) or to other types of data acquisition that present less macromolecular contamination such as TE-averaging [27], longer TE acquisition or to using localization in a third spatial dimension (3D-CSI).