PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Magn Reson Med. Author manuscript; available in PMC 2010 October 1.
Published in final edited form as:
PMCID: PMC2930376
NIHMSID: NIHMS219817

High Dynamic Range MRS Time-Domain Signal Analysis

Abstract

In the absence of water signal suppression, the 1H MRS in vivo water resonance signal-to-noise-ratio (SNR) is orders of magnitude larger than the SNR of all the other resonances. In this case, because the high SNR water resonance dominates the data, it is difficult to obtain reliable parameter estimates for the low SNR resonances. Herein, a new model is described that offers a solution to this problem. In this model, the time-domain signal for the low SNR resonances is represented as the conventional sum of exponentially decaying complex sinusoids. However, the time-domain signal for the high SNR water resonance is assumed to be a complex sinusoid whose amplitude is slowly varying from pure exponential decay and whose phase is slowly varying from a constant frequency. Thus, the water resonance has only an instantaneous amplitude and frequency. The water signal is neither filtered nor subtracted from the data. Instead, Bayesian probability theory is used to simultaneously estimate the frequencies, decay-rate constants, and amplitudes for all the low SNR resonances, along with the water resonance’s time-dependent amplitude and phase. While computationally intensive, this approach models all of the resonances, including the water and the metabolites of interest, to within the noise level.

Keywords: magnetic resonance spectroscopy, NMR, MRS, in-vivo, signal modeling, water suppression, water filter, parameter estimates, Bayesian, probability theory

INTRODUCTION

The signal-to-noise ratio (SNR) of the water signal in hydrogen (1H) magnetic resonance spectroscopy data (MRS) from intact biological systems can be orders of magnitude greater than the SNR of the metabolite resonances of interest. The water resonance lineshape in these systems can be broad and asymmetric due to tissue magnetic susceptibility gradients and insufficient magnet shim capabilities. Instrumental instabilities present during data acquisition can cause the water signal’s phase and amplitude to vary with time. The result is a water signal with a complicated, non-exponential, non-single frequency time-domain evolution. The water resonance’s high SNR and non-ideal behavior, combined with the modest resonance frequency dispersion found at static field strengths common to 1H MRS in vivo, make it challenging to obtain accurate parameter estimates for the frequencies, amplitudes, and decay-rate constants of the low SNR metabolite resonances.

Of particular note, in the frequency domain, the wings of the high SNR water lineshape extend across and beyond (thus, alias within) the spectral bandwidth, and underlie the low SNR resonances. Increasing the bandwidth minimizes aliasing effects but increases the noise. Analysis in the time domain avoids this problem because these effects are implicit in the model, but it requires a model for the water signal capable of modeling the water down to the noise level. However, it is impossible to model a high SNR water signal down to the noise level by exploiting analytical expressions for the water lineshape such as a pure exponential, Gaussian, or a combination thereof because these models ignore phase variations, and are inherently symmetrical in the frequency domain.

The development of pre-acquisition and post-acquisition methods for dealing with the extreme dynamic range of 1H MRS in vivo is an active area of research. Broadly considered, researchers take two approaches: (i) pre-acquisition selective suppression of the water resonance and (ii) post-acquisition modeling or filtering of the water resonance. Pre-acquisition water signal suppression offers the benefit that receiver gain can be increased to take advantage of the spectrometer’s full dynamic range capabilities while reducing frequency-domain lineshape overlap with the low SNR metabolite resonances. Post-acquisition modeling offers the benefit that the water resonance can serve as an internal amplitude, frequency, and phase reference to aid in the analysis of the metabolite resonances. Collecting data without water suppression also avoids possible confounding effects of magnetization transfer between water and metabolite resonances via chemical exchange and cross relaxation (1, 2). Importantly, both approaches seek to avoid perturbation of the signal components defining the metabolite resonances.

1H MRS techniques to either suppress the water magnetization (3-5) or enable data collection without water suppression (6-9) are well known and have been elegantly extended and exploited by many researchers, vide infra. This manuscript describes a new strategy for post-acquisition modeling of a single, high SNR water resonance simultaneous with modeling the metabolite resonances. While motivated by in-vivo experiments ongoing in our laboratory, the method is applicable to any MR free-induction decay (FID) where one resonance has a SNR that is orders of magnitude greater than the other resonances.

A variety of post-acquisition schemes are available for quantification of low SNR metabolite 1H resonances in vivo. While the water resonance can be modeled in the frequency domain (10, 11), there are advantages to time-domain modeling or filtering (12). Most time-domain methods either explicitly model the water resonance and subtract it from the experimental FID (12, 13) or mathematically filter the FID (14-17) to selectively remove the water resonance. These methods seek to minimize the effect of the water resonance on estimates of metabolite signal parameters. Wavelet transformation techniques (18, 19) and the Gabor transform (20) are time-scale and time-frequency analysis tools, respectively, that also seek to remove the time-domain water signal from the data.

Input assumptions for these various post-acquisition methods vary. In some instances, one must know the metabolite resonance frequencies or select an appropriate filter bandwidth. Almost all approaches assume the water-resonance decay profile is well-described as a sum of mono-exponentially decaying sinusoids (12). Cabanes et al., have explored this approach using a finite basis set of exponentials (21). A more complex analytical model, Voigt decay, treats the water signal decay as a convolution of a Lorentzian (exponential in t) and a Gaussian (exponential in t 2). The Voigt-decay model (22) assumes the signal’s evolution produces a frequency-domain lineshape that is symmetrical about the resonance frequency, a situation often violated in practice.

