Iterative Thresholding Algorithm for Multiexponential Decay Applied to

Jan 8, 2013 - Pulsed gradient spin echo (PGSE) is a well-known NMR technique for the determination of diffusion coefficients and the analysis of mixtu...
1 downloads 8 Views 2MB Size
Download PDF
Article pubs.acs.org/ac

Iterative Thresholding Algorithm for Multiexponential Decay Applied to PGSE NMR Data Mateusz Urbańczyk,† Diana Bernin,‡ Wiktor Koźmiński,† and Krzysztof Kazimierczuk*,§ †

Faculty of Chemistry, University of Warsaw, Pasteura 1, 02-093 Warsaw, Poland Swedish NMR Centre, University of Gothenburg, Box 465, 40530 Göteborg, Sweden § Centre of New Technologies, University of Warsaw, Zwirki i Wigury 93, 02-089 Warsaw, Poland ‡

S Supporting Information *

ABSTRACT: Pulsed gradient spin echo (PGSE) is a well-known NMR technique for determining diffusion coefficients. Various signal processing techniques have been introduced to solve the task, which is especially challenging when the decay is multiexponential with an unknown number of components. Here, we introduce a new method for the processing of such types of signals. Our approach modifies the Tikhonov’s regularization, known previously in CONTIN and Maximum Entropy (MaxEnt) methods, by using the S1-norm penalty function. The modification enforces sparsity of the result, which improves resolution, compared to both mentioned methods. We implemented the Iterative Thresholding Algorithm for Multiexponential Decay (ITAMeD), which employs the S1-norm minimization, using the Fast Iterative Shrinkage Thresholding Algorithm (FISTA). The proposed method is compared with the Levenberg−Marquardt-Fletcher fitting, Non-negative Least Squares (NNLS), CONTIN, and MaxEnt methods on simulated datasets, with regard to noise vulnerability and resolution. Also, the comparison with MaxEnt is presented for the experimental data of polyethylene glycol (PEG) polymer solutions and mixtures of these with various molecular weights (1080 g/mol, 11 840 g/mol, 124 700 g/mol). The results suggest that ITAMeD may be the method of choice for monodispersed samples with “discrete” distributions of diffusion coefficients.

P

In order to determine A(D) from experimental data, S(g), which is the inversion of Laplace Transform (ILT), must be applied. Unfortunately, direct numerical ILT is known for strong vulnerability to noise and numerical instability. To circumvent these problems, various approaches have been proposed.1,2,4 In general, they can differ in total band shape and singlefrequency methods. The most important examples of the total band shape approach are Direct Exponential Curve Resolution Algorithm (DECRA), 5 Speedy Component Resolution (SCORE),6 and Multivariate Curve Resolution (MCR).7 The simplest single-frequency method is plain curve fitting, such as, for example, the Levenberg−Marquardt algorithm8 in its default version. Others include CONTIN, 9 Spline Model (SPLMOD),10 and Maximum Entropy.11 Notably, an algorithm described in this paper can be also classified as a singlefrequency method; therefore, we will further consider only methods of this type for comparison.

ulsed gradient spin echo (PGSE) is a well-known NMR technique for the determination of diffusion coefficients and the analysis of mixtures.1,2 The technique is based on the measurement of a signal attenuation that is due to the diffusion of studied compounds during the time between coding and decoding magnetic field gradient pulses. A resulting signal attenuation was described by Stejskal and Tanner in 1965,3 using the following equation: 2 2 2

S(g ) = S(0)e−Dγ g

δ Δ′

(1)

where S(g) is a signal magnitude that is dependent on the magnetic field gradient strength g, S(0) the signal magnitude at zero magnetic field gradient strength, D the diffusion coefficient, γ the gyromagnetic ratio of the studied nucleus, δ the duration of magnetic field gradient pulse, and Δ′ the effective diffusion time. The attenuation can be also described as the Laplace transform of a distribution of diffusion coefficients (A(D)): Ψ=

S(g ) = S(0)

∫D

Dmax

2 2 2

A(D)e−Dγ g

δ Δ′

dD

min

Received: November 2, 2012 Accepted: January 8, 2013 Published: January 8, 2013

(2)

where Ψ is the signal attenuation. © 2013 American Chemical Society

1828

dx.doi.org/10.1021/ac3032004 | Anal. Chem. 2013, 85, 1828−1833

Analytical Chemistry

Article

