Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2276337

Formats

Article sections

Authors

Related links

J Magn Reson. Author manuscript; available in PMC 2008 March 28.

Published in final edited form as:

Published online 2005 December 1. doi: 10.1016/j.jmr.2005.11.004

PMCID: PMC2276337

NIHMSID: NIHMS17871

* Corresponding author. Fax: +1 410 614 1977. E-mail address: ude.uhj.irm@lmottob (P.A. Bottomley)

The publisher's final edited version of this article is available at J Magn Reson

See other articles in PMC that cite the published article.

Accurate and robust quantification of in vivo magnetic resonance spectroscopy (MRS) data is essential to its application in research and medicine. The performance of existing analysis methods is problematic for in vivo studies where low signal-to-noise ratio, overlapping peaks and intense artefacts are endemic. Here, a new frequency-domain technique for MRS data analysis is introduced wherein the circular trajectories which result when spectral peaks are projected onto the complex plane, are fitted with active circle models. The use of active contour strategies naturally allows incorporation of prior knowledge as constraint energy terms. The problem of phasing spectra is eliminated, and baseline artefacts are dealt with using active contours-snakes. The stability and accuracy of the new technique, *CFIT*, is compared with a standard time-domain fitting tool, using simulated ^{31}P data with varying amounts of noise and 98 real human chest and heart ^{31}P MRS data sets. The real data were also analyzed by our standard frequency-domain absorption-mode technique. On the real data, *CFIT* demonstrated the least fitting failures of all methods and an accuracy similar to the latter method, with both these techniques outperforming the time-domain approach. Contrasting results from simulations argue that performance relative to Cramer-Rao Bounds may not be a suitable indicator of fitting performance with typical in vivo data such as these. We conclude that *CFIT* is a stable, accurate alternative to the best existing methods of fitting in vivo data.

Accurate quantification of in vivo MR spectroscopic (MRS) data is often thwarted by low signal-to-noise ratio (SNR), overlapping resonances, and baseline artefacts arising from physiological motion, localizing gradients, and/or delayed acquisition. Standard quantification approaches apply numerical methods to estimate parameters that characterize a mathematical model of the MRS signal, and are generally classified as either time- or frequency-domain methods. Time-domain methods operate directly on the raw data and offer some convenient ways of addressing incomplete or corrupted data [1]. They can be further divided into interactive and black box methods. Interactive methods employ numerical optimization techniques to iteratively solve a non-linear least-squares optimization problem [1–3] based on user-defined initial parameter values. Prior knowledge about peak structure may be used to constrain the solution. The optimization can be performed either locally, or globally [4,5]. Global optimization does not require user interaction but is inherently computationally inefficient. Black box methods provide non-iterative fitting solutions employing linear prediction or state-space theory [1,3,6–8]. These are fast and require little user interaction, but perform poorly with in vivo data because they allow little or no prior knowledge. A major problem with time-domain techniques is that frequency-selective analysis is not straightforward, albeit soluble [9–12]. The frequency domain is certainly the natural domain for spectral analysis: it is in the frequency domain wherein user interactivity invariably occurs, from the input of initial values and prior knowledge, to the final judgment of acceptance of a fit.

Frequency-domain methods model the spectrum as a set of Lorentzian, Gaussian or Voigt peaks, and fit the data using least-squares methods. The real and imaginary parts or the complex signal [13,14] can be used in the fitting process, but quantification is more robust when the number of parameters being estimated is minimized. This can be achieved, for example, by first correcting the spectrum for zero and first-order phase, so that all of the peaks have zero phase. Then, only the pure absorption-mode signal need be fitted. However, the baseline is often severely distorted by the acquisition delays required for spatial encoding gradients which convolve the spectrum with sinc functions [15], by transient eddy currents from the gradients, and by physiological motion. This too must be corrected, either by modeling the baseline with a set of polynomial or sinusoidal basis functions, or by deconvolution.

While new solutions are being investigated [16–18], phase and baseline corrections are difficult to automate and usually cannot match the performance of an experienced manual operator when artefacts are large. This is evidenced, for example, by a 50% increase in the standard deviation of normal human cardiac phosphocreatine (PCr) to adenosine triphosphate (ATP) ratios quantified in phosphorus (^{31}P) spectra by automated Marquardt frequency-domain fitting and by non-linear least squares fitting in the time-domain [19], as compared to manual frequency-domain fitting of data acquired with the identical protocol and equipment [20]. A 50% increase in the uncertainty of a PCr/ATP measurement is equivalent to a 50% reduction in the SNR of the measurement from what could have been achieved. Frequency-domain techniques have been reviewed elsewhere [13], and other variants include the LCModel which fits the spectrum with a linear combination of in vitro spectra acquired from individual metabolites [21,22], and time-frequency analysis using wavelets which separates and subtracts each component from the raw signal [23].