Methods that filter the time-domain data reduce its information content. Filtered data are a weighted average of the filter-response function convolved with the spectrum of the original data. As desired, such averaging predominantly affects the water signal. However, all the other signals in the data are also affected by the filter. Filtering always affects the metabolite resonances in an undesirable way and all successful filter-based methods incorporate procedures that attempt to minimize these distortions.

This paper describes a time-domain method that simultaneously computes parameter estimates for a high SNR resonance (water) and low SNR resonances of interest (metabolites). In brief, a new model is proposed for the time-domain water resonance. The water signal is modeled as a complex sinusoid whose amplitude and phase are slowly varying functions of time (i.e., as a time-domain signal that deviates slowly from the ideal of an exponentially decaying, constant frequency resonance). The time-dependent component of the water signal’s phase and amplitude is assumed to be a non-harmonic, non-periodic modulation. The water-resonance model is combined with a conventional model for the low SNR metabolite resonances of interest. The water signal is neither filtered nor subtracted from the data. Instead parameter estimates for all the modeled resonances are computed concurrently within a Bayesian framework.

THEORY

The problem to be addressed is a parameter estimation problem in which there is a single high SNR water resonance and a known number of low-SNR metabolite resonances. We note that this approach could be made even more general with the additional estimation of the number of metabolite resonances in the data. Additionally, extension to more than one high SNR resonance would also be straightforward.

As in all Bayesian analysis problems, one begins by relating the data to the hypotheses of interest. In this problem, the hypotheses of interest involve the true, but unknown, values of the frequencies, decay-rate constants, and amplitudes of the metabolite resonances, as well as the parameters describing the water resonance. The relationship between the data and the parameters of interest is established through a model consisting of three parts: a model for the metabolite resonances, a model for the water resonance, and models for constant offsets that may be present in the data.

First, the metabolite resonances are modeled as a sum of complex, exponentially decaying sinusoids:

SM(ti)k=1mAMkΓk(ti)
[1]

where m is the number of known metabolites (23). The time-domain amplitude of the kth metabolite is denoted as AMk. Note, this amplitude is proportional to the integral over the metabolite lineshape in the frequency domain. The Γk (ti) are complex exponentially decaying sinusoids and are given by

Γk(ti)exp{2πfMki[ti+t0]αMkti+iθ}.
[2]

The frequency of the kth metabolite resonance is designated as fMk, its decay-rate constant is αMk. The phase offset of the sinusoid is θ (the zero-order phase) and t0 is a time delay that manifests itself in the frequency domain as a frequency-dependent (first-order) phase. The minus sign is a convention used to define what is meant by a positive frequency in the data acquisition scheme.

In our experience, the model described by Eqs. [1] and [2] is adequate for the relatively low SNR metabolite resonances in vivo. However, this model does not accurately describe the high SNR water-resonance signal. As the SNR increases, time dependencies in the signal phase and amplitude become relevant and the signal model described by Eqs. [1] and [2] can no longer model the water signal into the noise. To accommodate this signal behavior, an alternative model is proposed wherein the water resonance’s phase and amplitude vary in a nonlinear, time-dependent fashion. The challenge is to model only those amplitude and phase modulations supported by the data, without over-fitting the water. The metabolite and the water resonances experience the same non-periodic, non-harmonic modulations, but the metabolite signals’ relatively low SNR means that Eqs. [1] and [2] can model the metabolite signals to within the noise level.

The second model component, the model for the high SNR water signal, SW (ti), is

SW(ti)A(ti)ΓW(ti)exp{iϕ(ti)},
[3]

where A(ti) and ϕ(ti) describe the water signal’s time-dependent amplitude and phase. The water-signal’s complex exponentially decaying sinusoid component, ΓW (ti), is given by

ΓW(ti)exp{2πifWtiαWti},
[4]

where fW and αW would be the water resonance’s frequency and decay rate constant in the absence of a time-domain amplitude and phase modulation. A time delay, t0, is not needed in Eq. [4] because it is accounted for in ϕ(ti). The instantaneous water frequency, fW (t), is given by

f^W=ddt[fWtϕ(t)2π].
[5]

Using exp {iϕ(ti)}[equivalent] cos[ϕ(ti)]+ isin[ϕ(ti)], Eq. [3] is rewritten in a more convenient form:

SW(ti)=A(ti)(cos[ϕ(ti)]+isin[ϕ(ti)])ΓW(ti).
[6]

Both A(ti)cos [ϕ(ti)] and A(ti)sin [ϕ(ti)] are arbitrary functions of time, and replacing them with the functions Ac (ti) and As (ti) gives

SW(ti)=[Ac(ti)+iAs(ti)]ΓW(ti)
[7]

as the water-signal model. The two functions Ac (ti) and As (ti) express empirically how the actual water signal deviates from a pure complex exponentially decaying sinusoid. To physically model the amplitude and phase modulation of the water, these functions must have some simple characteristics. For example, Ac (ti) and As (ti) should vary continuously and they should not be highly oscillatory. Simple functions having these properties are low-order polynomials, and Ac (ti) and As (ti) can be expressed as simple polynomial expansions of the form:

Ac(ti)=j=1ncAcjLj(ti),
[8]

and

