MRS (magnetic resonance spectroscopy), also known as NMR (nuclear magnetic resonance) spectroscopy, is widely used to identify relative abundance of isotopes of atoms, with unpaired nuclear spin, in molecules. The ubiquitous biological tissue imaging technique, MRI (magnetic resonance imaging), is based on MRS. The isotopes of interest in biochemistry, biology and organic chemistry include hydrogen-1 (proton), which is the most predominant, carbon-13, oxygen-17, sodium-23, and phosphorus-31, which are spin-aligned in their lowest stable quantum states in the presence of an applied magnetic field. If exposed to a sweep of radio frequency (RF) waves of the electromagnetic spectrum, these nuclei can absorb energy from the electromagnetic field and hop (i.e., flip the spin orientation) to the next higher energy quantum state. MRS is a record of relative numbers of nuclei, which hop to the higher quantum state, as the frequency is swept across a range.
The frequency at which a nucleus flips to the higher state varies according to the effective magnetic field experienced by the nucleus which in turn depends on the atom and its functional group (neighboring atoms). The dependence of the RF absorption frequency on the functional group allows H-atoms (and others listed above) in a molecule to be separated according to functional group. For example in benzyl alcohol, the H-atoms in the benzyl group, alkyl group and hydroxyl group all can be identified separately using MRS. That is because the effective magnetic field experienced by the H-atoms depends slightly on the countervailing magnetic field—counter to the applied field—due to the electrons of the functional groups. This counter magnetic field differs among functional groups.
The frequencies for absorption depend on the strength of the applied magnetic field. In order to have a measure of chemical shift independent of the applied field, the frequencies are converted into commonly used chemical shift in ppm (parts per million).
Typically, a 1-D MRS experiment sends a 90° electric pulse (duration is for a quarter wavelength or 90°) to a sample and measures the decay signal (FID or Free Induction Decay as a function in time) after the end of pulse. The decay signal is a record (in time) of the nuclei, which have hopped to higher energy state, returning (i.e., spin relaxing) to their ground state after the end of the pulse.
2-D MRS is obtained by sending a second pulse a short time after the first and measuring the decay signal. The two signals, one after each pulse, can be plotted on a 2-D graph (e.g., as a contour plot). This kind of spectroscopy using 2-D MRS experiments is referred to as 2-D COrrelation SpectroscopY or 2-D COSY. The first 2D COSY experiment was proposed by Jeener in 1971, and carried out by Ernst and colleagues in 1975 (see Jeener, “The Unpublished Basko Polje (1971) Lecture Notes About Two-Dimensional NMR Spectroscopy”, NMR and More, in the Honour of Anatole Abragam, Les Editions de Physique, Paris, 1994; Müller, Kumar, and Ernst, “Two-dimensional carbon-13 NMR spectroscopy”, J. Chem. Phys. 63, 5490-1, 1975; Aue, Bartholdi, and Ernst, “Two-dimensional spectroscopy: Application to nuclear magnetic resonance”, J. Chem. Phys. 64, 2229-46 (1976).
2-D COSY can be used to identify and study coupled resonances that cannot be observed and distinguished using conventional 1-D spectroscopy. It can therefore provide a much more informative snapshot of the biochemistry of the brain than is achievable in 1-D. 2-D COSY is starting to be used for studying neurological brain disorders (see, for example, Watanabe, et al, “A Post-processing framework for localized 2-D MR spectroscopy in vivo”, Magnetic Reson Med Sci., 12(3):215-21, 2013; and Dufour, et al., “2-D COSY NMR Spectroscopy as a Quantitative Tool in Biological Matrix: Application to Cyclodextrins”, AAPS Journal, 17(6): 1501-1510, 2015).
Some general purpose methods for processing raw 2-D COSY MRS data can be found in Zhou, “2-D NMR Spectrum Processing with MNOVA”, 2012; and also in “Acquisition and Processing of 2-D NMR Spectra using TOPSPIN 3.5”, 2016.
The invention concerns a method and a system for pre-processing of noisy 2-D COSY MRS data to produce a clean 2-D COSY spectrum.
Raw 2-D COSY data needs to be processed to eliminate noise and enhance the signal-to-noise ratio (SNR) for visualization and precise identification of protons and other nuclei of the functional groups of the molecules. The result of processing raw data, dubbed pre-processing, is to create a clean 2-D COSY signal, which is a surface defined on a 2-D frequency axes, converted to ppm.
Metabolite concentrations extracted from the clean 2-D COSY signal produced as output can, among other applications, be used for diagnosing and characterizing neurological disorders of the brain, including PTSD (Post Traumatic Stress Disorder) and mTBI (mild-Traumatic Brain Injury). Since acquisition of raw MRS data are performed in vivo, the diagnoses of neurological disorders by examining MRS data has been dubbed “virtual biopsy”.
In general, according one aspect, the invention features a system and a method for pre-processing raw proton 2-D COSY MRS data to produce a clean 2-D COSY signal. Pre-processing steps include: singular value decomposition (SVD)-based channel weighting, spectral registration, apodization, zero-filling in the time domain, conversion to the frequency domain, and final adjustment of the peak locations using known metabolites with high SNR as landmarks.
In general, according to one aspect, the invention features a magnetic resonance spectroscopy (MRS) pre-processing system. The system comprises one or more MRS machines that produce raw two dimensional (2-D) correlation spectroscopy (COSY) MRS data obtained from a pool of subjects using the same MRS machine and a computer system executing a pre-processing tool that pre-processes the raw 2-D COSY MRS data to produce clean 2-D COSY signals by performing channel weighting, spectral registration, apodization, zero-filling in the time domain, conversion to the frequency domain, and adjustment of the peak locations using known metabolites.
Preferably, the channel weighting of raw data is singular value decomposition (SVD)-based channel weighting and the spectral registration is performed on multiple iterations of raw data of a subject to correct for drift across the iterations. The registered iterations can be averaged. Further, the spectral registration can be performed in the time domain.
The apodization can include applying a Tukey window filter.
The peak locations can be adjusted using known metabolites by mapping non-uniformly spaced piecewise ppm axes, and corresponding spectral values in the F1-F2 plane to uniformly spaced ppm axes using interpolation. The known metabolites might include N-acetylaspartate (NAA), creatine (Cre), and choline (Cho).
Further adjustment of the peak locations using known metabolites can also comprise comparing metabolites for different subjects and then correcting for individual differences in the axes for each subject caused by slight variations in scanner frequency values, and for subject-specific errors in the locations of landmark singlets.
In general, according to another aspect, the invention features a magnetic resonance spectroscopy (MRS) pre-processing method. This method comprises performing channel weighting on raw 2-D COSY MRS data, spectrally registering the 2-D COSY MRS data, apodizing the 2-D COSY MRS data, zero-filling the 2-D COSY MRS data in the time domain, converting the 2-D COSY MRS data the frequency domain and adjusting of the peak locations using known metabolites to produce clean 2-D COSY MRS data.
The above and other features of the invention including various novel details of construction and combinations of parts, and other advantages, will now be more particularly described with reference to the accompanying drawings and pointed out in the claims. It will be understood that the particular method and device embodying the invention are shown by way of illustration and not as a limitation of the invention. The principles and features of this invention may be employed in various and numerous embodiments without departing from the scope of the invention.
In the accompanying drawings, reference characters refer to the same parts throughout the different views. The drawings are not necessarily to scale; emphasis has instead been placed upon illustrating the principles of the invention. Of the drawings:
The invention now will be described more fully hereinafter with reference to the accompanying drawings, in which illustrative embodiments of the invention are shown. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein; rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.
Unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the relevant art and will not be interpreted in an idealized or overly formal sense unless expressly so defined herein.
To avoid repetitions the acronym “COSY” will denote “2-D COSY” hereinafter whether the prefix “2-D” is present or not.
Without the adjective “raw”, COSY or 2-D COSY will denote a pre-processed signal. The adjective “clean” will normally denote the final pre-processed signal although it will also denote the MRS signal at intermediate stages of pre-processing. The main objective is to distinguish MRS signal at various stages of pre-processing from the raw MRS data.
Ideally, because of lack of standardization of MRS machines, the pre-processing system tool 300 should have parameters that are set specific to MRS machine 110.
In embodiments, the MRS data was acquired from the same brain regions of subjects, a specific voxel of size 20 mm×20 mm×20 mm. The MRS data in this document are labeled by subjects, identified as “A00ijk” with i, j and k being digits 0 through 9.
From hereinafter, “iteration” will be used for the 8 iterations and “average” will be used in the traditional sense of “sum of values divided by the number of values”.
The raw 2-D COSY MRS data acquisition process is similar to that of 1-D MRS. Multiple iterations of the same pulse sequence are executed, and each coil in the MRS scanner collects the HD signal generated by each iteration. The word “coil” refers to the metal coil carrying electric current that generates a magnetic field along the coil axis. In a 2-D scan a second frequency axis is generated by introducing an evolution delay (T1) into the pulse sequence, and repeating the scan multiple times, each time adding an additional T1 delay. At each T1 increment, multiple iterations are collected.
To avoid confusion, the word “channel” is used for “coil”.
For typical values of the raw data parameters, number of channels=32, number of iterations=8, number of T1 increments=64 and number of T2 increments=2048, raw 2-D COSY MRS data matrix (synonymous with FID matrix) is a 4-dimensional matrix of complex numbers with dimensions [num increments of T1 (64)×num samples of T2 (2048)×num iterations (8)×num channels (32)]. The two dimensions in “2-D” COSY refer to T1 and T2, which will eventually be transformed by FFT (Fast Fourier Transform) to frequency values and then to chemical shifts in ppm, F1 and F2, respectively.
To clarify T1 and T2 are in time units. F1, corresponding to T1, and F2, corresponding to T2, are in frequency units after FFT and then in ppm (defined later).
Clean 2-D MRS signals are produced after pre-processing raw MRS data:
Raw 2-D COSY MRS data is typically too noisy and too high dimensional, 4, for visualization and quantification of the metabolite peaks and valleys, and for performing statistical analysis.
The relevant metabolites for diagnosing neurological abnormalities among combat veterans (e.g., Post-Traumatic Stress Syndrome, PTSD and mild-Traumatic Brain Injury, mTBD are glutamate (Glu), glutamine (Gln), gamma-amino butyric acid (GABA), creatine (Cre), lactate, N-acetylaspartate (NAA), myo-inositol (mI), threonine, valine and choline (Cho).
For the dataset used, the raw data consists of 8×32=256 surfaces, defined on 64×2048 points on the 2-D plane, which are impossible to visualize and difficult to analyze. Instead, a 2-D COSY scan, after pre-processing, is typically analyzed in the frequency domain, since that is where the metabolite locations are defined. It is visualized as a 3-D surface or a 2-D contour plot, like the one shown in
However, to get from 4-D time-domain data to a 3-D frequency domain surface requires several pre-processing steps. The pipeline described here is designed to mitigate underlying signal quality issues commonly found in raw 2-D COSY data, increasing the SNR of the final 3-D metabolite surface. This enables more accurate identification of metabolite locations and quantification of their concentrations.
SVD-based channel weighted averaging 215:
Instead of weighting each channel equally, an SVD-based algorithm, described in Rodgers and Robson, “Receive array magnetic resonance spectroscopy: whitened singular value decomposition (WSVD) gives optimal Bayesian solution”, Magnetic Resonance in Medicine, vol. 63, pp. 881-891, 2010, is used to derive channel weightings from the raw scanner data. The data for the 32 channels are summed with the weightings produced by the algorithm, creating a [64×2048×8] matrix from a [32×64×2048×8] matrix. The channel dimension thus drops out, reducing the raw data dimension by 1, i.e., from 4 to 3.
Spectral Registration of Iterations 220:
Typically, multiple iterations (8 in the dataset used here) of each of the T1 scans are acquired, and then averaged using channel weighting 215 to increase the SNR of the final averaged signal in the time domain. Often, these signals drift across iterations, which will decrease the SNR of the final averaged signal. Possible drift is corrected by applying the spectral registration algorithm described in Near, et al., “Frequency and phase drift correction of magnetic resonance spectroscopy data by spectral registration in the time domain”, Magnetic Resonance in Medicine, vol. 73, pp. 44-50, 2015.
Averaging Iterations 225:
The next step in
Apodization 230:
To increase the SNR of the signal, a filter is applied to each column of the [64×2048] FID (time domain) matrix, a process which is referred to as apodization. Each column of data contains the 64 T1 data points for one of the 2048 T2 values. A Tukey window filter is used, although there are many other filters types that can be used for apodization.
Zero-Filling 235:
The apodized time-domain data is then zero-filled in T1 so that it more closely matches the spectral resolution of the signal in T2 once it is converted to the frequency domain using 2-D FFT. The T2 sampling period in the dataset used to develop this technology was 0.25 ms, giving it a bandwidth of 2000 Hz spread across 2048 samples and a frequency domain resolution of 0.977 Hz/sample. The T1 sampling period was 0.8 ms, giving it a bandwidth of 625 Hz spread across 64 increments and a frequency domain resolution of 9.77 Hz/sample. To more closely match the frequency domain resolution in F2, one zero-pads the columns of the T1-T2 matrix to a size of [640×2048]. When converted to the frequency domain using the 2-D FFT, the frequency domain resolution in F1 will be 0.977 Hz/sample.
Transforming to Frequency Domain (FFT) 240:
The signal is then converted to the frequency domain using the 2-D FFT.
Conversion of frequency axes in Hz to chemical shifts in ppm 245:
The ppm axis is derived from the frequency axis, f, by: ppm=4.7−(f/freq), where 4.7 is the standard offset value (in ppm) of the residual water peak located around 0 Hz, and freq is the scanner frequency in MHz. This parameter is pulled directly from the raw scanner data file, and is somewhere in the range of 123 (MHz), although it is not identical for all subjects. Note that f and freq are in same units.
Alignment of 2-D Spectra to Landmark Singlet Peaks 250:
Due to issues related to data acquisition and relaxation of nuclei, the locations of metabolites in the 2-D spectrum may not end up where they theoretically should be in the chemical shift (ppm) plane after step 245 of the processing pipeline. These location shifts are most noticeable visually, and easiest to detect analytically by focusing on metabolites with the largest amplitude peaks in the spectrum, such as the singlets of N-acetylaspartate (NAA), creatine (Cre), and choline (Cho). These relatively large amplitude peaks should be located close to the diagonal of the spectrum, where ppm values are equal. For instance, the singlet NAA peak should be around [2.008, 2.008] ppm, and Cre should be located close to [3.027, 3.027] ppm in vivo 2-D COSY scans of the brain. If these peaks have drifted from their theoretical chemical shift locations, it is likely that the other metabolites have also drifted, although it is much harder to detect in low amplitude peaks. Therefore, high SNR singlets are used as landmarks to shift the 2-D spectra so that these landmarks are located closer to where they should be on the diagonal of the chemical shift plane. The procedure is as follows.
First, the values that lay as close to the diagonal as possible in the 2-D COSY matrix are extracted. The frequency locations corresponding to the maximum amplitude of the signal closest to a set of landmark peaks are identified. In the dataset used to develop this technology, the following peaks are used as landmarks: NAA (2.008 ppm), total Cre (tCre) (3.027 ppm), Cho (3.185 ppm), mI (3.5217 ppm), and Cre (3.913 ppm), although other peaks can be used depending on the characteristics of the dataset.
In F1 and F2, a new ppm axis is constructed piecewise between each pair of landmarks by finding the equation of the line intersecting the theoretical peak locations (in ppm) and the empirical peak locations (in Hz). The non-uniformly spaced piecewise ppm axes, and corresponding spectral values in the F1-F2 plane are then mapped to uniformly spaced ppm axes using interpolation. If there are multiple subjects in the population whose metabolites are being compared, this same process should be applied to all the subjects identically, using the same final uniformly spaced ppm axes for interpolation. This process therefore corrects for individual differences in the ppm axes for each subject caused by slight variations in scanner frequency values, and for subject-specific errors in the locations of landmark singlets after step 245 of the processing pipeline.
Since the re-alignment process only utilizes the diagonal (ppm-values) of the COSY matrix, and not the actual 2-dimensional peak magnitudes, it remains to be checked whether the final locations of the peaks had been adjusted properly. After the alignment step, 250, the [F2, F1.] coordinates of the largest magnitude peak in the region around each reference singlet peak (in this case NAA and total Cre, tCre, peaks), original and adjusted, were retrieved from the 2-D COSY spectrum from the entire subject population and plotted in
All of the processing steps in this pipeline were implemented in Python using components of the OpenMRSLab project's Suspect toolbox (https://github.com/openmrslab/suspect).
Overall, a 2-D COSY pre-processing method produces clean 2-D COSY signals from raw data that can be directly compared across a population of subjects.
While this invention has been particularly shown and described with references to preferred embodiments thereof, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the scope of the invention encompassed by the appended claims.
This invention was made with Government support under contract number W81XWH-10-1-0785, awarded by the U.S. Department of Defense. The Government has certain rights in the invention.