The simplest example of the first group is the Levenberg− Marquardt algorithm,8 where the number of components is assumed directly. SPLMOD,10 which is a more sophisticated method, requires only the maximum number of components to be set and seeks the optimal solution automatically. The first widely used algorithm from the second group was CONTIN,9 which was shown to be optimal for signals with multiple (and unknown) numbers of components. Then, in 1998, Delsuc and Malliavin implemented the Maximum Entropy (MaxEnt) algorithm to process diffusion spectra.11 Delsuc et al. reported that the method performs better than CONTIN in resolving similar diffusion coefficients. Importantly, both CONTIN and MaxEnt can be thought of as Tikhonov regularization12 of the following least-squares minimization problem: ΦA − Ψ

2 S2

+ τ Θ(A)

(3)

where Φ is the Laplace Transform Matrix, Θ(A) is a regularization function, and τ keeps a balance between the accordance of A with the experimental data and the regularization term. The Θ(A) is defined as follows for the MaxEnt method:

Figure 1. Comparison of the convergence rate of ITAMeD and ISTA.

In practice, to process experimental data S(g) and obtain A(D), one also needs to make certain assumptions about D. In general, one of two approaches can be applied. First, one can assume that only given (usually small) number of values of D correspond to signal components. Such a procedure allows one to find both D and A(D), but requires proper prediction of several components. The second approach is the representation of A(D) and D as vectors A and D, respectively. The array contains assumed values of diffusion coefficients that are possibly present in the signal. This signal is then “decomposed” into a sum of exponentials with decay rates from the vector D.

Θ(A) = −∑ i

⎞ ⎛ Ai Ai ⎟ ⎜ log⎜ ⎟ ∑j Aj ⎝ ∑j Aj ⎠

(4)

whereas, for CONTIN, it is

Θ(A) = A

S2

(5)

Alternatively, instead of a regularization term, one may use a non-negativity constraint to the least-squares problem (non-

Figure 2. Comparison between the results of (A) ITAMeD, (B) CONTIN, (C) MaxEnt, and (D) NNLS processing of the simulated datasets. For each method, 100 diffusion coefficient distributions were calculated (marked with differently colored lines). The black dashed line represents the set diffusion ratio and magnitude for simulation. 1829

dx.doi.org/10.1021/ac3032004 | Anal. Chem. 2013, 85, 1828−1833

Analytical Chemistry

Article

Figure 3. Comparison between the results of (A) ITAMeD, (B) CONTIN, (C) MaxEnt, and (D) NNLS processing of the simulated datasets with two peaks at varying ratio of diffusion coefficients. The red dashed line represents the reference value set in the simulation.

Θ(A) = A

Sp

(6)

(usually 0 < p ≤ 2). Various settings of p were extensively discussed in papers on signal processing.15,16 In 2004, Daubechies et al.17 discussed the advantages of setting p = 1. In such a case, the optimization problem can be described as follows: min ΦA − Ψ A

2 S2

+τ A S1

(7)

Importantly, the minimum can be found by using very robust iterative thresholding algorithms (ISTAs).18 The approach based on the S1 norm promotes sparsity of the result, i.e., among many solutions A that are in agreement with experimental data Ψ, the one with a smallest number of components is chosen.17 In other words, there is an implicit assumption that the number of components is as small as possible (but is not set directly!), and that corresponding diffusion coefficients belong to the vector D. In this sense, the S1 norm regularization can be thought of as a hybrid of the two groups of algorithms mentioned above.

Figure 4. Result of ITAMeD processing for a simulated signal having a “broad” distribution of diffusion coefficients.

negative least-squares, NNLS13). This method, although not commonly used in PGSE NMR, is important to us, because it was reported14 to reveal similar sparsity-enforcing properties as the algorithm introduced in the present work (see below). Interestingly, there are many more possibilities for the choice of Θ(A) that have not been widely exploited in the processing of diffusion signals. For instance, the regularization used in CONTIN can be generalized to an arbitrary norm S p : 1830

dx.doi.org/10.1021/ac3032004 | Anal. Chem. 2013, 85, 1828−1833

Analytical Chemistry

Article