As(ti)=j=1nsAsjLj(ti),
[9]

where nc and ns, are the number of the polynomials used in the expansion, Lj are polynomials, and Acj and Asj are the polynomial expansion coefficients.

Time-domain NMR data can contain constant offsets and in many cases, the first point of the complex data can be contaminated by instrumental artifacts. Additionally, spectrometer manufacturers often modify the first data point by setting its value to half its original complex data value. Regardless of the cause, the first data point in the NMR FID often needs to be modeled separately. Consequently, we introduce a third model component, and we refer to these models as constant models, although strictly this is not the case. This third model component is given by an offset in the real channel, an offset in the imaginary channel, a real first-point offset and an imaginary first-point offset. These four models are designated as SC (ti). We have indicated these models are time dependent because the first data point models are time dependent.

The full time-domain model is a sum of the three different model components, SM, SW and SC, and is given by

di=SM+SW+SC+ni
[10]

where a complex datum is represented as di and ni represents a complex noise value.

In Bayesian probability theory, all of the information about a parameter given the data and the prior information is summarized in a posterior probability density function. For instance, a plot of the posterior probability density function for a metabolite resonance’s frequency represents the degree of belief that a particular value of the parameter is the true, but unknown, value of the parameter. The frequency estimate can be found from the distribution’s mean or its peak value, while the uncertainty in the frequency estimate is related to the distribution’s width. Derivations of the posterior probabilities for the parameters of interest are detailed in the Appendix.

The Bayesian calculations are implemented using a Markov chain Monte Carlo (McMC) algorithm with simulated annealing to draw samples from the joint posterior probability for all the non-linear parameters, i.e., everything except the amplitudes and the constants. Monte Carlo integration is then used to generate samples from the posterior probability for the individual parameters. The reported parameter estimates and uncertainties are the means and standard deviations of these samples.

RESULTS AND DISCUSSION

This section discusses the application of the techniques described in this paper to data generated in silico (i.e., computer simulated) and to data obtained experimentally in vivo. The simulated data consists of a high SNR “water” signal, plus three low SNR signals from “metabolites”. The metabolite signals are simulated as complex exponentially decaying sinusoids, while the water signal is simulated as having both a time-dependent change in frequency (a chirp) and a Gaussian decay profile. Thus, the simulated data allow one to determine how well the metabolite parameters can be determined even though the model used to generate the simulated water signal does not have the same analytical form as that used in the analysis calculations. The in vivo data were obtained on a Varian NMR Systems (Palo Alto, CA) INOVA small-animal MR scanner operating at 4.7 T using localization by adiabatic selective refocusing (LASER) pulse sequences (21).

The estimates for the water and metabolite resonance parameters are made directly from the time-domain data. However, it is convenient to display these results in the frequency domain. To make these frequency-domain displays, the parameters from the Markov chain Monte Carlo simulation having maximum posterior probability were used to generate a time-domain model of the data. This time-domain model, the original data, and the residuals (i.e., the difference between the data and model) were then Fourier transformed, phased and displayed. In the text that follows, it is these Fourier transforms that we will be discussing.

Parameter Estimates from Simulated Data

This subsection describes parameter estimation results from two computer-simulated FIDs using the methods described in the Theory and Appendix Sections. Table 1 lists the frequency, decayrate constant and amplitude input values for the water and three metabolite signals used to generate the FIDs. In the first case, noise was added, σ=80, to give a simulated water SNR of 11,600. In the second case, the noise level was increased, σ=1,000, to give a water SNR of 930. The SNRs for the lowest amplitude simulated metabolite signal, metabolite 2, are 9.6 and 0.8 in the σ=80 and σ=1,000 cases respectively. The metabolite signals’ decay profiles are pure mono-exponential. The water signal’s decay profile is Gaussian. The total phase of the water signal has two components. One component varies linearly with time and the other varies quadratically. The quadratically varying phase changes by π over the duration of the signal. The amplitude of the first data value in each of the two quadrature channels was changed (halved) and constant offsets are present in both channels.

Table 1
Bayesian Parameter Estimates From A Computer-Generated FID

Table 1 summarizes the parameter estimates and uncertainties for the frequencies, decay rate constants, and amplitudes of the water and metabolite signals in the σ=80 and σ=1,000 FIDs. A water signal decay-rate constant estimate is not reported, as the water decay profile is Gaussian, not exponential. The water SNR is so high that the parameter estimates for its signal in both FIDs are nearly identical. The actual error in the metabolite parameter estimates for the σ=80 FID are well within one standard deviation.

Metabolite amplitude estimate quality is crucial for many MRS applications in vivo. Although not meant to be an exhaustive study, these results suggest that useful amplitude estimates can be obtained with metabolite resonance SNR of ~10:1. The estimated uncertainties for the σ=1,000 case are much higher. But even with a SNR of 1:1, the parameter estimates for Metabolites 1 and 3 are within one standard deviation of the input values. For metabolite 2, whose SNR is less than 1:1, the decay-rate constant and amplitude estimates are within two standard deviations, while the frequency estimate is very close to the input value.