In this paper, we present a new frequency-domain technique for MRS quantification that overcomes phasing problems and addresses baseline issues with current techniques as applied to low-SNR, high-artifact regimes. Using a Lorentzian model, we show that the real-imaginary trajectory of each peak in a spectrum scribes a circle in the complex plane with a diameter equal to the peak height. The complex spectrum is fitted using a set of active circle models. The use of both real and imaginary components as a continuous circular trajectory with this circle fitting (CFIT) approach improves the stability of the estimation process, effectively eliminates zero and first-order phase-correction, and reduces the baseline problem to one of completing circles for peaks whose loops are open in the complex plane. Prior knowledge is easily incorporated as external energy terms [24]. The CFIT technique is tested with simulated brain ^{31}P spectra, as well as with real in vivo localized human chest and heart ^{31}P spectra in which problems of low SNR, phasing, baseline, and overlapping peaks are endemic. The stability and accuracy of CFIT are compared with our current standard analysis tool, a frequency-domain absorption mode fitting routine, *CSX* [25,26], and with a widely available time-domain method, *AMARES* [6].

The complex time-domain MRS signal, *s*(*t*), is most commonly modeled as the sum of exponentially damped sinusoids:

$$s(t)=\sum _{k=1}^{K}{a}_{k}{\text{e}}^{j{\phi}_{k}}{\text{e}}^{(-{d}_{k}-j2\pi {f}_{k})t}+n(t),$$

(1)

where *K* is the number of sinusoidal components in the signal, *f _{k}* is the frequency of the

$$S(f)=\sum _{k=1}^{K}\frac{{a}_{k}}{{d}_{k}}{\text{e}}^{j{\phi}_{k}}(\frac{1}{1+{(2\pi (f-{f}_{k})/{d}_{k})}^{2}}+j\frac{2\pi (f-{f}_{k})/{d}_{k}}{1+{(2\pi (f-{f}_{k})/{d}_{k})}^{2}})+N(f).$$

(2)

