Abstract
Two-dimensional (2D) NMR relaxometry has been widely used as a powerful new tool for identifying and characterizing molecular dynamics. Various inversion algorithms have been introduced to obtain the versatile relaxation information conveyed by spectra. The inversion procedure is especially challenging because the relevant data are huge in 2D cases and the inversion problem is ill-posed. Here, we propose a new method to process the 2D NMR relaxometry data. Our approach varies from Tikhonov regularization, known previously in CONTIN and Maximum Entropy (MaxEnt) methods, which need additional efforts to compute an appropriate regularization factor. This variety improves Butler–Reeds–Dawson algorithm to optimize the Tikhonov regularization problem and the regularization factor is calculated alongside. The calculation is considerably faster than the mentioned algorithms. The proposed method is compared with some widely used methods on simulated datasets, regarding algorithm efficiency and noise vulnerability. Also, the result of the experimental data is presented to test the practical utility of the proposed algorithm. The results suggest that our approach is efficient and robust. It can meet different application requirements.
1. Introduction
Relaxometry refers to the measurement of relaxation time in NMR. Relaxation time is known to correlate with many macroscopic material characteristic parameters like viscosity, crystallinity, concentration, and lots of microscopic parameters like length of polymer chain or kind of radical group, radius of molecule cluster [1]. This makes relaxometry analyze a competitive tool in analytical chemistry.
Generally, the spin-spin relaxation (T2) spectra and the spin-lattice relaxation (T1) spectra should be combined to boost the analysis. The distribution information of T2 and T1 conveyed by spectra was used to determine the relative amounts of free Mn2+ ions and chelated manganese ions when both species are present in the same aqueous solution [2]. Water-soluble manganese porphyrins can serve as MRI contrast enhancement agents, and manganese ions bound to enzymes or proteins are useful paramagnetic probes to elucidate the structures and functions of proteins and enzymes [3]. Given the information of T2 and T1, the values of 1/T2, 1/T1, and T1/T2 can be easily calculated. The information of T1/T2 and T2 was very useful to obtain information about the reaction of nanoparticles with molecular targets [4], based on the fact that biotinylated nanoparticles probes reacted with an avidin molecular target to form stable clusters permitting T2 and T1 to be measured as a function of cluster size. In these previous studies, T2 spectra and T1 spectra were measured by independent experiments.
In comparison with classical one-dimensional (1D) techniques, such as applying Carr-Purcell-Meiboom-Gill (CPMG) sequence for measuring T2, inverse-recovery sequence for T1 and pulsed gradient spin echo for diffusion coefficients (D), 2D NMR relaxometry technique measures two parameters at a single experiment. When using 1D spectra, there is a crucial risk that different structures in the real world would be grouped into the same one if their peaks overlapped with each other. In addition, when most of the peaks overlap heavily, 1D spectra will be much too complex for interpretation. By the introduction of an appropriate additional spectral dimension, 2D relaxometry produces 2D spectra. These spectra are simplified and lots of extra information is obtained from a single experiment. The overlapped peaks existing in 1D NMR relaxometry can be distinguished easily in a suitable 2D spectrum. Recently, 2D NMR relaxometry has been widely used as a powerful tool for identifying and characterizing molecular dynamics, especially in materials consisting of many complex components [5–7].
2. Methods
Generally, the physical and mathematic model of a 2D NMR experiment is described using the following equation [8]:where and are kernels modeling a relevant parameter for each, is a 2D spectrum, the symbol ‘T’ denotes matrix transpose, (,) is the sampling time, and stands for the experimental noise which is assumed to be white Gaussian noise with zero-mean. It is noteworthy that the kernels can vary in different forms depending on the type of experimental pulse sequences. In inversion-recovery-CPMG experiments, for example, the forms of the kernels are as follows [9]:
The kernels indicate the influence of a certain factor (T1, T2, or D) to the measured decay data. However, there is no closed form for . And data is collected at discrete values in practice. Consequently, a numeric approximation to is made. Equation (1) has a discretized form as follows:
The noise term is ignored or integrated into because noise should be incomparable to valuable data.
If the kernels are independent to each other, (3) can be easily restated into 1D form using Kronecker product:where , , . Here the operator represents vectorization that makes a vector by stacking columns of a matrix. The symbol ‘’ stands for Kronecker product of two matrices.
To determine from (4), an Inverse Laplace Transform (ILT) should be applied to experimental data. Up till now, the 2D ILT problem has been transformed into a 1D one. A careful reader may have a suspicion whether existing 1D ILT algorithms can be applied to solve the problem directly. Unfortunately, the data size of the transformed 1D problem is too large for 1D inversion algorithms.
In 2D case, is produced by two kernels using Kronecker product or a lexicographically order. Taking an inversion-recovery-CPMG experiment for example, just 8 CPMG echo trains with different polarization time are collected, and only 2048 echoes are stored at each CPMG train. Suppose that the preset T1-T2 grid is 128128 sized, then reading all the entries of will take 1GiB (2048812812832/8 Byte, in a 32-bit system) in memory. The procedure to calculate inverse matrix of is both memory and time-consuming. And any derivation involved optimization method would produce some other large temporary variables causing long computations. These problems may not be handled properly with an ordinary computer.
Additionally, the 2D ILT problem is notoriously ill-posed and ill-conditioned. We can get a solution even with a zero-fitting-residual, but in most cases, this accurate and precise solution will not be the exact one which contains relevant information about the true world. For example, in an ideal case, we set and the noise-free measured data , then we will get the “real” spectrum . However, in practice experimental data is inevitably collected with noise. If the measured data changes to due to noise, we will get immediately. The fitting residual is zero but the resultant spectrum does not have any practical significance. Small perturbations in measured data change the results by a large margin. The standard solution of the 2D inversion problem always produces a spectrum severely contaminated by noise.
To reduce the data size, a singular value decomposition (SVD) based compressing method [10] is widely used. Data compressing is performed with each kernel individually before combining two kernels into one matrix. The fundamental idea of this method is to truncate all the small singular values of to considerably reduce the data size. A trade-off between size and accuracy is made by the retaining number of singular values. The key point and the difficulty of this method are how to set the retaining number adaptively.
To circumvent the ill-posed problem, regularization techniques should be applied. In general, these methods can be divided into iterative regularization and direct regularization. The iterative regularization can be implemented fast and cheaply. Early iterations will produce good approximations and approximations degrade at later iterations. The extent of regularization corresponds to stop-ping criteria.
The method based on truncated SVD (TSVD) is one of the most popular iterative methods. TSVD should be performed iteratively to produce a nonnegative solution. An earlier version of TSVD [11] is provided to handle data with a good signal-to-noise ratio (SNR). To overcome the discontinuity produced by the former method, the algorithm proposed by JIANG et al. [12] in 1D inversion can also be used in the 2D case. In another way, TAN et al. [13] combine LSQR and TSVD together to get a more continuous spectrum. Lin et al. [14] studied the relationship between SNR of the sampling data and the optimal value of retaining number when handling a certain logging data. All the above-mentioned variants of TSVD method tend to be time efficient but may cause some artificial peaks due to singular value truncation. This drawback is evident when processing low SNR data.
A recently proposed iterative regularization method is based on the Trust-Region Algorithm for the Inversion (TRAIn) [15]. In TRAIn, the nonnegativity of the resultant spectrum is imposed by setting and the minimizing procedure is implemented with respect to . Meanwhile, the extent of regularization is controlled by a termination threshold relevant to noise level. The termination threshold is equivalent to the product of the estimated standard deviation of noise and a factor a little bit larger than 1. This method is reliable for 1D inversion, especially suitable for samples with nonsymmetrical distributed peaks in spectrum. Unfortunately, in 2D inversion cases, the spectrum shape is totally unknown due to the procedure of vectorization in (4). Additionally, 2D inversion is a large-scale problem and the truncated conjugate gradient method used to solve the trust-region subproblem will cause extensive computations for a given termination threshold. What is worse, several termination thresholds should be tried to select an optimal solution.
Other regularization methods that are widely used include Constrained Regularization (CONTIN) [16], Maximum Entropy (MaxEnt) [17], and other variants of Tikhonov regularization. All these variants of Tikhonov regularization lead to minimizing the following least-squares problem:where is a penalty function imposed to overcome the ill-posedness, is the so-called regularization factor included to balance the fidelity-to-data and regularity. Berman et al. [18] proposes a method to solve this problem, and the determination of regularization factor is a crucial issue when adopting these approaches. If is too small, the penalty term contributes little, then the solution will still be contaminated by noise a lot. If is too large, the solution will be a poor approximation of the original problem.
To choose an appropriate , many methods are devised, among which L-curve method [19], generalized cross validation (GCV) [20], and discrepancy principle [21] are well-known. The L-curve method supposes that the curve is shaped as an “L” when plotting () in log-log scale. Here is the corresponding solution to a given . Obviously, the value of should be varied in a wide range to form a curve. The optima is defined as the one that produces the pointer locating at the corner of the L-shaped curve. This method is easy to understand and implement, but time-consuming. Meanwhile, in some cases, the curve may not be a well L-shaped one and the locating algorithm based on curvature will fail. Besides the L-curve method, GCV is also time-consuming because many iterations over on the solution space are required. Unlike these two methods, discrepancy principle is much simpler and the computational cost will not be extravagant. If the noise level is known or can be estimated, the discrepancy principle chooses the that makes the residual equal to the noise level.
Despite regularization factor, penalty function is also important to direct regularization methods. The penalty term can take on different forms as follows:
For CONTIN, is defined as (see (6a)), where L is a low-pass operator (a second derivative operator is used as default in CONTIN) to enforce spectra smoothness because is believed to be mostly continuous. Unfortunately, due to this unequal penalization introduced by the second derivative operator, CONTIN tends to over-smooth weak peaks and under-smooth strong ones [22, 23]. Consequently, the resolution of spectrum may degrade.
Many other algorithms chose L as the identity matrix, if the real spectrum is smooth with a smaller Frobenius norm. In 2D inversion problem, L should not be chosen as second or other high-order derivative operator, because the vectored is a large vector (typically has 6464 rows or more) and high-order derivative operator will slow down the algorithm sharply.
The penalty function of MaxEnt is the entropy of the solution [17] (see (6b)), which is imposed to make the distribution of spectrum obeying Bayesian statistics. Then by employing the entropy function, information recovery is promoted. However, since the entropy term is undefined for negative values, all the negative values should be changed to be nonnegative after each iteration. This will cause a poor computational efficiency.
If there is a priori information that the spectrum consists of some “discrete” distributions, can be set to be L1-norm (see (6c)). The Iterative Thresholding Algorithm for Multiexponential Decay (ITAMeD) [24] used the Fast Iterative Shrinkage Thresholding Algorithm (FISTA) [25] to employ L1-norm minimization. ITAMeD enforces sparsity of the spectrum, i.e., among many solutions in agreement with experimental signal , the one with smallest number of components will be chosen as the resultant spectrum. Obviously, spectra produced by ITAMeD may rise sharply from the baseline. The character of enforced sparsity has made it a powerful tool to recover narrow peaks. But for broad distribution profiles, ITAMeD will not be a suitable alternative. Zhou et al. [26] extended this method to 2D cases to produce sparse T1-T2 spectra.
Unlike the variant methods reviewed above, Berman et al. [18] impose two penalty terms to confine more characteristics of the spectrum. The PDCO (Primal-Dual interior method for Convex Objectives) solver is employed to solve the complicated minimization problem. The main drawback of this approach employing multiple penalty terms is the extent allocation of each penalty function. Recently, Campisi-Pinto et al. [27] featured a method based on numerical experiments and simulations that maximize PDCO reconstruction accuracy under signal-to-noise conditions. They successfully put this method in 2D case.
Based on these arguments, we propose the application of Butler–Reeds–Dawson (BRD) algorithm to process the 2D NMR relaxometry data. BRD algorithm is designed to estimate solutions of first kind integral equations with nonnegative constraints and optimal smoothing [28]. Advantages like regularity, high efficiency, and adaptability have made this method a competitive alternative to solve ill-posed problems.
A priori knowledge of noise variance is needed because the quality of solution is connected to noise level in accordance with the Morozov’s discrepancy principle. To estimate the variance of noise from polluted data, Venkataramanan et al. [10] assume that the decay time is longer than the echo spacing. Hence any three successive data points should be located at the same straight line. The deviate extent is used to estimate the noise variance. Obviously, their assumption is bad when sampling fast-decay signals. To avoid this disadvantage, a wavelet-based method is proposed. The proposed noise estimation method employs a wavelet-based denoising filter to smooth the original CPMG echo trains. The filter uses near symmetric wavelets and calculates a universal threshold, and then soft thresholds the first level coefficients. Generally speaking, the soft-threshold approach involving wavelet transform packs most of the original signal energy into a few significant coefficients. The reconstructed CPMG train will be smoother. After a simple subtraction between the smooth version and the original version of measured data, the errors causing discontinuity and oscillation are extracted. To quantify the noise, the histogram of the amplitude of these extracted errors is calculated. The noise in CPMG train is assumed to be addictive white Gaussian with zero-mean in most relevant references. Accordingly, the histogram is fitted to the expression in (7). Immediately, the noise variance is evaluated by the standard derivation of the noise from fitting parameter .
Figure 1 shows the result of an experimental CPMG echo train. During the noise estimating, what we care about is the histogram of the subtraction of the original signal and the filtered signal. The histogram in the top-right subplot in Figure 1 shows that the noise is Gaussian.
After noise estimation, an adaptive SVD based compressing method is implemented to reduce the problem scale. The detailed procedures of the proposed compressing method, which sets the retaining number equal to the rank of kernel matrix, is given in supporting information. Then two kernels are rearranged into one using Kronecker product or a lexicographically order.
Up till now, the 2D inversion problem has been rewritten into an identical 1D form (see (4)) and the variance of noise has been estimated. This a priori information will be used by BRD algorithm to perform the standard Tikhonov regularization (see (8)).
The above minimum problem should be solved under the nonnegativity constraint. If we solve the problem directly, all the negative components calculated by optimizing iterations should be rewritten to nonnegative values. This will degrade the monotonicity of optimization algorithm. To overcome this drawback, BRD algorithm transforms the constrained problem to an unconstrained convex optimization one (see (9)).where
Here denotes Heaviside function and makes a diagonal matrix using the specified values. The nonnegativity constraint is imposed implicitly by a Heaviside function. Obviously, is positive definite. Immediately, we will get the first and second derivation of with respect to . A more detailed review of the BRD algorithm can be found in supporting information. In absence of any constraint, the minimization problem (see (9)) can be solved by any Newton-like method if regularization factor is given.
Now the only remaining problem is to choose an appropriate regularization factor . To connect the quality of spectra to the noise level of measured data, the norm of the fitting residual is assumed to be equal to the noise variance (see (10)). This assumption is easy to understand. Because there will be no fitting error if using noise-free data, then the fitting residual comes from noise.
In (12), is the solution of problem in (9) with a given ( is initialized to (15)); stands for the number of measured data. If was initialized with a certain value, we can easily get an optimal with respect to the given . Then the ideal value of should be calculated using the optimal . After updating , we will get the next optimal and then the next . After enough iterations, will converge to a fixed value and the corresponding will be output as the optimal solution. This update searching strategy used by BRD algorithm is advantageous over L-curve method and GCV, because it avoids lots of time-consuming unnecessary inversions.
Compared to other algorithms mentioned above, the BRD algorithm increases the computation efficiency greatly. However, there is a distinct drawback of this method in that the optimal solution is sensitive to the accuracy of the estimated noise variance. The update searching loop will not converge if the SNR of the measured data is too high or the estimated noise level is smaller than the actual one. To overcome this defect, some additional terminal criteria are imposed to improve the robustness of the BRD algorithm.
To demonstrate the improved BRD algorithm for the Inversion (BRDAIn) of two-dimensional NMR relaxometry data, we define the fit error as
The slope of the log-log graph of is proved [29] to be in the interval of :
Therefore, the larger the regularization factor is, the larger the fitting residual becomes. Additionally, based on numerical simulations, we find that the optimizing algorithm to solve the minimization problem (see (9)) converges much faster for larger .
Based on these assertions, is initialized to be a large value (at least larger than the optimal ) to further improve algorithm efficiency in the proposed method. The corresponding will be calculated with little computational effort. Once is calculated, the next is readily computed according to (12) and (13). As the update searching goes, decreases. If arrives at the estimated noise floor, the algorithm will stop. If changes too slowly, i.e., the slope of the log-log plot of is smaller than a preset tolerance TOL, the update searching should terminate because the extra time-consuming iterations make little sense to the output result. If becomes too small, the extent of regularization will be too low, consequently the algorithm terminates if is less than a preset minimum value . All the above-mentioned terminal criteria are shown in Figure 2.
It is noteworthy that the above-mentioned update searching method is first introduced by the original extension of the BRD method, but the algorithm efficiency of this original version can vary in a wide range. If the start point of is smaller than the optima one , the first several steps will cause long time computing. What is worse, it will take many more iterations to converge because may oscillate within a small interval containing . In our improved method, is initialized towhich is much larger than .
Considering our terminal criteria and initializing strategy together, when the update searching starts, will not increase. Consequently, no oscillating happens and the algorithm will converge rapidly.
In the following sections, the noise vulnerability, effectiveness, and spectra resolution of the proposed algorithm are analyzed. The comparison of our algorithm to other methods is also presented.
3. Materials and Experiments
The simulated and experimental data were processed using MATLAB on Windows 8 (64-bit edition), configured with Intel Core i5 3470 CPU running at 3.2GHz. The simulations were designed to illustrate the efficiency and robustness of the proposed algorithm.
3.1. Simulations
Simulation 1 (SIM-1): Algorithm Efficiency. For SIM-1, the forward-modeling spectrum used to mimic the inverse-recovery-CPMG experiment is unimodal. The spectrum contains a Gaussian peak centered at (T2, T1) = (10,100) microseconds. Some white Gaussian noise with zeros-mean were added to the simulated data to make the SNR of data equal to 100. Here SNR is defined as the ratio between the maximum value of signal and the standard derivation of noise.
Several methods including TSVD, CONTIN, TRAIn, ITAMeD, and the proposed BRDAIn method were applied to process the synthetic data. To compare the efficiency of each algorithm, time elapsed to invert the same data was recorded. For TSVD, the MATLAB implementation was shipped from Hansen’s Regularization toolbox [30]. CONTIN calculations were implemented using a MATLAB version [31] (, calculated by L-curve method). The code of TRAIn and ITAMeD we used can be found in the supporting information of the original articles. PDCO was implemented from PDCO method [18, 27]. For all these methods, the preset T1-T2 mesh grid is 64 rows and 64 columns.
Simulation 2 (SIM-2): Noise Vulnerability. SIM-2 used the same forward-modeling spectrum to simulate the inverse-recovery-CPMG experiment. Four sets of noisy data with different SNRs were generated by randomly adding white Gaussian noise to the simulated data. The SNRs of these synthetic data vary in a broad range, including 50, 25, 10, and 5. These datasets were processed by TRAIn, ITAMeD, and BRDAIn method.
3.2. Experiments
We have run inverse-recovery-CPMG experiments to test the practical utility of BRDAIn. The sample consists of two bottles of MnCl2 doped water with different concentrations (1.0g/L and 1.5g/L, respectively). Waiting times of inverse-recovery-CPMG sequence are 0.2ms, 0.57ms, 1.627ms, 4.642ms, 13.24ms, 37.765ms, 107.722ms, 307.267ms, 876.451ms, and 2500ms. Echo spacing and echo number of each CPMG are 0.235ms and 8000, respectively. The experiment was performed at 32 degrees Celsius on a 23.5MHz NMR spectrometer (MesoMR, produced by Suzhou Niumag Analytical Instrument Corporation, China). The entire dataset of this experiment can be found in supporting information. And the inversion spectrum using BRDAIn is shown in Figure 8.
4. Results and Discussion
Results of simulated and experimental data produced by afore-mentioned algorithms allowed us to evaluate the performance of BRDAIn in comparison to others.
For data with a high SNR in SIM-1, all methods yield a reasonable T1-T2 spectrum, but BRDAIn and PDCO finish the calculation in less than 1 minute. The BRDAIn has slightly less running time than PDCO.
It is noteworthy that none of the above-mentioned algorithms except BRDAIn and PDCO can handle the original T1-T2 inversion problem in SIM-1 in a relatively limited time. Therefore, the same data compressing procedure as BRDAIn was involved before CONTIN, TRAIn, ITAMeD, TSVD and PDCO calculations. Besides, all these five methods need prior calculations to get an appropriate regularization factor. In our experiment, the regularization factor was calculated using L-curve. The data in Table 1 just record the elapsing time with a certain value of the relevant factor. But for TSVD method, Hansen has implemented a fast L-curve method to calculate the truncating location. So TSVD can finish the whole calculation rapidly.
For CONTIN, the second derivative operator slows the calculations down sharply because the operator matrix is huge in 2D cases. TRAIn procedure is much faster in comparison with CONTIN. However, both CONTIN and TRAIn need enough times of calculations to form an L-curve to get a proper parameter. For ITAMeD, when the number of iterations increases, it will take much more time to converge.
It may be briefly summarized that only BRDAIn, PDCO, and TSVD can process the 2D relaxometry data in a short time. But as can be seen in Figure 3, there are lots of artificial small peaks in the TSVD result. Meanwhile, the whole CONTIN calculation takes about 26 hours (100 cases with different values of regularization factor were calculated to plot an L-curve). And the resultant spectrum is over-smoothed obviously. Therefore, both TSVD and CONTIN method were not adopted in the later experiments.
(a) TRAIn
(b) ITAMeD
(c) BRDAIn
(d) TSVD
(e) CONTIN
(f) PDCO
As the SNR of the simulated data decreases in SIM-2, the peak in the inverted spectra becomes broader for all inverting algorithms except PDCO. This is easy to understand because a higher SNR means a higher accuracy of sampling and an accurate peak bound and vice versa. In comparison with other algorithms except PDCO, this broadening effect is obviously suppressed when using BRDAIn. The PDCO method has a good ability to control the peak, but it may narrow the spectral peaks. Of course, it is also possible that I did not choose good parameters.
For results of ITAMeD (Figure 5) and BRDAIn (Figure 6), the noise vulnerability is about the same. But for TRAIn method (Figure 4), there is a notable feature that all the resultant spectra are symmetric. This is the combined effect of the setting in TRAIn method and the vectorization procedure to restate 2D inversion problem into 1D form. Consequently, if there is no prior knowledge about peak shape, TRAIn is not a good alternative. For the results of PDCO (Figure 7), as the SNR decreases, many pseudopeaks will appear in PDCO results. But if the SNR is high, the PDCO method is a good choice to obtain reliable results, especially when the SNR is over 1500 [27].
(a) SNR=50
(b) SNR = 25
(c) SNR = 10
(d) SNR = 5
(a) SNR=50
(b) SNR = 25
(c) SNR = 10
(d) SNR = 5
(a) SNR=50
(b) SNR = 25
(c) SNR = 10
(d) SNR = 5
(a) SNR=50
(b) SNR = 25
(c) SNR = 10
(d) SNR = 5
The above simulations have demonstrated the efficiency and robustness of BRDAIn. And the practical utility of our method will be shown from result of experimental data.
Mn2+ MRI contrast enhancement agents were commonly used in clinical. The concentration of water-soluble Mn2+ can affect the relaxation rate obviously. T1-T2 spectra can be used to rapidly determine whether the concentration was optimal. The experiment simply demonstrated the relation between concentration and relaxation time and it is only the beginning of what you can do using 2D relaxometry. Generally speaking, the more paramagnetic ions dissolved into water, the shorter the proton relaxation time is. Therefore, two bottles of MnCl2 doped water with different concentrations (the concentration gap should be large enough) will appear as two isolated peaks along the diagonal in T1-T2 spectrum.
5. Conclusion
Two-dimensional inversion algorithms with NMR had been intensively studied all over the world. However, existing methods are either too sensitive to noise or pay an exhausting importance to estimated noise variance. To overcome these shortcomings, a new method based on Butler–Reeds–Dawson algorithm is proposed. We demonstrate the performance of the proposed method by simulation and experiments. The results show that BRDAIn is both efficient and robust. Comparing to other methods mentioned in this paper, BRDAIn seems to be an effective extension of the existing algorithms for the inversion of 2D NMR relaxometry data.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by National Natural Science Foundation of China (no. 60972122), National Key Scientific Instrument and Equipment Development Projects (Grant 2013YQ170463), and Innovation Program of Shanghai Municipal Education Commission (no. 14ZZ135).