A visual summary of results from the computer-simulated FID analyses is shown in Fig. 1. Figures 1a and 1b depict the σ=1,000 FID’s DFT. The base of the water signal completely overlaps two of the metabolite signals. Figures 1c and 1e show DFTs of FIDs computed from the Bayesian parameter estimates: Fig 1c is the total signal and Fig 1e is just the metabolite signals. Figure 1d is the absorption-mode residual signal obtained by subtracting Fig. 1c from Fig. 1b. The residual spectrum shows how the model described in Eq. [10] can quantitatively characterize a large-amplitude signal with non-exponential decay, estimating parameters for the low SNR metabolite signals in Eq. [1] while accounting for the non-exponential decay of the large amplitude water signal.

Figure 1
The results from the Bayesian analysis of a computer-simulated FID. The simulated FID mimics a localized MRS voxel of brain white matter from the Corpus callosum in vivo without water suppression. There are three metabolite resonances. The water and metabolites’ ...

Parameter Estimates from Mouse-Brain White Matter in Vivo

Monitoring N-acetylaspartate, NAA, choline containing compounds, CHO, and total creatine, CRN, levels in vivo has relevance to the study of neurological disease states involving neuronal loss, axonal dysfunction and malignancy. The lower and middle traces in Fig. 2a show a mouse-brain white-matter spectrum from the Corpus callosum acquired in vivo without water suppression using a TE of 100 ms and a TR of 1s. The voxel size was 1 × 2 × 4 mm3. The high SNR water resonance dominates the spectrum. The vertical scale for the middle trace in Fig. 2a is increased 200 × to display the metabolite peaks, which are clearly overlapped with the water resonance. For data analysis, the signal model employed is water plus three metabolite resonances – NAA, CHO, and CRN. The upper trace in Fig. 2a depicts the sum of the individual model resonances for water, NAA, CHO, and CRN metabolite methyl signals. Figure 2b displays the 1-6 ppm region of Fig. 2a. The bottom and middle traces in Fig. 2b are the models for the three metabolite resonances and the water resonance, respectively. Note the asymmetric nature of the water-model frequency-domain lineshape. Figure 2c is the residual power spectrum. The residual’s absorption-mode display typically contains out-of-phase signals, so the residual power spectrum is a useful means to determine reasonable metabolite resonance candidates in subsequent analyses. The lack of water signal in the residual demonstrates how well Eq. [10] models the metabolite signals, constants, and an asymmetrical water resonance. Table 2 summarizes the Bayesian parameter estimates of interest. The computed (25, 26) NAA, CHO and CRN 1H concentration estimates based on a water proton equivalent 1H concentration of 80 M (7) are also listed in Table 2. While these data were not collected under rigorous quantitative conditions, the metabolite amplitude estimates are consistent with literature reports (27-29).

Figure 2
The DFTs of the experimental and model FIDs from a 1 × 2 × 4 mm3 voxel of mouse-brain white matter from the Corpus callosum in vivo. The 512 scan FID was acquired without water suppression. The traces in Figs. a and b are displayed with ...
Table 2
Bayesian Concentration and Parameter Estimates for Mouse-Brain White-Matter In Vivo.

Parameter Estimates from Mouse-Kidney Tissue in Vivo

The concentration of trimethylamines, TMA, osmolytes within kidney tissue extracts, is a marker for polycystic kidney disease in Han:SPRD rats (30). The TMA methyl resonance found at ~ 3.3 ppm in vivo primarily represents the osmolytes betaine and glycerophosphorylcholine. Osmolyte levels fluctuate with the degree of kidney-tissue hydration (31, 32). Collecting kidney MRS spectra in vivo without water suppression is useful, as the amplitude of the water signal can serve to normalize the TMA concentration. Unfortunately the TMA resonance is found within 1.5 ppm of the water resonance frequency, and the water signal can interfere with TMA amplitude estimates.

Figures 3a and 3b show a respiratory-gated spectrum of mouse-kidney obtained in vivo without water suppression using a TE of 30 ms and a TR of 1.5 s. The voxel size was 3.5 × 3.5 × 3.0 mm3. In this spectrum the wide base of the water resonance lineshape completely overlaps the TMA resonance, even though their resonance frequencies are separated by 300 Hz at 4.7T. Figure 3c depicts the model for the water, TMA and lipid (CH2)n resonances. The water and TMA amplitude estimates are 1780 ± 10 and 8.5 ± 2.1, respectively, giving a water:TMA amplitude ratio of ~210.

Figure 3
The DFTs of the experimental data, model and residual FIDs from a 3.5 × 3.5 × 3 mm3 mouse-kidney voxel. The 4 scan FID was collected in vivo with respiratory gating. Figures 3a-3c are displayed with a small vertical offset for clarity. ...

The experimental, Fig. 3b, and residual, Fig. 3d, spectra have an in-phase feature at ~ 6.8 ppm and the residual spectrum contains an anti-phase feature at ~ 2.7 ppm. These features are water-resonance sidebands caused by harmonic vibration in the gradient coils (33). The gradient-coil motion induces a frequency modulation of B0 (34). Neither Eq. [6] nor Eqs. [1] and [2] model vibration-induced harmonic signal modulation. While not pursued here, a model for gradient modulation (33) could be added to Eq. [10] to identify and quantify the gradient sideband signals.

Parameter Estimates from Mouse-Muscle Tissue in Vivo

