1 Introduction
In research studies or in clinical investigations, images are used not only for evidencing the inner structures of the body but also for quantifying the available NMR parameters, like T1, T2 or the apparent diffusion coefficient, D. Especially, T2 evolution on the treatment time-course is more and more considered as a possible early apoptosis marker, hence an important NMR parameter to be taken into consideration in research or in clinical treatment [1,2].
In these conditions, the accurate determination of the relaxation time constants characteristic for targeted tissues, determined from a specific region of interest (ROI) on the reconstructed image becomes of major importance. Moreover, the complexity of the molecular motions leading to the nuclear relaxation in biological systems confers multi-exponential characteristics to the investigated systems. The need to quantitatively determine the relaxation time constants from non-exponential decays leads to some mathematical manipulations in order to extract more accurately the transformed signal amplitudes. All these reasons imply a certain degree of accuracy when relaxation calculations are performed using MRI data, close to the requirements usually encountered in spectroscopic techniques.
But what has become obvious in spectroscopy and almost unnoticed, becomes cumbersome and difficult when data are retrieved from images amplitude. The post-acquisition data processing proposed here implies phase correction of the reconstructed images, smoothing the retrieved decay curves by linear prediction – forward and fitting by the singular value decomposition method.
1.1 Phasing the images
The great majority of reconstruction methods implemented on the imaging machines are using the absolute amplitude of the FFT coefficients of the recorded signals, in order to obtain the final image. This method is simple to implement and has the great quality of avoiding the necessary phase corrections (used in spectroscopy). The absolute value reconstruction method uses the pair of images representing the real and imaginary components of the FFT transform to calculate the magnitude, I, of each corresponding complex pair:
This positive bias increases as the SNR decreases or as the nuclear sensitivity of the resonant nuclei decreases, as it is the case for sodium MRI. The presence of the positive bias level is affecting the accurate determinations of the NMR relaxation parameters such as T1, T2 as well as any result of pixels algebra. Nevertheless, even in the case of high SNR proton images, multi-exponential relaxation may be covered up by the positive noise level. So far, a few attempts were done to overcome the noise problem of the reconstructed images in MRI.
The use of power magnitude values, I2 [3] is reducing the noise effect over the region of interest, allowing the substraction of noise contribution. However, this power routine is only valid for exponential relaxation decays, while most of the biological samples are heterogeneous and thus non-exponential.
Another way to avoid the undesired biased noise produced by magnitude calculation is to phase-correct the images. If real, phase-corrected images are to be produced, they have the advantage of preserving the same noise characteristics as the original acquired signals. Ideally, the phase correction is achieved by rotating all the FFT coefficients in the complex plane, such that all information are transferred to the real component while the imaginary one tends to zero (noise level). An attempt to obtain phase-corrected images was done by van der Weerd et al. [4]. Their algorithm is based on the phase calculation of every pixel for the first and second echo images and the resulting correction is then applied for all echoes in the sequence. The first two echoes are thus becoming magnitude images while the rest of echo images are real, phase corrected ones. Nevertheless, when extracting the decay curve from the whole echo images train, a small bias is introduced for the first two points, characterised by the highest SNR values in the sequence.
Ideally, one should phase-correct all the images, using all available, acquired information. We have developed a novel and accurate method for phase correcting all images, by phasing each pixel separately, using all the echo images available in the sequence [patent pending]. The phase correction procedure is briefly shown in Fig. 2. Following the FFT procedure, the echo decays represent a mixing of real and imaginary intensities, as shown in Fig. 2a. The phasing procedure acts on each pixel of the images in order to maximise the real intensities while minimising the imaginary ones. Following the phase-correction procedure, the real phase-corrected images are obtained, allowing the retrieval of corrected amplitudes for the decay curves, Fig. 2b.
The phasing procedure applied on the spin–spin relaxation decay is shown in Fig. 2c, as compared to the absolute value decay. As expected, the phase corrected decay tends to the Gaussian noise level, while the absolute value decay tends to the constant value given by the averaged absolute value of the noise. The decays correspond to a solute NaCl sample characterised by an exponential spin–spin relaxation process. The SVD fit of an exponential decay is characterised by one singular value detaching from noise, as obtained for the phase corrected decay, Fig. 3. The presence of the constant positive value for the absolute decay induces a second singular value detaching from noise, indicating completely erroneous biexponential decay for the same sample.
The procedure allows both the accurate T2 determination for any ROI defined on an image obtained by a multiple-echo experiment, like multiple-slice multiple echoes (MSME), and correct image algebra, if required. We have applied the phase correction routine on sodium images, characterised by rather poor SNR values, in agarose-gel samples as well as in vivo mouse liver. The exponential decays thus obtained were fitted using the singular value decomposition method in order to obtain objective fitting parameters.
1.2 Data analysis
The relaxation decays retrieved from the MSME images are expected to represent a sum of, at least, two exponentials:
The data points, {si} obtained from the T2 decay are rearranged in a matrix form having a Hankel structure:
The requirements for using the SVD analysis of a decay curve are a constant sampling of the decay and the total description of the curve up to noise level. The noise plays an important role in this analysis since it represents the limit for the rank detection of the Hankel matrix. This objective fitting method depends on the Gaussian, unbiased noise level, since every data point is carrying the additive information regarding both the pure signal and the noise existing along the detection channel. In order to make clear the role of the noise in the SVD algorithm, we simulated several biexponential decays, similar to those determined for sodium in living tissue (liver) [7]. Hence we generated a sum of two exponentials (with time decays and amplitudes close to experimental values) for which the noise level was increased up to a level where the fitting started to be erroneous and the singular values were no longer detached from noise. Fig. 4 shows the detection of the two singular values characterising the biexponentiality of the simulated decay as a function of the noise level. The SNR factor is defined as [6]:
SVD fit results on simulated relaxation decays for different SNR values, without LP smoothing
SNR | T2long (ms) | T2short (ms) | Along | Ashort | Chi-square |
14 | 28.7 | 16.5 | 1896 | 200 | 0.885 |
20 | 28.9 | 16.2 | 1884 | 226 | 0.87 |
40 | 34.8 | 10.6 | 1491 | 765 | 0.792 |
80 | 35.1 | 10.7 | 1476 | 767 | 0.78 |
Simulated | 35 | 10 | 1500 | 750 | – |
Chi-square values are defined as [8]:
In order to increase the fitting accuracy we had to smooth the data. Again, unlike the spectroscopical CPMG sequence, the MRI inter-pulse distance is limited as well as the number of echoes due to the power deposition condition. For small animal imaging, the number of echoes may vary between 16 and 32, for clinical investigation usually being less. It is obvious that usual smoothing routines (interpolation, average) that are successfully applied in spectroscopy for 200–300 echoes cannot be applied here, since we cannot afford losing any data point information. A simple and quite satisfactory smoothing routine is based on linear prediction theory.
1.3 Data smoothing using linear prediction —forward
Usually, the linear prediction procedures are used in spectroscopy as a replacement of FFT procedure, in order to fully determine the characteristics of a dumped sum of sinusoids (frequency, phase, amplitude and width at half height) [9,10].
The mathematical model for NMR signals (i.e. the periodic, time-dependent free induction decay or exponential decay) is nonlinear. The nonlinearity can be circumvented by invoking the principle of linear prediction (LP) modelling that assumes that each data point can be expressed as a linear combination of a set of M consecutive data points. Using M previous data points to calculate subsequent ones is called forward linear prediction since the prediction is forward in time:
Usually, linear prediction is used to extend in time insufficient collected signals in both directions [7]. This is the case for missing points at the beginning of the acquisition array due to the spectrometer dead-time or rather long echo-time or at the end of the acquisition array for missing data towards null. The backward LP completes missing points at the beginning of the acquisition array, while forward LP completes missing points at the end.
A “side effect” of the LP is noise cleaning, since the experimental data are digitally reconstructed using the first data points having the highest signal-to-noise ratio.
The following example illustrates the use of forward linear prediction routine. First, six linear predictors f1 … f6 are calculated using a linear system of, for example, 10 equations as below.
The effect of forward linear prediction procedure on the SV levels is quite impressive. Fig. 5 shows again the dependence of SV (for the same generated curves as in Fig. 4), after performing the LP procedure.
Firstly, the singular values characterising the noise level are lowered by the linear prediction procedure, for all SNR cases. Consequently, the second SV rises above the noise level for the SNR = 40 case, becoming noticeable. Secondly, the procedure greatly improves the fitting results, especially for lower SNR cases (Table 2).
SVD fit results on simulated relaxation decays for different SNR values, following linear prediction – forward procedure
SNR | T2long (ms) | T2short (ms) | Along | Ashort | Chi-square |
14 | 31.7 | 6.72 | 1752 | 691 | 0.822 |
20 | 33.1 | 9.0 | 1619 | 705 | 0.81 |
40 | 34.3 | 10.3 | 1525 | 728 | 0.697 |
80 | 34.8 | 10.2 | 1499 | 749 | 0.69 |
Simulated | 35.00 | 10.00 | 1500 | 750 | – |
A crucial point on using a LP procedure is the choice of the number of the predictors and the number of equations used for calculation [10]. On the one hand, these numbers will impose the quantity of noise carried through the signal. On the other hand, increasing the number of predictors (and equations) increases the accuracy of prediction too. There is an optimum, depending on each particular case that should be looked for.
When processing very low SNR systems, the number of predictors should be larger than the number of exponentials expected to be present in the decay. For our study, taking into consideration the lowest SNR, we used six predictors and 10 linear equations.
Table 2 presents the fitting results on the noisy simulated decays (as in Table 1), after linear prediction – forward smoothing.
It is interesting to notice that the most beneficial LP treatment is obtained for the lowest SNR curves, for which the fitting results as well as the chi-square values are improved.
2 In vivo T2 relaxation results on mouse liver
In vivo experiments were performed on 17 nc/z female mice, normal and tumoral (Institut Curie breeding facility), on a Bruker Biospec system operating at 4.7 T. The RF coil was a home-made sodium (quadrature), proton (linear) birdcage coil, operating at 53/200 MHz. Single slice, multi-echo (16–32 echoes) images were recorded for sodium studies (TR = 350 ms, TE = 6.04 ms, FOV = 6.8 cm, matrix 64 × 64, slice thickness 6 mm) using 160 averages. Firstly, the images were phase corrected in order to avoid the rician noise created by modulus reconstruction. The Na relaxation curves, measured in regions of interest (ROI), were forward linear projected in order to reduce the noise amplitude and to reach the base line, required for a proper multi-exponential SVD fitting. Error estimations were performed on 100 runs by Monte Carlo simulations.
Fig. 6 shows the sodium MRI results obtained on a tumoral mouse liver as compared to control. ROI1 is a chosen area in the liver tumor, while ROI2 is chosen in the normal tissue. The sodium contrast is principally due to the changes in local cellular organization. Due to the quadrupolar characteristics of sodium ions relaxation in biological systems, the electrical changes of the environment are immediately reflected on its relaxation parameters. Generally, tumors are characterised by loose cellular packing, especially at the centre, where necrotic regions are developing. The necrotic regions are not empty spaces, but are full of extracellular juice of Na ions at 140 mM concentration. These circumstances allowed us to distinguish different areas in the liver. A biexponential fit of the decay could be obtained in the tumoral region due to the specificity of sodium ions relaxation and to the locally increased SNR. In the normal region, the decay appears monoexponential due to the following reasons. Firstly the cellular packing is more compact, allowing less extracellular space and consequently the fast relaxing component T2fast is expected to be close to the detection limits imposed by the time delay between the pulses of the CPMG sequence. Secondly, the total sodium content of the tissue is expected to be decreased, resulting in a lower SNR, preventing the second SV to rise above the noise level.
The relaxation parameters characterising the tumoral and normal tissue regions obtained after phase correction and linear prediction procedures are presented in Table 3.
Experimental results on tumoral mouse liver (ROI1) obtained in vivo as compared with control (ROI2)
ROI | T2long (ms) | T2short (ms) | Along/Ashort | SNR |
1 | 39.8 ± 4.2 | 4.21 ± 2.2 | 1.78 | 40 |
2 | 24.78 ± 3.7 | – | – | 25 |
They are clearly distinguishable on the basis of the corresponding relaxation properties, reflecting mostly the cellular packing of the tissue.
In conclusion, the use of a spin–echoes sequence in conjunction with post-processing techniques described here allowed to evidence the multi-exponential behaviour in tumoral regions of the mouse liver where the SNR is sufficient (about 40 or greater). Sodium ions' spin–spin relaxation holds therefore a great potential to become an endogeneous marker for tumor diagnosis and possibly for therapy response due to its sensitivity to the environmental changes.