Edited31P brain spectra using maximum entropy data processing

MAGNETIC RESONANCE IN MEDICINE 4,

385-392 (1987)

Edited 31PBrain Spectra Using Maximum Entropy Data Processing M. L. WALLER* AND PAUL s. Kings College, Strand, London, WClA 2LS, United Kingdom, and ?Institute of Neurology, Queen Square, London, WClN 3BG, United Kingdom Received September 17, 1986; revised January 14, 1986 Strategies for estimating peak areasin "P brain spectra are discussed,and the widespread convolutiondifference method is analyzed. The Skilling maximum entropy method (MEM) algorithm is applied to in vivo spectra and provides an estimate of the spectrum that is operator-independent, although at the expense of some negative bias. It may be possible Inc. to overcome this bias. 0 1987 Academic b,

INTRODUCTION

31PNMR spectroscopy (NMRS) of the human brain in vivo was first reported in 1983 (I).It was shown to be capable of safely assessing the state of cerebral cell metabolism in neonatal infants and thus assisting in clinical management. Semiquaniitativeinformation on metabolite concentrations was obtained from the relative peak areas; however, there were two problems in doing this: (1) many of the peaks overlapped and (2) a large broad signal from membrane phospholipid was superimposed on the metabolite signals, which had the effect of placing the metabolite peaks on a sloping baseline. Many of the published spectra were filtered to suppress the broad signal and enhance the resolution; however, both of these processes distort the relative peak areas. We wished to develop improved methods of estimating peak areas in human spectra. Several publications (2-4) on the application of a maximum entropy method (MEM) algorithm to NMR spectra suggested that this would be worth investigating. In the published sp:ctra resolution was apparently improved and background noise was suppressed. We therefore undertook a study of the usefulness of the MEM reconstruction in the estimation of peak areas in biological 31Pbrain spectra, starting with rat spectra which are similar to those of humans. The unedited spectrum from a rat brain, collected in vivo at 1.89 T, is shown in Fig. 1. It can be seen that the broad lipid line is not symmetrical. This is attributed to relaxation b y chemical-shift anisotropy (5). It has been shown that narrowband irradiation, cenlered between metabolite peaks, saturates all of the broad part of the lipid resonance uhile having only a minimal effect on the other lines (6, 7). So, in this sense, the contributions from all the lipids present may be regarded as a single, in-