While MR spectroscopy is a convenient method for monitoring muscle-tissue triglyceride content in vivo, obtaining muscle-tissue spectra near bone is problematic. The difference between the magnetic susceptibilities of bone and muscle tissue produces field gradients that degrade the local field homogeneity. Figure 4a (solid line) shows the spectral region between 3-6 ppm, focusing on the unsuppressed water resonance. Data were obtained in vivo from a 2.5 × 2.5 × 3 mm3 voxel of mouse muscle tissue located close to the spine. The TE was 30 ms and TR was 1.5 s. The susceptibility induced field inhomogeneity is so severe that the frequency-domain water resonance appears characteristic of three overlapping resonances, two high-SNR, high-frequency peaks and at least one weaker, low frequency peak. The dashed line, slightly offset for clarity, in Fig. 4a is the result of modeling the water resonance using Eq. [7]. The residual is shown in the top trace of Fig. 4a.

Figure 4
The water spectral region of the experimental data, signal model, and residual from a 2.5 × 2.5 × 3 mm3 mouse-muscle tissue voxel in-vivo. The voxel is adjacent to the spinal chord. The 16 scan FID was collected without water suppression ...

In the frequency domain, the presence of broad resonances from lipids and macromolecules can confound parameter estimates for nearby secondary metabolite resonances. The method described here can directly incorporate rapidly decaying signals from lipid or macromolecular resonances. When the residual, the difference between the model and the data, displays features consistent with rapidly decaying resonances, these resonances are simply added to the model using appropriate frequency and decay-rate constant priors. If lipid or macromolecular signals are not visible in the residual, they have no effect on the parameter estimates.

Figure 4b depicts the full spectrum of the muscle tissue with an expanded vertical scale. The dashed line is the model spectrum for the water and the two highest amplitude lipid resonances. The water-resonance lineshape in the model spectrum has a broad high-frequency component. This type of broad asymmetrical lineshape distortion can be confused with an improperly phased absorption-mode spectrum. The approach described here completely avoids reliance on frequency-domain phase adjustments or ad-hoc baseline manipulations to obtain resonance amplitude and frequency estimates for the metabolite resonances. It is possible to obtain lipid and water resonance parameter estimates (and their uncertainties) for muscle-tissue voxels adjacent to bone even when a high-degree of water-resonance lineshape distortion is present.

It is instructive to examine the mouse muscle FID water resonance and its model signal. The water spectral region in the DFT of the mouse-muscle FID contains two high SNR signal maxima. Nonetheless, a single complex, decaying sinusoid with a nonlinear time-dependent phase and amplitude models the water resonance to within the noise. Figure 5a depicts the real component of the experimentally acquired quadrature FID during the initial 28 ms of data acquisition. The line that intersects zero on the ordinate in Fig. 5a is the residual from the experimental data and the water and lipid models. Figures 5b - 5d graphically illustrate how the water model, Eq. [7], differs from a pure exponentially decaying sinusoid. The fractional change in the real component of the signal model’s complex time varying amplitude Ac(ti) is shown in Fig. 5b. Recall from Eqs. [4] and [7] that the effects of the complex time varying amplitude are independent of the changes in signal intensity described by the transverse decay constant, αW. Thus, a time-independent amplitude would produce a horizontal line with the ordinate equal to zero (i.e., zero fractional change). Instead, the amplitude changes significantly during the first 10 ms. Figure 5c shows that the total accumulated phase of the water model also has a significant, nonlinear, time-dependent component. Figure 5d tracks the instantaneous frequency (more precisely the change in frequency), the derivative of the total phase of the water signal (see Eq. [5]). Figure 5d is, thus, a graphical representation of how the water signal’s frequency deviates from the ideal case of a single constant frequency.

Figure 5
The initial 28 ms of the mouse-muscle tissue FID and its model signal depicted in graphical form. (a) The experimental real-channel data and the time-domain residual signal from the water and lipid resonance models shown in Fig. 4. The residual signal ...

Assuming that water-suppression procedures are not employed, the water resonance can serve as an internal amplitude reference. Given the water signal model, the amplitude estimated at t = 0 is an appropriate reference parameter. With respect to a reference frequency, one could chose either the frequency estimate at t = 0 or the estimated amplitude-weighted frequency integrated over time. One might imagine that the water signal model profile could serve as a more precise metabolite resonance model than the pure exponential decay model described in Eq. [1]. While true in principle, in practice the metabolite SNR is generally insufficient to support such a complex signal model.

CONCLUSIONS

This work illustrates a new approach for obtaining parameter estimates from 1H spectra acquired in vivo. The method differs from previous approaches in several significant ways. Values for the parameters of interest, including those describing the water signal’s time-dependent phase and amplitude, are estimated using Bayesian probability theory. The high-SNR water resonance is modeled as a single complex sinusoid whose amplitude and phase are slowly varying functions of time. The original data are not modified by post-acquisition preprocessing because the model for the data contains all the parameters necessary to model the resonances to within the noise. The water signal is not filtered or subtracted from the FID; rather, the water and the metabolite signals are simultaneously modeled during the calculations that generate the joint posterior probability for all of the parameters. While the McMC simulation algorithms are computationally intensive, they are well-suited for parallel computing systems. Computational times are highly dependent on the number of resonances in the model and the degree of resonance overlap. A typical analysis requires 5 to 30 minutes using a Dell PowerEdge 6850 server equipped with four Dual Core 3 GHz Xeon processors.

Acknowledgments

We thank Dr. Jeffrey J. Neil for his encouragement, support, and comments. We gratefully acknowledge our colleagues S.-W. Sun and S.-K. Song for providing the mouse-brain MRS data. This work was supported by the Small Animal Imaging Resource Program (SAIRP) of the National Cancer Institute, grant U24 CA83060; and by the National Institute of Neurological Disorders and Stroke, grants NS35912 and NS41519.