Simulations. We compared our approach with the Levenberg−Marquardt-Fletcher exponential fitting, NNLS, CONTIN, and MaxEnt, with regard to the noise vulnerability and achievable resolution. Tests for Noise Vulnerability. The signal used for simulations was chosen to mimic one of the experiments described below. It consisted of three-exponential decay with the following diffusion coefficients: 1.6 × 10−11 m2/s, 6.3 × 10−11 m2/s, and 2.3 × 10−10 m2/s. The component magnitudes were in the ratio of 3:1:2. The simulation was repeated 100 times, each time with random Gaussian noise at four different levels, from 0.1% to 1% of the first data point (generated for each run independently). Sixty four (64) exponentially distributed gradients were used (see the SI for details). The vector D contained 32 equally spaced elements (from 10−11 m2/s to 10−9.5 m2/s). NNLS calculations were performed using built-in MATLAB function lsqnonneg. CONTIN calculations were done using MATLAB emulation of the CONTIN program20 (alpha = 5 × 10−5), MaxEnt using the NMR Processing Kernel (NPK) (DOSY2D module with option me_preset_value = 5),21 and the Levenberg−MarquardtFletcher (LMF) fit using the LMFnlsq MATLAB script.22 The LMF fitting was done only for diffusion coefficients, not for component amplitudes, which were fixed to their proper values (3:1:2). The ITAMeD procedure was performed with 10 000 iterations for all simulations and experiments described in the text. Figure 1 presents the convergence for one of our datasets. Other datasets revealed very similar convergence rates, which ensures that 10 000 iterations is generally sufficient. For comparison, the convergence curve of regular ISTA is also depicted in Figure 1. Importantly, the S1 norm obtained after 10 000 iterations is still higher than for FISTA processing. ISTA requires many times more iterations to reach the FISTA result. For CONTIN, MaxEnt, and ITAMeD, the positions and amplitudes of the peaks were estimated with the Gaussian fitting (MATLAB script Gaussian curve f it23). In the case of NNLS, peak parameters were estimated from the position of the highest peaks. The obtained distributions of diffusion coefficients for a noise level of 0.1% (generated for each run independently) are presented in Figure 2. The full results of simulations can be found in the SI. Tests for Component Resolution. The important aspect of the processing method is its ability to resolve diffusion coefficients that do not differ significantly. We tested the resolution of all mentioned methods by performing series of simulations with pair of peaks of equal amplitude, but varying difference of diffusion coefficients. We changed the relative diffusion coefficients ratio from 1:1.5 to 1:3. Each simulation was performed with additional Gaussian noise on the level of 0.1% of the first data point. The results of the simulation are shown in Figure 3. In addition, to demonstrate the possible negative aspects of sparsity enforcing methods, we simulated the ITAMeD reconstruction of a “broad” diffusion coefficient distribution having Gaussian shape. The expected diffusion profile and the result of ITAMeD calculations are shown in Figure 4. Experiments. We tested our approach on experimental data of polyethylene glycol (PEG) polymers (Agilent Technologies, United Kingdom), solutions, and mixtures of these with various molecular weights (1080 g/mol, 11 840 g/mol, and 124 700 g/ mol) of a narrow molecular weight distribution (polydispersity index of 1.05) and compared the results with those from

Figure 5. Comparison between the (A) ITAMeD, (B) CONTIN, and (C) MaxEnt for the experimental data: three-component mixture of PEGs with molecular weights of 1080 g/mol, 11 840 g/mol, and 124 700 g/mol (red line); single-component solution of PEG 124 700 g/ mol (gray line); a single-component solution of PEG 11 840 g/mol (dark blue line); and a single-component solution of PEG 1080 g/mol (light blue line). The dotted line represents the reference data obtained by the monoexponential LMF fitting of single-component solutions with the following diffusion coefficients: 1.6 × 10−11 m2/s, with a sum of squared residuals (ssq) of 1.1 × 10−3 for PEG124700; 6.3 × 10−11 m2/s, and ssq = 2.5 × 10−4 for PEG11840; 2.3 × 10−10 m2/s, and ssq = 2.8 × 10−4 for PEG1080.

The disadvantage of ISTA is a slow convergence rate, which may be particularly inconvenient in the case of ill-conditioned problems, such as ILT. As one of the solutions, the fast iterative shrinkage-thresholding algorithm (FISTA) was proposed.18 In the present paper, we propose the application of S1 norm regularization with FISTA to process experimental PGSE datasets. We also present the comparison of the method referred to as Iterative Thresholding Algorithm for Multiexponential Decay (ITAMeD) with commonly applied algorithms.

■

METHODS The performance of ITAMeD was evaluated in the simulations and experiments described below. For the detailed description of the algorithm and the MATLAB19 code with example data, see the Supporting Information (SI). 1831

dx.doi.org/10.1021/ac3032004 | Anal. Chem. 2013, 85, 1828−1833

Analytical Chemistry

Article

Figure 6. Comparison between the CONTIN (green dashed line), MaxEnt (blue dotted line), and ITAMeD (red line). The processed experimental dataset is a mixture of all three PEGs in the concentration ratio of 3:1:2. Reference diffusion coefficients marked with dotted lines were obtained by monoexponential fitting of signals obtained for single-compound samples.

result.24 Our simulation suggests that a similar effect can be observed in MaxEnt reconstruction as well. On the other hand, NNLS succeeds in resolving components, even for diffusion coefficient ratios as low as 1:1.5. Unfortunately, it is also heavily disturbed by noise, which is clearly visible in Figure 2. All the approaches are comparable in terms of numerical efficiency. Processing time was on the order of seconds on the Linux 64bit personal computer (PC) with an Intel i5 3.3 GHz processor. We found that the PEG mixtures are well-suited for the tests of multiexponential fitting methods, because all of the protons in PEG molecules, independent of molecular weight, resonate at the same frequency. Thus, the only way of distinguishing between them using NMR is to perform a PGSE experiment and resolve their diffusion coefficients. The experiments show that both MaxEnt and ITAMeD successfully recover the diffusion coefficients of PEG mixture, reproducing values obtained from measurements of singlecomponent samples, as shown in Figure 5. For CONTIN, the deviations are slightly higher. Note that, for the PEG124700 sample, the diffusion profile obtained by ITAMeD is, in fact, a two-component profile, with a small peak at ∼9 × 10−10 m2/s. This is consistent with the result of LMF monoexponential fitting, where the sum of squared residuals is significantly higher than that in the case of PEG1080 and PEG11840 (see the caption for Figure 5). The effect can be thus explained by the presence of rapidly diffusing compounds, probably due to polymer hydrolysis. Notably, other methods, i.e., CONTIN and MaxEnt, do not give a significant “second peak” for PEG124700, but the peak is slightly shifted toward higher diffusion coefficients. Similarly, as in the simulations, the MaxEnt disturbs the magnitude of the weakest component of the 3:1:2 sample, while ITAMeD successfully establishes all three components (Figure 6). For CONTIN, the weakest component appears, but at a significantly shifted position, which is also consistent with the simulation results. The good resolution of the ITAMeD result allows easy peak integration, for comparison of the component magnitudes. For the three-component mixture from Figure 5, we obtained

MaxEnt and CONTIN processing. The measurements were acquired for each PEG solution and mixtures of all PEGs with 0.1% (w/w) of each (see Figure 5). We tested the mixture with different ratios (relative concentration of 3:1:2, PEG124700: PEG11840: PEG1080) which is presented in Figure 6. 1H PGSE NMR experiments were carried out on a 600 MHz Bruker spectrometer (Bruker, Germany) equipped with a Model Diff30 diffusion probe with GREAT40 amplifiers, providing a maximum gradient strength of 12 T/m. A 5-mm sample tube was used with a sample volume of 400 μL. A stimulated spin-echo sequence, including spoilers, was used with Δ = 50 ms and a gradient pulse duration δ of 2 ms. Gradient strengths were varied exponentially in 64 steps from 0.07 to 5 T/m, and the number of scans was 16 (see SI for details). All measurements were performed at 25 °C, using the PEG signal at 3.7 ppm for evaluation.

■

RESULTS AND DISCUSSION Simulated and experimental data processed with aforementioned algorithms allowed us to draw conclusions on the performance of ITAMeD, in comparison to other approaches. Only univariate approaches were used for comparison, since these were the most relevant for a signal that has only one component in a frequency dimension. First, the simulation shows that NNLS, CONTIN, and LMF methods are much more prone to instability caused by noise than MaxEnt and ITAMeD (Figure 2 and Tables 1 and 2 in the SI ). Also, it seems that ITAMeD is the most successful among the methods in recovering the position and magnitude of the weakest component. For other components, the reconstruction accuracy (standard deviation of peak position) is comparable to that obtained by CONTIN and MaxEnt only in the low noise case (0.1%). For higher noise levels, it is slightly worse; nevertheless, the total result seems to be superior, i.e., the error for all three peaks is equilibrated and acceptably low. Both CONTIN and MaxEnt successfully recover the diffusion coefficients of the two strongest peaks, but the amplitude of the weakest component seems to be disturbed. This problem is related to the effect known as “oversmoothing”. CONTIN has been already reported to oversmooth the 1832

dx.doi.org/10.1021/ac3032004 | Anal. Chem. 2013, 85, 1828−1833

Analytical Chemistry

Article

(2) Callaghan, P. T. Translational Dynamics and Magnetic Resonance: Principles of Pulsed Gradient Spin Echo NMR; Oxford University Press: Cambridge, MA, 2011. (3) Stejskal, E. O.; Tanner, J. E. J. Chem. Phys. 1965, 42, 288. (4) Johnson, C. S. Prog. Nucl. Magn. Reson. Spectrosc. 1999, 34, 203− 256. (5) Antalek, B. Concepts Magn. Reson. 2002, 14, 225−258. (6) Nilsson, M.; Morris, G. A. Anal. Chem. 2008, 80, 3777−3782. (7) Van Gorkom, L. C. M.; Hancewicz, T. J. Magn. Reson. 1998, 130, 125−130. (8) Press, W. H.; Teukolsky, S.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: New York, 1992. (9) Provencher, S. Comput. Phys. Commun. 1982, 27, 229−242. (10) Roush, F. W. Math. Social Sci. 1984, 7, 298−300. (11) Delsuc, M. A.; Malliavin, T. E. Anal. Chem. 1998, 70, 2146− 2148. (12) Tikhonov, A. N. Soviet Math. Dokl. 1963, 4, 1035−1038. (13) Lawson, C. L.; Hanson, R. J. In Solving Least Squares Problems, 3rd Edition; Lawson, C. L., Hanson, R. J., Eds.; Prentice−Hall Series in Automatic Computation; Prentice−Hall: Englewood Cliffs, NJ, 1995; pp 1−337. (Originally published in 1974.) (14) Slawski, M.; Hein, M. http://arxiv.org/abs/1205.0953 2012, (15) Ramlau, R.; Teschke, G. Numer. Math. 2006, 104, 177−203. (16) Zarzer, C. Inverse Problems 2009, 25, 025006. (17) Daubechies, I.; Defrise, M.; De Mol, C. Commun. Pure Appl. Math. 2004, 57, 1413−1457. (18) Beck, A.; Teboulle, M. SIAM J. Imaging Sci. 2009, 2, 183−202. (19) MATLAB, version 7.12.0 (R2011a); The MathWorks, Inc.: Natick, MA, 2011. (20) Marino, I.-G. File ExchangeMATLAB Central, 2004. (21) Tramesel, D.; Catherinot, V.; Delsuc, M.-A. J. Magn. Reson. 2007, 188, 56−67. (22) Balda, M. File ExchangeMATLAB Central, 2007. (23) Sivan, Y. File ExchangeMATLAB Central, 2006. (24) Stock, R. S.; Ray, W. H. J. Polym. Sci.: Polym. Phys. Ed. 1985, 23, 1393−1447.

relative magnitudes of 0.92:0.85:1 (PEG124700:PEG11840:PEG1080), whereas for the mixture from Figure 6, it was 3:0.75:1.76. The observed enhanced resolution of the ITAMeD processing can be explained as a result of using the sparsitypromoting S1 norm in a penalty function. Nevertheless, the sparsity-promoting approach reduces the number of potential applications of the method to the samples with “discrete” diffusion profiles. For instance, it may disturb the reconstruction in the case of a polydisperse sample, where a continuous distribution would be wrongly represented as a set of “peaks” (see, e.g., the two “peaks” in Figure 4). Thus, one should not draw conclusions concerning the polydispersity of a sample from distribution shapes obtained by ITAMeD, since they can be artificially narrowed, because of the sparsity-enforcing constraint. As mentioned in the beinning of the paper, we consider the concept of ITAMeD as a hybrid idea combining (implicit) sparsity assumption (as in LMF) and the lack of assumption on the number of components (CONTIN and MaxEnt). Thus, the method may be considered as an alternative for CONTIN/ MaxEnt processing of PGSE data optimal for monodispersed samples.

■

CONCLUSION We have discussed the advantages and disadvantages of the fast iterative shrinkage thresholding algorithm (FISTA) applied to NMR diffusometry. The simulations and experiments, as well as theoretical background, suggest that ITAMeD is a method of choice in the case of monodispersed samples, i.e., those with diffusion coefficient distribution dominated by zeros. However, sparsity-enforcing procedures may disturb the broad diffusion coefficient distributions of polydispersed samples.

■

ASSOCIATED CONTENT

S Supporting Information *

This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by Bio-NMR Project No. 261863, funded by European Commission’s Framework Program 7 (FP7), and by VINNOVA, through the VINN Excellence Centre SuMo Biomaterials (Supramolecular Biomaterials Structure dynamics and properties). M.U., W.K., and K.K. thank the Foundation for Polish Science for the support with TEAM Programme. K.K. thanks the Polish Ministry of Science and Higher Education for the support with Iuventus-Plus Programme. The Swedish NMR Centre is acknowledged for spectrometer time.

■

REFERENCES

(1) Price, W. S. NMR Studies of Translational Motion: Principles and Applications (Cambridge Molecular Science); Cambridge University Press: New York, 2009. 1833

dx.doi.org/10.1021/ac3032004 | Anal. Chem. 2013, 85, 1828−1833