* Present address: Department of Instrumentation and Analytical Science, University of Manchester Institute of Scien(:e and Technology, P.O. Box 88, Manchester, M60 IQD, UK. 385

0740-3194/87 $3.00 Copyright 0 1987 by Academic b, Inc. AU rights of reproduction in any form rserved.

386

COMMUNICATIONS PCr

20

10

0

-10

-20

ppm

FIG. 1. Unedited "P rat brain spectrum obtained at 1.89 T. Estimation of the metabolite peak areas is hindered by the peak overlap and the underlying broad resonance attributable mostly to lipids. The seven lines of interest are due to the 0, a,and y peaks of adenosine triphosphate (ATP), phosphocreatine (PCr), inorganic phosphates (Pi), phosphomonoester (PM), and phosphodiester (PD).

dependent, unwanted spectral line. (The use of saturation in human studies may be ruled out by the radiofrequency power deposition that would be required.) Two strategiescan be adopted for editing out the lipid contribution to the spectrum. The first is to collect the free induction decay (FID) using a spin-echo experiment. The relaxation time of the lipid is much shorter than the lines of interest and hence the formation of an echo will discriminate against it (8).However, this will introduce T2weighting, which cannot be corrected without an accurate knowledge of the lineshapes involved. The second strategy is to accept the limitations of the data acquired with a single surface-coil,single-pulse experiment and to use appropriate signal processing techniques to obtain an edited spectrum. This approach is already in use (I, 6), but the means that are currently employed to effect it lead to some quantitative inaccuracies in the spectra. Although it has not yet given improved quantitation, our proposed alternative is shown to be favorable in several respects. TRADITIONAL SIGNAL PROCESSING IN NMRS

The signal processing methods that have been traditionally employed in NMRS involve multiplying the collected FID by some apodizing function (9)(e.g., a decaying exponential for noise suppression). Multiplicative weighting in the time domain is, by a property of the Fourier transform (IO),equivalent to convolution in the frequency domain. Many shapes of filter have been proposed that trade off the conflicting requirements of resolution enhancement and noise suppression. They are generally unimodal functions that are specified by a small number of parameters. These are adjusted so that the resolution enhancement (due to the rising part of the function) is subordinated to noise suppression (due to the falling part) as the signal-to-noise ratio (SNR) deteriorates along the FID. The usefulness of this general approach is determined by the validity of assuming a single lineshape for all the spectral lines of interest (usually Lorentzian). Although this may be reasonable for high-resolution NMRS, it is not reasonable for in vivo studies where a variety of lineshapes and linewidths are encountered.

387

COMMUNICATIONS

The following time-domain filter function, U(t),has been used for brain spectroscopy

(a

U(t)=

1

1 + A exp(-Bt) ’

where A and 13 are free parameters to be set by the user, and is shown in Fig. 2. How this reduces the contribution of the lipid to the spectrum can be understood by expressing U(t)iis U(t)= 1 - V(t). [21 V(t)is thus given by the unimodal function V(t)=

A exp(-Bt) 1 + A exp(-Bt) .

[31

Since the convolution operation is distributive, it can be seen from Eq. [2] that the spectrum of a FID that has been multiplied by U(t) is the difference between the unfiltered spet:trum and its convolution with the Fourier transform of V(t),shown in Fig. 3. Thus the effect of U ( t )is to attenuate lines of width > B. Suitable choices of A and B will preferentially reduce the contribution of the lipid resonance (B !v 1000 Hz N 30 ppm at 32.5 MHz) or enhance resolution (B N 50 Hz _N 1.5 ppm). However, distortions of relative area will be introduced due to the variety of lineshapes in the spectrum. For the same reason, attempts to separate overlapping peaks by the use of multiplicative, time-domain filtering will also make peak area estimates inaccurate. REGULAR[ZED SPECTRUM ESTIMATION BY THE MAXIMUM ENTROPY METHOD

It is difficult to predict theoretically the shape of the lipid line that arises in 31P NMRS of the brain in viva It can be seen from Fig. 1 that it is difficult to assess it

3.00

2.00

4.00

6.00

8.00 10.00 12.00 14.00 16.00 18.00 Timelms

1.00

RG.2. Timedomain 6lter U(t)used to suppress the phospholipid resonance in brain spectra. The constants A and B (Eq. [ 11) are 4 and 400 Hz, respectively.

388

COMMUNICATIONS

3

FIG.3. Discrete Fourier transform of V ( t )(Eq.[2]). Parameters are the same as those for Fig. 2.

empirically. Hence an attractive alternate means of spectrum estimation is to discard those FID samples that are significantly affected by the lipid line and to attempt to estimate the spectrum from the remaining data that pertain only to the metabolite lines. In practice this requires us to reject the samples collected during a period of approximately 0.63 ms after excitation. The difficulty with this approach is that the spectrum estimation problem is now ill-posed since the reduced data set is incomplete. An infinite number of different spectra are consistent with the limited data set and so in order to make a single reliable choice (i.e., to regularize the problem), extra constraints must be imposed. Two pieces of prior knowledge are available to assist in this respect: (i) The absorption mode spectrum is nonnegative. (ii) The absorption and dispersion mode spectra are a Hilbert transform pair (11). In estimating the spectrum it is reasonable to demand that we choose the leaststructured spectrum that is consistent with both the experimental data and our prior knowledge. It has been shown that meeting this demand is equivalent to selecting the normalized spectrum, F(v), that has the constrained maximum of configurational entropy, H , (12) where H=-

s

F(v)ln[F(u)]du.

[41

It is important to distinguish between this definition of entropy and the one first proposed by Burg (13)who defined entropy as S where S= -

s

ln[F(v)] dv.

[51

Finding the spectrum with the data-consistent maximum of S has been shown to be

COMMUNICATIONS

389

equivalent to fitting the spectrum with an inverse polynomial (14). In our application, where the lineshapes are not known, there are no grounds for assuming that this model is appropriate. Finding th': data-consistent spectrum that maximizes the entropy H, constrained by conditions (i) and (ii), is thus the appropriate regularized estimation technique for the problem ,at hand. Note that it does not impose any particular lineshapes (we do not know what they are). It has been shown that lineshape information can be readily incorporated into the MEM algorithm (2-419, but it is not necessary for it to function. The implementation of our estimation technique was based on the Skilling Cambridge maxinium entropy data processing package (16). SIMULATION STUDY

To assess the inherent capability of MEM to estimate quantitatively accurate spectra, a simple spectrum was simulated comprising two Lorentzian-shaped lines (Fig. 4). They were chosen to be well separated so that errors in numerical integration would not obscure the effectiveness of the algorithm. Accuracy was assessed over a range of SNR values. SNR was defined, in the time domain, as the ratio of the signal magnitude to the root mean square noise level, evaluated at 1he first sampling point. The number of samples and the sampling interval for the FIDs were chosen to ensure that accurate numerical integration of spectral peaks could be carried out. The range of integration for each line was 10 times its full width at halCmaximum (FWHM). Most (94%)of the total area lies between these limits. Without dLscarding any of the simulated FID values, our MEM estimation technique was used to obtain each spectrum and the peak areas were obtained by numerical integration. 'The results, normalized to the true area that lies in a range of 10 times the FWHM, are shown in Fig. 5.

FREQ. / PPM

FIG.4. Simulated spectrum used to test MEM algorithm for accuracy and precision of area estimation.

390

COMMUNICATIONS 1A

m

12

c

. .2

1.0

m

c

*

g 0.8

U L

E 0.6 ‘5 0

0 line at 1Oppm 0 l i e at -1Oppm 0.4

20

40

60 80 FID SNR at T= 0

100

120

FIG.5. Estimated peak areas of the spectrum shown in Fig. 4 using MEM over a range of SNR values. The areas are normalized to the true area that lies in a range of 10 times the FWHM.

The so-called conservative behavior of the MEM algorithm causes peak areas to be underestimated at low SNR. Its reluctance to confirm the presence of any feature that could be due to noise alone leads to the suppression of the ‘‘wings’’ of lines. The overestimation at high SNR is due to the algorithm failing to completely suppress the noise in the spectrum and then forcing unsuppressed noise to be positive. Hence the estimated peak areas are systematically biased, in this case. The accuracy of interactive estimation techniques cannot be usefully quantified because they are strongly affected by the skill of the operator. RESULTS

The two spectrum recovery methods described above (i.e., conventional use of V(t) in Eq. [ l ] and MEM) were evaluated using the FID which, by means of manual phasing and the discrete Fourier transform, yields the spectrum shown in Fig. 1. This FID was obtained from a rat brain, in vivo, using a surface-coil probe. Single-pulse excitation was used and sampling, at 8 kHz,was commenced 0.125 ms after excitation, to avoid frequency-dependent phase errors (I7). After manual phasing, the filter function V (t )together with decaying exponentials for noise suppression was applied. Both of these types of filter were applied several times during an interactive session, with the free parameters being varied to achieve the preferred compromise among baseline flatness (i.e., removal of the lipid line), noise suppression, and broadening of the peaks of interest. The result is shown in Fig. 6. The spectrum in Fig. 7 was obtained by MEM processing. After manual phasing of the FID,the only requirements were (a) to estimate the noise level in the time domain, which is simply achieved by analyzing the data collected after the signal has dropped well below the noise level, and (b) to specify the number of leading FID values to be discarded. No further interaction was required to complete the processing.

39 1

COMMUNICATIONS

10

0

-10

-20

ppm

FIG. 6. Edited rat brain spectrum obtained interactively using conventional multiplicativetimedomain filtering (see text). DISCUSSION

The well-regularized incorporation of prior knowledge, (i) and (ii) above, using MEM, perm1ts a high degree of noise suppression, without the attendant loss of resolution that 1 s unavoidable with traditional methods of noise reduction. The MEM spectrum can be numerically integrated without further processing because the baseline is automatically placed at zero amplitude. This feature only fails in the (practilxlly unlikely) event of excellent SNR (>50). Spectra obtained by discrete Fourier transform will always have a nonzero baseline due to the effective discontinuity between the last and first points of any FID that is produced without spin-echo formation. Given speo%rathat are composed of lines of a known shape (e.g., Lorentzian), several groups (3, 1.5) have demonstrated impressive deconvolved spectra using algorithms that produce the data-consistent maximum of configurational entropy H as defined in [4]. Successful estimation of NMR spectra has also been achieved by maximizing Burg entropy S as defined in [ 5 ] under circumstances where the lineshape could reasonably be expected to conform to an autoregressive model (14). Our studies of simulated data have shown that our configuration of the MEM algorithm, which does not assume the presence of specific lineshapes, is not quanti-

10

0

-10

-20

ppm

FIG.7. Edited rat brain spectrum obtained automatically using MEM without any assumed lineshapes.

392

COMMUNICATIONS

tatively accurate. Note that the conventional interactive method (Fig. 6) also produces systematic bias through the use of the filter of Eq. 111. In a clinical environment, however, routine studies will produce many ostensibly identical spectra and it may only be their departures from some well-known “template” spectrum that are of interest. (Routine NMRS analysis in vitro often poses a similar problem.) It may be possible to use a template spectrum as a nonzero default in the MEM algorithm and to maximize the entropy with respect to this solution. This would provide estimates biased toward the normal peak area, rather than toward zero, and therefore might be clinically useful. It might be possible to calibrate, and hence remove, the bias altogether as follows: By simulating controlled variations of the areas of one or more peaks, in turn, and observing the corresponding changes in the estimated areas of the peaks of interest, the bias of the MEM algorithm could be measured. It would be necessary to determine the minimum FWHM of the broad lipid resonance that could be expected to occur in brain. This would indicate how much time elapses after excitation before the FID contains only signal from the metabolites. We might expect to produce an automated, operator-independent, quantitatively accurate technique for estimating edited spectra. Automated techniques for phasing a FID have been reported (e.g. (2)). ACKNOWLEDGMENTS We are indebted to Dr. John Skilling for the provision of the Cambridge maximum entropy method software and we would like to thank him and Dr. Ernest h u e for valuable discussions. The spectrum was collected on the spectrometer at University College London, with the assistance of Dr. Susan Wray, while one of us (P.T.) was in employment there. M.W. acknowledges support from the SERC. REFERENCES 1. E. B. CADY,A. M. DE L. COSTELLO,M. J. DAWSON,D. T. DELPY,P. L. HOPE,E. 0. R. REYNOLDS, P. S. TOFTS, AND D. R. WILKIE, Lancet 1 1059 (1983).

2. S. SIBISI, Nature (London) 301, 134 (1983). 3. E. D. LAUE,J. SKILLING,J. STAUNTON,S. SIBISI,AND R. G. BRERETON, J. Mugn. Reson. 62, 437 (1985). 4. E. D. LAUE,J. SKILLING,AND J. STAUNTON,“Fourth Annual Meeting of the Society for Magnetic Resonance in Medicine,” 191, 1985. 5 . A. C. MCLAUGHL~N, P. R. CULLIS,M. A. HEMMINGA, D. I. HOULT,G. K. RADDA,G. A. RITCHIE, AND R. E. RICHARDS, FEBS Left. 57,213 (1975). P. J. SEELEY, 6. P. s. TOFTS AND s. J. WRAY, J. PhpYiOl. 359, 4 17 (1985). 7. J. J. H. ACKERMAN, J. L. EVELHOCH, B. A. BERKOWITZ,G. M. KICHURA,R. K. DEUEL,AND K. S. LOW, J. Mugn.Reson. 56,3 18 ( 1984). 8. K. M. BRINDLE,M. B. SMITH,B. RAJAGOPALAN, AND G. K. RADDA,J. Mugn. Reson. 61,559 (1985). 9. J. C. LINDONAND A. G. FERRIDGE, Prog. NMR Spectrosc. 14,27 (1980). 10. R. N. BRACEWELL, “The Fourier Transform and Its Applications,” p. 208, McGraw-Hill, New York, 1978. 11. R. R. ERNST,J. Mugn. Reson. 1,7 ( 1969). 12. E. T. JAYNES, Proc. IEEE 70,939 (1982). 13. J. P. BURG,“Presented at 37th Meeting of Society of Exploration Geophysicists, Oklahoma City, 1967.” 14. V. V I ~ IP., BARONE,L. GUIDONI,AND E. MASSARO,J. Mugn. Reson. 67,91 (1986). 15. F. NI, G. C. LEVY,AND H. A. %HERAGA,J. Mugn. Reson. 66,385 (1986). 16. J. SKILLING AND R. K. BRYAN,Mon. Not. R. Astr. SOC.211, 111 (1984). 17. D. I. HOULT,C-N CHEN,H. EDEN,AND M. EDEN,J. Mugn. Reson. 51, 110 (1983).