APPENDIX

In this appendix, we present the details of the Bayesian calculations. First, the model, Eq. [10], is split into models for the real and imaginary data channels (i.e., the 0°and 90° phase-shifted, timedomain data channels) and organized in a way that facilitates the marginalization of the amplitudes and constants. If dRi and dIi represent the real and imaginary parts of the complex data value sampled at time ti, then the model of the complex data is given by

dRi=MR(ti)+nRi
[11]

dIi=MI(ti)+nIi
[12]

where nRi and nIi represent the noise in the real and imaginary channels. The functions MR (ti) and MI (ti) are the real and imaginary parts of Eq. [10].

Bayesian probability theory is used to compute the joint posterior probability for each nonlinear parameter in the model, Eq. [10]. The Bayesian calculations use a Markov chain Monte Carlo (McMC) simulation to approximate the joint posterior probability distribution of the nonlinear parameters, independent of the amplitudes in Eqs. [10]. The target distribution of the McMC simulation is the joint posterior probability for the nonlinear parameters given all the real and imaginary data, DR and DI, so DR [equivalent] {dR1dRN} and DI [equivalent] {dI1dIN}. This probability is represented symbolically by P[mid ] DRDII), where Ω stands for all of the nonlinear parameters:

Ω{fM,αM,fW,αW,nc,ns,θ,t0}
[13]

where fM and αM are all of the metabolite resonance frequencies and decay rate constants, so fM [equivalent] {fM1fMm}, etc. The symbol I serves as a reminder that all probabilities are conditional on the prior information I.

The posterior probability for the nonlinear parameters is computed by application of Bayes’ theorem

P(ΩDRDII)=P(ΩI)P(DRDIΩI)P(DRDII)
[14]

where P[mid ] I) is the joint prior probability for the parameters, P (DRDI [mid ] ΩI) is the direct probability for the data given the nonlinear parameters, and P (DRDI [mid ] I) is a normalization constant. Dropping this normalization constant, one obtains

P(ΩDRDII)P(ΩI)P(DRDIΩI).
[15]

The direct probability for the data, P (DRDI [mid ] ΩI), is a marginal probability because it does not depend on any of the amplitudes or noise standard deviations. Reintroducing these parameters, and using the sum rule, one obtains

P(ΩDRDI)P(ΩI)dAdσP(DRDIAσΩI),
[16]

where A is the collection of all the amplitudes appearing in Eq. [10]. If Am represents the amplitudes of all the metabolite resonances, then Am [equivalent] {Am1Amm}. Similarly, if AC and AS represent all of the amplitudes in the cosine and sine expansions respectively, then, AC [equivalent] {AC1ACmC}.and AS [equivalent] { AS1ASmS}. The amplitudes associated with the constant model are designated Aε [equivalent] { ARFP, AIFP, ARC, AIC}, where ARFP is the real first-point amplitude, AIFP is the imaginary first-point amplitude, and ARC and AIC are the real and imaginary constant offsets respectively. The total number of amplitudes, nε , is given by

nε=m+nC+nS+4,
[17]

where m is the number of metabolite resonances, nC and nS are the number of the polynomials in the cosine and sine expansions Eqs. [8,9], respectively, and there are four constant amplitudes.

The right-hand side of equation [16] may be further expanded by repeated application of the sum rule,

P(ΩDRDI)P(θI)P(t0I)P(fWI)(αWI)P(ncI)P(nsI)×k=1m[P(fMkI)P(αMkI)]×dAdσ[P(σI)j=1nεP(AjI)P(DRΩAσI)P(DIΩAσI)]
[18]

where probabilities of the form, P[mid ] I) are prior probabilities for the respective parameters1, Aj represents the jth amplitude, and A represents all of the amplitudes appearing in the model. The two terms that are not of this form, P(DR [mid ] ΩAσI) and P(DR [mid ] ΩAσI), are direct probabilities for the real and imaginary data and are essentially likelihood functions. To make this factorization, P(DRDI [mid ] Ω AσI) = P(DR [mid ] Ω AσI)P(DI [mid ] ΩAσI), one must assume logical independence, i.e., the probability we assign for the real data does not depend on the imaginary data.

Next, numerical values are assigned to represent the various probabilities in Eq. [18]. We assign them in order, starting with the prior probability for the phase:

P(θI)={1360if0θ3600otherwise.
[19]

The prior probability for the delay time is assigned as a Gaussian:

P(t0I)=(2πσt02)12exp{t022σt02},
[20]

where σt0 = 3ΔT and ΔT is the sampling rate. This prior effectively assumes the real, but unknown, value of [mid ] t0 [mid ] < [mid ] 10ΔT[mid ]. Values outside this range are allowed, but highly improbable.

All of the prior probabilities for frequencies, including the water frequency, are assigned using the same functional form. For example, if γ is a resonance frequency, its prior probability is assigned from a low and high parameter value. If Lγ and Hγ represent the low and high bounds on γ, the prior probability is assigned as

P(γI){exp{(γ^γ)22σγ2}ifLγγHγ0otherwise
[21]

where

γ^(Lγ+Hγ)2
[22]

and

σγ(HγLγ)3.
[23]

So the entire low-high range represents a three standard deviation interval, with the prior tending to keep the frequency in the middle of its known range. Additionally, we specify the condition that the metabolite frequencies must be ordered so fm1 < fm2… < fmN.

The prior probability for all the decay rate constants is also assigned as a Gaussian. If δ represents one of these decay rate constants, then the prior probability assigned is:

P(δI){exp{δ22(3lb)2}if0δ10lb0otherwise
[24]

where lb (line broadening) is an average metabolite linewidth.

The prior probabilities for the polynomial expansion orders are assigned using a Laplacian:

P(ncI)exp{nc}
[25]

and

P(nsI)exp{ns}.
[26]

The prior probability for the standard deviation of the noise prior probability is assigned using a Jeffreys’ prior (35):

P(σI)1σ.
[27]

The prior probabilities for the amplitudes are assigned using a Gaussian of the form:

P(AjI)=(2πσ2β2Gjj)12exp{β2GjjAj22σ2}
[28]

where σ is the standard deviation of the noise prior probability. The quantity Gjj is the squared length of each complex model function considered as a vector and it ensures that each amplitude’s prior scales properly. The hyperparameter β is set to 0.01 in the program that implements this calculation.

The likelihood for the real part of the FID data is assigned as

P(DRΩAσI)(2πσ2)N2exp{i=1N[dR(ti)MR(ti)]22σ2}
[29]

where N is the number of complex data values in the FID. Similarly, the likelihood for the imaginary data is assigned as

P(DIΩAσI)(2πσ2)N2exp{i=1N[dI(ti)MI(ti)]22σ2}.
[30]

Assignment of all the probabilities in Eq. [20] is now complete.

Equations [19]-[30] are substituted into the joint posterior probabilities for the non-linear parameters in Eq. [18] and the integrals over the amplitudes and the noise standard deviation are evaluated analytically. The integrals over the amplitudes are Gaussian quadrature integrals. The integrals over the noise standard deviation is a Gamma function. Examples of how these integrals were evaluated are provided in (36) and (37).

Footnotes

1In parameter estimation, the priors provide bounds on the various estimated parameters, supplying reasonable ranges over which the Markov chains can search. In addition, the priors serve to define the labels (e.g., ‘1’, ‘2’) associated with each parameter (e.g., frequencies). For high signal-to-noise data, priors have little additional influence on the analysis beyond supplying reasonable bounds for the Markov chains and assuring consistent parameter labeling. In model selection, the priors serve an additional function, since the probability for a given model is a function of the product of likelihood and priors. Thus, each prior assignment acts to lower the posterior probability for a model, thereby serving as a penalty function against more complicated models.

References

1. McLean MA, Simister RJ, Barker GJ, Duncan JS. Magnetization transfer effect on human brain metabolites and macromolecules. Magn Reson Med. 2005;54(5):1281–1285. [PubMed]
2. Leibfritz D, Dreher W. Magnetization transfer MRS. NMR Biomed. 2001;14(2):65–76. [PubMed]
3. Haase A, Frahm J, Hanicke W, Matthaei D. 1H NMR chemical shift selective (CHESS) imaging. Phys Med Biol. 1985;30(4):341–344. [PubMed]
4. Ogg RJ, Kingsley PB, Taylor JS. WET, a T1- and B1-insensitive water-suppression method for in vivo localized 1H NMR spectroscopy. J Magn Reson, Series B. 1994;104(1):1–10. [PubMed]
5. Dreher W, Leibfritz D. New method for the simultaneous detection of metabolites and water in localized in vivo1H nuclear magnetic resonance spectroscopy. Magn Reson Med. 2005;54:190–195. [PubMed]
6. Hurd RE, Gurr D, Sailasuta N. Proton spectroscopy without water suppression: the oversampled J-resolved experiment. Magn Reson Med. 1998;40(3):343–347. [PubMed]
7. Clayton DB, Elliott MA, Lenkinski RE. In vivo proton spectroscopy without solvent suppression. Concepts Magn Reson. 2001;13:260–275.
8. van Der Veen JW, Weinberger DR, Tedeschi G, Frank JA, Duyn JH. Proton MR spectroscopic imaging without water suppression. Radiology. 2000;217(1):296–300. [PubMed]
9. Dong Z, Dreher W, Leibfritz D. Toward quantitative short-echo-time in vivo proton MR spectroscopy without water suppression. Magn Reson Med. 2006;55(6):1441–1446. [PubMed]
10. Zhu G, Smith D, Hua Y. Post-acquisition solvent suppression by singular-value decomposition. J Magn Reson. 1997;124(1):286–289. [PubMed]
11. Mierisova S, Ala-Kopela M. MR spectroscopy quantitation: a review of frequency domain methods. NMR Biomed. 2001;14:247–259. [PubMed]
12. Vanhamme L, Sundin T, Van Hecke P, Van Huffel S. MR spectroscopy quantitation: a review of time-domain methods. NMR Biomed. 2001;14:233–246. [PubMed]
13. de Beer R, van Ormondt D. Analysis of NMR data using time-domain fitting procedures. In: Rudin M, editor. In-vivo Magnetic Resonance Spectroscopy I: Probeheads, Radiofrequency Pulses, Spectrum Analysis, NMR Basic Principles and Progress series. Vol. 26. Berlin, Heidelberg: Springer-Verlag; 1992. pp. 201–248.
14. Pijnappel WWF, van den Boogaart A, de Beer R, van Ormondt D. SVD-based quantification of magnetic resonance signals. Journal of Magnetic Resonance. 1992;97(1):122–134.
15. Poullet J-B, Sima DM, Van Huffel S, Van Hecke P. Frequency-selective quantitation of short-echo time 1H magnetic resonance spectra. Journal of Magnetic Resonance. 2007;186(2):293–304. [PubMed]
16. Poullet J-B, Sima DM, Simonetti AW, De Neuter B, Vanhamme L, Lemmerling P, Van Huffel S. An automated quantitation of short echo time MRS spectra in an open source software environment: AQSES. NMR in Biomed. 2007;20(5):493–504. [PubMed]
17. Zhengchao D, Dreher W, Leibfritz D. Toward quantitative short-echo-time in vivo proton MR spectroscopy without water suppression. Magn Reson Med. 2006;55(6):1441–1446. [PubMed]
18. Barache D, Antoine J-P, Dereppe J-M. The Continuous Wavelet Transform, an Analysis Tool for NMR Spectroscopy. J Magn Reson. 1997;128(1):1–11.
19. Serrai H, Senhadji L, Clayton D, De Certaines JD. Evaluation of the continuous wavelet transform as a potential method for postacquisition water suppression in biomedical proton magnetic resonance spectroscopy. Spectr Lett. 2000;33(1):47–67.
20. Antoine JP, Coron A, Dereppe JM. Water peak suppression: time-frequency vs time-scale approach. J Magn Reson. 2000;144(2):189–194. [PubMed]
21. Cabanes E, Confort-Gouny S, Le Fur Y, Simond G, Cozzone PJ. Optimization of residual water signal removal by HLSVD on simulated short echo time proton MR spectra of the human brain. J Magn Reson. 2001;150(2):116–125. [PubMed]
22. Whittenburg SL. Bayesian analysis of time-domain Voigt functions. Spectrochimica Acta, Part A: Molecular and Biomolecular Spectroscopy. 1996;52A(10):1169–1174.
23. Abragam A. Principles of nuclear magnetism. Chapter 2 Oxford, UK: Clarendon Press; 1961.
24. Garwood M, DelaBarre L. The return of the frequency sweep: designing adiabatic pulses for contemporary NMR. J Magn Reson. 2001;153(2):155–177. [PubMed]
25. Tong Z, Yamaki T, Harada K, Houkin K. In vivo quantification of the metabolites in normal brain and brain tumors by proton MR spectroscopy using water as an internal standard. Magn Reson Imaging. 2004;22(7):1017–1024. [PubMed]
26. Simmons A, Smail M, Moore E, Williams SC. 1998 Serial precision of metabolite peak area ratios and water referenced metabolite peak areas in proton MR spectroscopy of the human brain. Magn Reson Imaging. 16(3):319–30. [PubMed]
27. Jenkins BG, Klivenyi P, Kustermann E, Andreassen OA, Ferrante RJ, Rosen BR, Beal MF. Nonlinear Decrease over time in N-acetyl aspartate levels in the absence of neuronal loss and increases in glutamine and glucose in transgenic Huntington’s disease mice. J Neurochem. 2000;74(5):2108–2119. [PubMed]
28. Choi CB, Kim HY, Han DY, Kang YW, Han YM, Jeun SS, Choe BY. In vivo1H MR spectroscopic findings in traumatic contusion of ICR mouse brain induced by fluid percussion injury. Eur J Radiol. 2005;55(1):96–101. [PubMed]
29. Schwarcz A, Natt O, Watanabe T, Boretius S, Frahm J, Michaelis T. Localized proton MRS of cerebral metabolite profiles in different mouse strains. Magn Reson Med. 2003;49(5):822–827. [PubMed]
30. Ogbron MR, Sareen S, Prychitko J, Buist R, Peeling J. Altered organic anion and osmolyte content and excretion in rat polycystic kidney disease: an NMR study. Am J Physiol-Renal Physiol. 1997;272(1):F63–69. [PubMed]
31. Dixon RM, Frahm J. Localized proton MR spectroscopy of the human kidney in vivo by means of short echo time STEAM sequences. Magn Reson Med. 1994;31(5):482–7. [PubMed]
32. Avison MJ, Rothman DL, Nixon TW, Long WS, Siegel NJ. 1H NMR study of renal trimethylamine responses to dehydration and acute volume loading in man. Proc Natl Acad Sci USA. 1991;88(14):6053–6057. [PubMed]
33. Clayton DB, Elliott MA, Leigh JS, Lenkinski RE. 1H spectroscopy without solvent suppression: characterization of signal modulations at short echo times. J Magn Reson. 2001;153(2):203–209. [PubMed]
34. Serrai H, Clayton DB, Senhadji L, Zuo C, Lenkinski RE. Localized proton spectroscopy without water suppression: removal of gradient induced frequency modulations by modulus signal selection. J of Magn Reson. 2002;154(1):53–59. [PubMed]
35. Jefferys H. Theory of probability. Chapter 3 Oxford, UK: Clarendon Press; 1983.
36. Sivia DS. Data analysis: a Bayesian tutorial. Oxford, UK: Oxford University Press; 1966.
37. Bretthorst GL. Bayesian spectrum analysis and parameter estimation. In: Berger J, editor. Lecture notes in statistics. Vol. 48. New York: Springer-Verlag; 1988.