Eq. (2) describes a helix in the three-dimensional (3D) space comprised of the real-imaginary complex plane and the chemical shift/frequency axis. The bracketed terms are of the form, *x* = 1/(1 + *t*^{2}), *y* = t/(1 + *t*^{2}), which is a parametric equation for a circle. The projection of the helix onto the complex plane thus forms a set of parametric equations that, for Lorentzian lines, describe circles with centers at the real-imaginary coordinates (*a _{k}*/2

The spectrum of a single peak is a helix in the 3D space comprised of the real, the imaginary, and the frequency axes. The projection of the helix onto the complex plane is a circle (A), while the projection of the helix onto the plane of the real and **...**

(A) The real part of the simulated ^{31}P human brain spectrum with noise *σ* = 25. Part (B) is the projection of the 3D spectrum trajectory onto the real-imaginary plane. Part (C) is a 3D plot of the spectrum showing the helical trajectory of the **...**

Because the trajectory of each peak in the complex plane is circular, we can model them with active circles. The active circles are contours that adaptively deform to best fit the model spectrum to the measured spectrum by minimizing a measure of the fit error while preserving their circular shape in the complex plane. The fit error, *e*(*f*), is the complex difference between the measured spectrum and the sum of the model defined in Eq. (2) plus the baseline in the complex plane. Initial values for the parameters are either provided by the user, estimated from the spectrum, or set to fixed values. The peak center frequency, *f _{k}*, needs to be set reasonably accurately, especially when the peaks partially overlap as is typically the case with

$$E=\int e{(f)}^{\ast}e(f)\phantom{\rule{4pt}{0ex}}\text{d}f$$

(3)

using a gradient descent approach leads to an iterative solution in which the parameters change under the effect of the fit error. In this sense, *e*(*f*) acts as an external error *force* [24].

The error force acting on the circumference of an active circle can be resolved into a radial-like component emanating from the origin of the complex plane, *e _{rk}*(

We now describe how the model circles evolve under the effect of the error forces. Details of how the minimization of E leads to the following results are given in Appendix A. Consider the *k*th circle, the incremental changes in the circle radius *r _{k}*, the center frequency

$$\mathrm{\Delta}{r}_{k}={c}_{r}\int {w}_{rk}(f-{f}_{k}){e}_{rk}(f)\text{d}f,$$

(4)

$$\mathrm{\Delta}{f}_{k}={c}_{f}\int {w}_{fk}(f-{f}_{k}){e}_{tk}(f)\text{d}f,$$

(5)

$$\mathrm{\Delta}{\phi}_{k}={c}_{\phi}\int {w}_{\phi k}(f-{f}_{k}){e}_{ok}(f)\text{d}f,$$

(6)

$$\mathrm{\Delta}{d}_{k}={c}_{d}\int {w}_{dk}(f-{f}_{k}){e}_{tk}(f)sgn(f-{f}_{k})\text{d}f,$$

(7)

where *e _{rk}*(

Baseline distortions manifest as a displacement of points on the circle such that the ends of each peak no longer complete a circle in the complex plane. To account for baseline distortion, the baseline is modeled as a free-form active contour: a snake [24]. The snake does not take the form of any specific function such as a polynomial or the like (Appendix A). While constrained by its elastic properties, the snake is adjusted to minimize an energy term representing the residuals of the fit between the measured data and the model spectrum. The elastic properties of the contour are classically defined in terms of its length and curvature. The baseline is simultaneously fitted with the peaks by updating both the active circles and the baseline in each iteration of the fitting process.

The CFIT technique can thus be summarized by the following steps:

- Initialize
*K*active circles by peak picking or by a black-box method; - Calculate the error force between the measured data and the model spectrum;
- Adjust the baseline snake parameters to minimize the fit error;
- Repeat steps 2–4 until either equilibrium or the maximum number of iterations is reached;
- Output peak frequencies, areas, and heights. Display the fit.

CFIT is applied to data prepared in the frequency-domain in the usual manner with exponential filtering, zero-filling (ZF), and FT of the time-domain signal [27]. The spectrum is typically normalized by the maximum value of the magnitude of a standard spectrum so that the same set of parameters, *c _{r}*,

Results of quantifying the simulated ^{31}P spectra comprised of 256 complex points using CFIT with an 8-fold (ZF = 8) and a 4-fold zero-filling (ZF = 4), and AMARES. The success rate (A), the frequency RMSE (B), the amplitude RRMSE (C), and the damping **...**

Although phasing the spectrum is not required in CFIT, a basic zero- and first-order phasing function is included for the sole purpose of facilitating a user’s ability to interactively pick peak parameters. This is standard for most current user-interactive approaches. The fitting process actually uses the initial un-phased spectrum to avoid any added distortions that may arise from first-order phase correction [15]. Phasing may also be applied to the post-processed fitted spectrum, if desired, for display purposes.

Prior knowledge is easily incorporated in CFIT by assigning an energy term to each constraint or feature that one wishes the solution to stick to. The constraints need not be linear. Minimization of the constraint energy is also carried out using the energy descent method. The external constraints affecting the model spectrum can be represented by spring forces that drive the contour towards the features of interest. Each of the active circles adapts under the influence of the error forces and the constraints until it is stabilized. In general, a parameter *p* that is to be estimated, is adjusted in each iteration under the influence of the two sets of forces via

$$p\leftarrow p+{k}_{p}({p}_{d}-p)+\int {e}_{p}(f)\text{d}f,$$

(8)

where *p _{d}* is the desired value that

CFIT was first tested with a simulated in vivo ^{31}P human brain spectrum consisting of 256 complex time-domain data points, zero-filled eight-times and comprised of 11 exponentials, defined elsewhere [6] with parameters summarized in Table 2. The spectrum in Fig. 2 is illustrated with complex white Gaussian noise of standard deviation (SD) *σ* = 25 added, and a 10-Hz exponential filter applied. The sampling frequency, *f _{s}* is 3 kHz, and the fine peaks in the ATP multiplets are split by 16 Hz. The technique was tested on 100 different simulated signal-plus-noise realizations, for each of 20 SNR levels in the range 5 ≤

Prior knowledge of the ATP doublets (α-ATP and γ-ATP peaks) and triplet (β-ATP) was implemented via the following constraints; *d*_{1} = *d*_{2} = *d*_{3}, *d*_{4} = *d*_{5}, *d*_{6} = *d*_{7}, *a*_{1} = *a*_{2}/2 = *a*_{3}, *a*_{4} = *a*_{5}, *a*_{6} = *a*_{7}, _{1} = _{2} = _{3}, _{4} = _{5}, _{6} = _{7}, *f*_{3} = *f*_{2} + 16 = *f*_{1} + 32, *f*_{5} = *f*_{4} + 16, *and f*_{7} = *f*_{6} + 16. The parameters *c _{r}*,

For comparison, the same data was quantified using AMARES as provided in the standard jMRUI package, available from http://www.mrui.uab.es/mrui/upon request. AMARES [6] is a time-domain fitting technique that has been demonstrated to have stable and accurate performance [6,8]. This program requires the user to pick the peaks providing initial values for the peak frequencies and linewidths, and to set the constraints in a similar fashion as in CFIT. The constraints used for AMARES (*d _{k}*,

As reported in [8], the robustness of each method is evaluated by computing its success rate, i.e., the number of times, out of the total number of simulation runs, that the method is able to resolve the 11 peaks within specific intervals lying symmetrically around the true frequencies of the peaks. The half-widths of these intervals are set to 8 Hz, i.e., half the separation of the closest peaks in the data. The peak error is quantified by the mean RRMSE defined as

$$\text{meanRRMSE}=\frac{100}{K}\sum _{k=1}^{K}\sqrt{\frac{1}{J}\sum _{j=1}^{J}\frac{{({p}_{kj}-{p}_{k})}^{2}}{{p}_{k}^{2}}}$$

(9)

for the amplitude and damping estimates, and as the mean RMSE

$$\text{mean\hspace{0.28em}RMSE}=\frac{1}{K}\sum _{k=1}^{K}\sqrt{\frac{1}{J}\sum _{j=1}^{J}{({p}_{kj}-{p}_{k})}^{2}}$$

(10)

for the frequency estimates, where *K* = 11, is the total number of peaks, *J* is the number of simulation runs in which the method was able to find every peak within the corresponding frequency interval, *p _{kj}* denotes the parameter estimate of the

Fig. 3 shows that the success rate for CFIT with high zero-filling and AMARES using simulated data are comparable. Both perform perfectly when SNR is high (*σ* < 50), but deteriorate as SNR declines further. Figs. 3(B–D) plots the mean RRMSE values for the amplitude and the damping constant, and the mean RMSE values for the frequency of the ATP complex (the first seven peaks in the simulated spectrum) as a function of noise level for the simulated data. The Cramer-Rao bound (CRB), a measure of the best attainable precision of an unbiased estimate in the presence of random noise [28], is also plotted for each of the estimated parameters. The results of these simulations are consistent with AMARES’ reputation as one of the best performing methods relative to the CRB [8]. Even so, relative performance can reverse in the presence of baseline artifacts. For example, perturbing the first four time points of the model brain spectrum by multiplication by random numbers produces a mild baseline roll as shown in Fig. 4. The success rate of CFIT is comparable to the rate without perturbation, albeit with some increase in the amplitude error. However the performance for uncompensated AMARES is decimated. AMARES performance does improve with initial filtering to reverse the effect of the imposed distortion, but CFIT remains the more successful in this instance (Fig. 4A). Fig. 5 shows a sample output of CFIT for one of the processed simulated spectrum for the case *σ* = 25. At this level of noise, the SNR of the β-ATP peak is about 25.

CFIT and AMARES were tested on 98 in vivo localized ^{31}P spectra acquired in 5–7 min from the human chest and heart at 1.5 T using a cardiac-gated, quantitative, one-dimensional chemical shift imaging (CSI) protocol employing adiabatic 60° excitations and an 0.7 ms acquisition delay as described elsewhere [29]. The spectra derive from normal subjects (30%), patients with dilated cardiomyopathy and congestive heart failure (40%), and patients with left ventricular hypertrophy (30%), participating in ongoing studies of cardiac creatine kinase flux wherein scan times are limited by the need to complete six MRS acquisitions in a single exam [29,30]. Despite its performance in simulations with respect to the CRB, AMARES is not used for in vivo cardiac ^{31}P MRS quantification in our laboratory due to an unacceptable failure rate. Instead, we routinely use a frequency-domain interactive fitting program, CSX [25,26], wherein peak height, frequency, and linewidths are fitted by a simplex routine after the user phases the spectrum, selects points for baseline correction, and identifies initial frequencies of peaks selected for fitting. The performance of CFIT was therefore additionally compared to CSX, our reference standard, for which fit results for all 98 spectra had already been documented. Spectra were thus selected for analysis based only on the chronological order in which studies were performed, and the availability of pre-existing CSX measurements with a minimum fitting threshold for SNR of about 5 for β-ATP.

Spectra to be analyzed for CFIT were prepared with 8-fold zero-filling, and a 10-Hz exponential filter. The adapting parameters and the windowing function were the same as in the simulations. The spectra were roughly phased by an experienced in vivo spectroscopist (PAB) who identified the initial peak locations. As in the simulations, the β-ATP peak was fitted to a triplet and the α- and γ-ATP peaks were fitted to doublets. One-to-three Lorentzian peaks were used to model the PCr peak since its shape was not Lorentzian: the use of multiple Lorentzian peaks to approximate other non-Lorentzian peak shapes (Gaussian, Voigt, etc) has been reported before [7]. The multiple PCr peaks are linked so that they have the same phase. This constraint is essential for eliminating solutions which have peaks that cancel. The in vivo spectra were quantified with AMARES in the same way by the same user. AMARES requires additional prior knowledge estimates of linewidths for the ATP multiplets, which are provided interactively. The spectra were analyzed 0.5–3 years earlier in the same way by the same user with the CSX program, except that they were more carefully phased (zero and first-order), baselines were corrected with user-defined points between the peaks, and the ATP peaks were fitted only as singlets but with Gaussian-shaped lines. These differences reduced the number of parameters that needed to be estimated for the CSX method compared to CFIT and AMARES.

Because the true values of in vivo parameters are not known, fitting failures cannot in practice be defined by comparison with true values. However, we can record a failure when the fit misses one peak completely, or when anomalously small or anomalously large areas of a peak are returned. The number of fit failures involving one or more of the ATP or PCr peaks in the ^{31}P MRS cardiac data was 4/98 for CFIT. This compares with 29/98 for AMARES, and 8/98 failures for CSX, the latter involving ATP only. The total number of cases in which at least one of the three fitting methods fails was 31.

The ratio of the peak areas of the β-ATP and γ-ATP provides an independent gauge of the robustness of the different fitting routines. This ratio must be approximately unity since the nucleoside triphosphate sources are the same, and additional contributions to the γ-resonance from other phosphates are small. The ratio obtained using the three fitting methods is plotted in Fig. 6, and the means and SD for the ratio obtained by each method are listed in Table 3. While AMARES yields erroneous results that greatly deviate from unity in many cases, CFIT performs comparably with CSX, our current standard.

The ratio of the peak areas β-ATP/γ-ATP for the in vivo ^{31}P spectra from the human chest and heart, as estimated using CSX, AMARES, and CFIT. The plot excludes all of the 31 fit failures. The true value for the ratio is ~1.

Mean and SD of the ratio of the peak area of the β-ATP/γ-ATP estimated using different methods for the in vivo ^{31}P MRS chest and heart data

The PCr/ATP ratio of peak areas is a physiologically important metabolic measure [20]. Fig. 7 plots the PCr/β-ATP ratio measured by AMARES and CFIT against our current standard, CSX. With perfect agreement, all the data would fall on the 45° line and lie in the known range of about 0.2–5 for normal and sick human chest and heart. The CSX data show reasonable agreement with CFIT but much larger disparity for AMARES, with values lying outside the physiological range. Table 4 reports the Pearson correlation coefficient analysis of the PCr/β-ATP, PCr/α-ATP, and PCr/γ-ATP ratios measured using AMARES, CSX, and CFIT. To avoid degrading the regression with the much larger number of outliers produced by AMARES, PCr/ATP ratios greater than 10 were eliminated from this analysis. Clearly, the agreement between CSX and CFIT is the highest for all three ratios. Again even with all of the AMARES fit failures excluded, assuming CSX *is* our standard, the mean RRMSE of the peak areas of β-ATP, α-ATP, and PCr peaks show greater scatter when measured with AMARES than with CFIT, especially for PCr, as demonstrated in Table 5. A sample CFIT analysis of one of the processed spectra is shown in Fig. 8.

The ratio of the peak areas PCr/β-ATP for the in vivo ^{31}P chest and heart spectra as estimated using AMARES (circles) and CFIT (points) vs. the CSX measurements (horizontal axis). AMARES and CFIT points that agree with CSX should fall on the 45° **...**

Analysis of a real ^{31}P spectrum from a female heart patient using CFIT. The SNR of the β-ATP peak is about 8. From bottom up, the real part of the measured spectrum, the real part of the fitted spectrum, the real part of the individual circles, **...**

Correlations (*r*) between AMARES, CSX, and CFIT for the peak area ratios PCr/β-ATP, PCr/α-ATP, and PCr/γ-ATP for the in vivo ^{31}P MRS signals

Mean RRMSE values and minimum (m) and maximum (M) error relative to CSX for the amplitude estimates obtained using AMARES and CFIT from the in vivo ^{31}P MRS signals after exclusion of the 31 fit failures

The computational time in CFIT is longer than AMARES but comparable or faster than the simple CSX frequency-domain fitting routine. Implementation on a PC (Intel Pentium M processor 1.5 GHz and 512 MB of memory) operating on the Microsoft Windows XP platform in MATLAB v6.5 (Mathworks, Natick, MA) required an average of 87 ± 47 iterations to achieve a stable solution in 2.9 ± 1.6 s.

We have introduced a new frequency-domain spectral analysis technique, CFIT, that utilizes *active circles* to directly fit the circular trajectory that Lorentzian resonances scribe in the complex plane. The method has the advantage of eliminating the effect of zero- and first-order phase correction on peak quantification. Baseline roll and overlapping peaks affect the continuity of circles in the complex plane, and are dealt with using *active snakes*. Prior knowledge is incorporated in the usual way, and non-Lorentzian peaks accommodated using multiple Lorentzian peaks. The ability of CFIT to accommodate large first-order phase errors and baseline artefacts suits it well to in vivo analyses where such problems are endemic. It is especially suited to quantifying low-SNR high-artifact free-induction decays in which acquisition is delayed by gradient-encoding, such as cardiac ^{31}P CSI data [29,30].

It is disappointing that the performance of CFIT is inferior to that of the standard time-domain fitting technique, AMARES, with simulated data, as indexed by the Cramer-Rao bounds for key fitting parameters in the absence of baseline artefacts. However, AMARES is not used in our laboratory for analyzing human cardiac ^{31}P CSI data due to its unacceptably high failure-to-fit rate, an anecdotal observation not previously documented. Therefore, the performance comparison of CFIT on about a hundred real in vivo data sets analyzed by our standard absorption-mode frequency-domain technique, CSX, as well as by AMARES is, in our view, the more relevant comparison. Despite the excellent performance of AMARES on simulated noisy data free of baseline artefacts, we document with simulations that the failure rate of AMARES can become excessive in the presence of distortions that cause baseline artifacts. Moreover, for the in vivo data its failure rate is indeed unacceptably high compared with the other two methods, and its accuracy much lower. Conceivably, better adjustment of AMARES constraints and filtering could enhance the performance with our in vivo data, just as better filtering that does not distort the circles could also benefit CFIT. In our view, the technique of choice is the one that provides the best fit to the worst spectrum from which measurements *must* be made, within the limits of detection permitted by the SNR. CFIT had the lowest failure rate with real human cardiac ^{31}P CSI data, and performed comparable to CSX with respect to measures of accuracy (β-ATP/γ-ATP; Fig. 6) and scatter (RRMSE; Tables 3–5, Figs. 6 and and7)7) on peaks with SNR as low as ~5. By these criteria CFIT is a suitable replacement for CSX, our current standard for cardiac ^{31}P MRS.

New quantification techniques are often accompanied by claims of accuracy and robustness [6,14]. However, caution is warranted when advantages are demonstrated only in cases of high SNR, idealized line-shapes, or signals containing only Gaussian noise, such as is common with simulations. In particular, many in vivo signals contain significant artefacts that cannot be modeled with Gaussian noise. An estimator that deals with in vivo data assuming Gaussian noise may thus be misled in situations where such artefacts arise. This observation together with our performance measurements on simulated and real data, render questionable the suitability of the Cramer-Rao bound as a measure of the performance of estimators of in vivo MRS signals in such settings.

Two key factors help CFIT overcome the effects of artefacts, resulting in better in vivo performance. The first is that the determination of step size during parameter adaptation involves an integration process, which tends to cancel random errors. The second factor is the bell-shaped window (implemented here as Lorentzian) used for weighting the error force that affects each active circle. Such weighting greatly reduces the contamination of nearby artefacts and also allows for selective frequency analysis. In this case, residuals from outside the quantified area are detected using the baseline snake. The ability to impose prior knowledge is an additional major factor that stabilizes the performance of fitting routines in general, including CFIT. Incorporating prior knowledge has been shown to benefit accuracy [28], and a wide range of prior knowledge including non-linear constraints, can be incorporated to guide quantification with CFIT as desired.

In conclusion, CFIT is a new frequency-domain spectral fitting technique that has been tested on real in vivo cardiac ^{31}P MRS data, warts-and-all. It provides accurate stable performance, and appears to be a suitable alternative to existing standards.

We thank Roger Eastman of Loyola College in Baltimore and Teresa Laudadio, Department of Electrical Engineering, ESAT-SCD(SISTA), Katholieke Universiteit Leuven, Leuven (Belgium) for help and discussions. The data derive from ongoing human heart studies performed with Robert G Weiss, for which Peter Barker’s CSX analysis tools have served us well and long: both are at Johns Hopkins.

Consider the model spectrum *S*(*f*) given in Eq. (2) including the baseline model. We define:

$${p}_{k}(f)=2\pi (f-{f}_{k})/{d}_{k}\phantom{\rule{12pt}{0ex}}\text{and}\phantom{\rule{12pt}{0ex}}{r}_{k}=\frac{{a}_{k}}{2{d}_{k}}.$$

(A.1)

Then

$$\begin{array}{ll}S(f)=& {\displaystyle \sum _{k=1}^{K}}2{r}_{k}{\text{e}}^{j{\phi}_{k}}\left(\frac{1}{1+{p}_{k}{(f)}^{2}}+j\frac{{p}_{k}(f)}{1+{p}_{k}{(f)}^{2}}\right)\\ & +\phantom{\rule{4pt}{0ex}}\text{baseline}(f)\end{array}$$

(A.2)

The fit error is the difference between the measured and the model spectra,

$$e(f)={S}_{\text{measured}}(f)-S(f),$$

(A.3)

with the error energy, *E*, given by Eq. (3). *E* is minimized using the gradient descent method.

The gradient of *E* with respect to *r _{k}*,

$$\frac{\partial E}{\partial {r}_{k}}\int \frac{\partial}{\partial {r}_{k}}e(f)\ast \phantom{\rule{12pt}{0ex}}e(f)\text{d}f$$

(A.4)

reduces to

$$\begin{array}{l}\frac{\partial E}{\partial {r}_{k}}=-2\int \langle e(f),{u}_{k}(f)/{r}_{k}\rangle \text{d}f\phantom{\rule{4pt}{0ex}}\text{with}\phantom{\rule{4pt}{0ex}}{u}_{k}(f)\\ \equiv 2{r}_{k}{\text{e}}^{j{\phi}_{k}}\left(\frac{1}{1+{p}_{k}{(f)}^{2}}+j\frac{{p}_{k}(f)}{1+{p}_{k}{(f)}^{2}}\right),\text{or}\end{array}$$

(A.5)

$$\frac{\partial E}{\partial {r}_{k}}=-4\int \langle \frac{1}{1+{p}_{k}{(f)}^{2}}e(f),{\widehat{u}}_{k}(f)\rangle \text{d}f,$$

(A.6)

where .,. is the inner product; *u _{k}* (

The gradient of *E* with respect to * _{k}*,

$$\frac{\partial E}{\partial {\phi}_{k}}=-4\int \langle \frac{{r}_{k}}{1+{p}_{k}{(f)}^{2}}e(f),{\stackrel{\u2322}{u}}_{k}^{\perp}(f)\rangle \text{d}f,$$

(A.7)

$$\frac{\partial E}{\partial {f}_{k}}=-8\pi \int \langle \frac{{r}_{k}/{d}_{k}}{1+{p}_{k}{(f)}^{2}}e\left(f),{\stackrel{\u2322}{\rho}}_{k}^{\perp}(f)\right)\rangle \text{d}f,$$

(A.8)

and

$$\frac{\partial E}{\partial {d}_{k}}=-8\pi \int \langle \frac{({r}_{k}/{d}_{k}){p}_{k}(f)}{1+{p}_{k}{(f)}^{2}}e(f),{\stackrel{\u2322}{\rho}}_{k}^{\perp}(f)\rangle \text{d}f,$$

(A.9)

where
${u}_{k}^{\perp}(f)\equiv j{u}_{k}(f)$ is perpendicular to the radial vector *u _{k}*(

We define the components of *e*(*f*) in the radial-like direction *u _{k}* (

$$\frac{\partial E}{\partial {r}_{k}}=-4\int {w}_{rk}(f){e}_{rk}(f)\text{d}f,$$

(A.10)

$$\frac{\partial E}{\partial {\phi}_{k}}=-4\int {w}_{\phi k}(f){e}_{ok}(f)\text{d}f,$$

(A.11)

$$\frac{\partial E}{\partial {f}_{k}}=-8\pi \int {w}_{fk}(f){e}_{tk}(f)\text{d}f,$$

(A.12)

and

$$\frac{\partial E}{\partial {d}_{k}}=-8\pi \int {w}_{dk}(f)sgn(f-{f}_{k}){e}_{tk}(f)\text{d}f,$$

(A.13)

where the weighting functions *w _{rk}* (

$${w}_{rk}(f)=\frac{1}{1+{p}_{k}{(f)}^{2}},$$

(A.14)

$${w}_{\phi k}(f)=\frac{{r}_{k}}{1+{p}_{k}{(f)}^{2}},$$

(A.15)

$${w}_{fk}(f)=\frac{{r}_{k}/{d}_{k}}{1+{p}_{k}{(f)}^{2}},$$

(A.16)

and

$${w}_{dk}(f)=\frac{({r}_{k}/{d}_{k})\cdot \mid {p}_{k}(f)\mid}{1+{p}_{k}{(f)}^{2}}.$$

(A.17)

These weighting functions are Lorentzian with the same linewidth as the processed peak. In our implementation, we used this same weighting function with half the line width because it yielded better accuracy and stability than with the full linewidth. Other weighting functions with different linewidths and properties can be substituted to achieve preferred levels of artefact rejection.

We refer to this as an *active circle model*. It can be considered an example of the more general set of analytical deformable models that are contours with a general predefined shape. In these models, prior knowledge about the shape is used directly by defining a curve described by a set of parameters that encode the general shape characteristics [31–33].

Active contour models, also known as snakes, are energy-minimizing contours that are constrained by their internal elastic forces, while driven by external forces towards desired external features. In image analysis applications, where active contour models were introduced [24], the external forces are provided by the image brightness and edge characteristics as well as higher-level information. The contour energy functional is constructed so that local minima correspond to desired features. The fitting problem is formulated as an energy minimization problem, and is solved using a suitable minimization technique. A contour is described parametrically by *v*(*s*) = (*x*(*s*), *y*(*s*)) where *x*(*s*), *y*(*s*) are *x* and *y* coordinates along the contour and *s* in [0,1] is the contour length. The contour model defines the energy of a contour *v*(*s*), the snake energy *E*_{snake}, as,

$${E}_{\text{snake}}={\int}_{0}^{1}{E}_{int}(v(s))+{E}_{\text{ext}}(v(s))+{E}_{\text{con}}(v(s))\text{d}s,$$

(A.18)

where *E*_{int} is the internal energy of the contour, imposing continuity and curvature constraints, *E*_{ext} is the external energy constructed to attract the snake to desired features; and *E*_{con} is another energy term which allows various constraints to be applied during the minimization process. An initial contour evolves by minimizing the snake energy using a minimization technique like gradient descent, as used here. In this work, the baseline is modeled as a snake constrained by its smoothness and is driven by the fit error. The baseline is simultaneously fitted with the peaks by performing the active circles updates and a baseline update in each iteration of the solution.

^{}NIH Grant 2RO1 HL 56882.

1. Vanhamme L, Sundin T, Hecke PV, Huffel SV. MR spectroscopy quantitation: a review of time-domain methods. NMR Biomed. 2001;14(4):233–246. [PubMed]

2. van der Veen JWC, de Beer R, Luyten PR, van Ormondt D. Accurate quantification of in vivo 31P NMR signals using the variable projection method and prior knowledge. Magn Reson Med. 1988;6(1):92–98. [PubMed]

3. Vanhamme L, Van Huffel S, Van Hecke P, van Ormondt D. Time domain quantification of series of biomedical magnetic resonance spectroscopy signals. J Magn Reson. 1999;140(1):120–130. [PubMed]

4. DiGennaro FS, Cowburn D. Parametric estimation of time-domain NMR signals using simulated annealing. J Magn Reson. 1992;96:582–588.

5. Metzger GJ, Patel M, Hu X. Application of genetic algorithms to spectral quantification. J Magn Reson B. 1996;110(3):316–320. [PubMed]

6. Vanhamme L, van den Boogaart A, Van Huffel S. Improved method for accurate and efficient quantification of MRS data with use of prior knowledge. J Magn Reson. 1997;129(1):35–43. [PubMed]

7. van den Boogart A. Quantitative data analysis of in vivo MRS data sets. Magn Reson Chem. 1997;35:S146–S152.

8. Laudadio T, Selen Y, Vanhamme L, Stoica P, Van Hecke P, Van Huffel S. Subspace-based MRS data quantitation of multiplets using prior knowledge. J Magn Reson. 2004;168(1):53–65. [PubMed]

9. Sandgren N, Selén Y, Stoica P, Li J. Parametric methods for frequency-selective MR spectroscopy—a review. J Magn Reson. 2004;168(2):259–272. [PubMed]

10. Knijn A, de Beer R, van Ormondt D. Frequency-selective quantification in the time domain. J Magn Reson. 1992;97:444–450.

11. Vanhamme L, Sundin T, Van Hecke P, Van Huffel S, Pintelon R. Frequency-selective quantification of biomedical magnetic resonance spectroscopy data. J Magn Reson. 2000;143(1):1–16. [PubMed]

12. Romano R, Motta A, Camassa S, Pagano C, Santini MT, Indovina PL. A new time-domain frequency-selective quantification algorithm. J Magn Reson. 2002;155(2):226–235. [PubMed]

13. Mierisova S, Ala-Korpela M. MR spectroscopy quantitation: a review of frequency domain methods. NMR Biomed. 2001;14(4):247–259. [PubMed]

14. Martin YL. A global approach to accurate and automatic quantitative analysis of NMR spectra by complex least-squares curve fitting. J Magn Reson A. 1994;111:1–10.

15. Allman T, Holland GA, Lenkinski RE, Charles HC. A simple method for processing NMR spectra in which acquisition is delayed: applications to in vivo localized 31P NMR spectra acquired using the DRESS technique. Magn Reson Med. 1988;7(1):88–94. [PubMed]

16. Golotvin S, Williams A. Improved baseline recognition and modeling of FT NMR spectra. J Magn Reson. 2000;146(1):122–125. [PubMed]

17. Soher BJ, Maudsley AA. Evaluation of variable lineshape models and prior information in automated 1H spectroscopic imaging analysis. Magn Reson Med. 2004;52(6):1246–1254. [PubMed]

18. Stoch G, Olejniczak Z. Missing first points and phase artifact mutually entangled in FT NMR data-noniterative solution. J Magn Reson. 2005;173(1):140–152. [PubMed]

19. Bottomley PA, Atalar E, Weiss RG. Human cardiac high-energy phosphate metabolite concentrations by 1D-resolved NMR spectroscopy. Magn Reson Med. 1996;35(5):664–670. [PubMed]

20. Bottomley PA, Weiss RG, Hardy CJ, Baumgartner WA. Myocardial high-energy phosphate metabolism and allograft rejection in patients with heart transplants. Radiology. 1991;181:67–75. [PubMed]

21. Provencher SW. Estimation of metabolite concentrations from localized in vivo proton NMR spectra. Magn Reson Med. 1993;30(6):672–679. [PubMed]

22. Provencher SW. Automatic quantitation of localized in vivo 1H spectra with LCModel. NMR Biomed. 2001;14(4):260–264. [PubMed]

23. Serrai H, Senhadji L, de Certaines JD, Coatrieux JL. Time-domain quantification of amplitude, chemical shift, apparent relaxation time T2, and phase by wavelet-transform analysis. application to biomedical magnetic resonance spectroscopy. J Magn Reson. 1997;124(1):20–34. [PubMed]

24. Kass M, Witkin A, Terzopoulos D. Snakes: active contour models. Intl J Comput Vision. 1988;2(3):321–331.

25. Soher BJ, van Zijl PCM, Duyn JH, Barker PB. Quantitative proton spectroscopic imaging of the human brain. Magn Reson Med. 1996;35(3):356–363. [PubMed]

26. Bonekamp D, Barker PB. MR spectroscopic imaging software for the neuroradiologist. ASNR, 43rd Annual Meeting; Toronto, CA. 2005.

27. in’t Zandt H, van Der Graaf M, Heerschap A. Common processing of in vivo MR spectra. NMR Biomed. 2001;14(4):224–232. [PubMed]

28. Cavassila S, Deval S, Huegen C, van Ormondt D, Graveron-Demilly D. Cramer-Rao bound expressions for parametric estimation of overlapping peaks: influence of prior knowledge. J Magn Reson. 2000;143(2):311–320. [PubMed]

29. Weiss RG, Gerstenblith G, Bottomley PA. ATP flux through creatine kinase in the normal, stressed, and failing human heart. Proc Natl Acad Sci USA. 2005;102(3):808–813. [PubMed]

30. Smith C, Bottomley PA, Gerstenblith G, Schulman SP, Weiss RG. Reduced ATP synthesis through creatine kinase in failing and non-failing hypertrophied human myocardium. Circulation. 2004;110(III610)

31. Cootes TF, Taylor CJ. Active shape models—‘smart snakes’. Proc British Machine Vision Conference; Springer-Verlag. 1992. pp. 266–275.

32. Lakshmanan S, Grimmer D. A deformable template approach to detecting straight edges in radar images. IEEE Trans Pattern Anal Mach Intell. 1996;18(4):438–443.

33. Staib LH, Duncan JS. Boundary finding with parametrically deformable models. IEEE Trans Pattern Anal Mach Intell. 1992;14(11):1061–1075.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |