Improving the Resolution of MPM Recovered Relaxometry Parameters with Proper Time Domain Sampling

Abstract

The matrix pencil method (MPM) is a powerful tool for processing transient nuclear magnetic resonance (NMR) relaxation signals with promising applications to increasingly complex problems. In the absence of signal noise, the eigenvalues recovered from an MPM treatment of transient relaxometry data reduce to relaxation coefficients that can be used to calculate relaxation time constants for known sampling time ∆t. The MPM eigenvalue and relaxation coefficient equality as well as the resolution of similar eigenvalues and thus relaxation coefficients degrade in the presence of signal noise. The relaxation coefficient ∆t dependence suggests one way to improve MPM resolution by choosing ∆t values such that the differences between all the relaxation coefficient values are maximized. This work develops mathematical machinery to estimate the best ∆t value for sampling damped, transient relaxation signals such that MPM data analysis recovers a maximum number of time constants and amplitudes given inherent signal noise. Analytical and numerical reduced dimension MPM is explained and used to compare computer-generated data with and without added noise as well as treat real measured signals. Finally, the understanding gleaned from this effort is used to predict the best data sampling time to use for non-discrete, distributions of relaxation variables.

1 Introduction

Nuclear magnetic resonance (NMR) spectroscopy has a rich history in the study of solid, liquid, and gas phase samples in chemistry, biology, physics, medicine, etc [1]. High-resolution NMR, typically performed at high magnetic field with superconducting magnets using chemical shifts, scalar J couplings [2, 3], and long-range dipolar couplings spectrally manifesting as sample-ordered line-splitting or relaxation-induced line broadening [4,5,6], has emerged as the premiere method to study three-dimensional chemical structure in solid and liquid phases. Lower magnetic field bench top, permanent magnets are beginning to become more popular in these high-resolution studies [7] although they have been used for years in magnetic resonance imaging and relaxometry [8, 9]. The imaging of materials [10] and living objects [11] is an entirely separate active research area that provides critical information regarding the preparation of new devices with engineered properties and of course, in the case of humans, insight into health and illness. NMR relaxometry, often performed on protonated samples with Larmor frequencies in the 1–20 MHz range, commonly uses low-resolution permanent magnets [8, 9]. Here the spectral resolution is poor, chemical shifts and scalar J couplings are not resolved, and typically the spin–lattice, longitudinal and spin–spin, transverse relaxation times, and the macroscopic diffusion coefficient as well as their respective two-dimensional correlations are the core measured parameters. The application of these low-field relaxometry measurements to complex mixtures often reveals several relaxation times, and, therefore, a unique signal that is described by a weighted sum of exponentials [12, 13]. Provided these parameters can be measured, and more importantly extracted from the time domain, transient signals, they can be used as a fingerprint for a given substance or mixture, much like the chemical shifts and J couplings form fingerprints for complex macromolecules in high-resolution NMR work [7].

The Fourier transform (FT) used to extract frequencies and amplitudes from oscillatory high-resolution time domain NMR-free induction signals is not capable of recovering relaxation time values from exponentially damped, non-oscillatory NMR relaxometry signals. As relaxometry signals are not complex and oscillatory, the FT produces a spectrum with one broad peak at zero frequency. In these situations, the inverse Laplace transform (ILT) has emerged as the current industry standard for recovering relaxation time values from multicomponent NMR relaxometry data [13]. The ILT takes exponentially damped, real, and non-oscillatory transient signals as input and provides distributions of relaxation rates or times as output. Discretization of the ILT for application to data is accomplished by routines such as the Lawson–Hanson algorithm [14], a non-negative least squares (NNLS) method. As is common for inversion problems, ILT is an inherently ill-posed problem when unmodified. The consequences of this ill-posed nature are that several solutions may satisfy the problem, and that these solutions are easily influenced by minute differences in the data. Using a set of pure exponential decay functions to model a data set that is a mixture of noise with pure exponential decay functions results in severe ILT instability. To restrict the number of possible solutions, L2 regularization is often applied at the expense of increased computation time and of course added artificial output distribution broadening [15].

In recent years, some alternative methods have found success in bypassing some of the stability complications of the L2 regularized NNLS algorithms. One such method is Modified Total Generalized Variation regularization which conjoins L1 and L2 regularizations allowing variability in the discernment of discrete and continuous distribution features [16]. Moving away from the NNLS problem, the q-exponential non-linear least squares method offers a statistics-driven approach to multiexponential modeling [17].

It is for these same reasons that the matrix pencil method (MPM) was developed and used to treat a wide range of relaxometry data from the study of textbook two-component samples with different spin–lattice and spin–spin relaxation times to unique emerging materials with new properties and vaccine-binding biomolecules [18]. The work builds on earlier uses of the MPM in high-field solid-state NMR problems [19], speech analysis [20], and remote radar sensing [21]. The MPM is an algebraic way to treat transient, exponential signals. As it is not a mathematical transformation, the MPM sidesteps some of the complications encountered when using ILT, providing clear, discrete solutions that make it well-suited for NMR relaxometry studies.

Initial work in this area was designed to advertise the usefulness of MPM to the broader NMR community by including applications from a wide range of disciplines as well as measurements of the core relaxation parameters [18]. A more recent, purely experimental paper considered the ultimate resolving capability of the conventional MPM algorithms introduced in the initial work [22]. Here two separate sample containers housing water with different paramagnetic impurity concentrations and known spin–lattice relaxation times were simultaneously placed on top of a commercially available Magritek PM-25 single-sided NMR instrument. The analysis of the two-component saturation recovery transient signals revealed that the resolving power of MPM was roughly twice that of ILT. In some cases, operation with MPM offered even more gain over ILT. However, in these cases, perfect amplitude fidelity or reproduction of accurate fractional relaxation coefficient contributions to the signal was sacrificed.

It was recognized in that recent experimental work that signal noise is the primary factor limiting MPM resolution. The greater the signal noise, the lower the resolution. The work reported here and described below acknowledges that the MPM reports eigenvalues that, in the absence of noise, are relaxation coefficients λ_m that depend on the inherent relaxation rate and signal sampling time ∆t. Thus, the ability to resolve two relaxation times or rates reduces to resolving two relaxation coefficients that depend on the choice of ∆t. This fact suggests that there should be an optimum ∆t value that provides the best chance of resolving relaxation coefficients with similar relaxation time values. The next section describes how to maximize the sum of square difference of relaxation coefficients SS_λ to obtain the best data sampling time ∆t_max and motivates a performance metric or noise tolerance N_T that can be used to verify the ∆t_max value theoretically as well as to compare to actual experimental measurements.

2 Theory

The application of both one- and two-dimensional MPM to NMR relaxometry data is described in greater detail elsewhere [18]. In the one-dimensional case considered here, the measured non-oscillatory, damped, relaxometry signal s((n – 1)∆t) sampled at the time (n–1)∆t serves as MPM algorithm input. All index counters like n, m, and q in the following equations are understood to begin at value one. Central to the MPM approach is representation of the signal in terms of a linear combination of damped exponential functions as

$$s((n - 1)\Delta t) = \sum\limits_{m = 1}^{mpts} {A_{m} \lambda_{m}^{n} } + N((n - 1)\Delta t)$$

(1)

where the relaxation coefficient is λ_m = exp(–R_m∆t), the relaxation rate R_m = 1/T_m is inversely proportional to the time constant T_m, the amplitude of the m^th relaxing signal component is A_m, and the total number of relaxing components is mpts. Standard Gaussian white noise represented by the zero average and cross-correlation, 〈N(n∆t)〉 = 〈N(n∆t)N(m∆t)〉 = 0, non-zero mean square, 〈N(n∆t)²〉 = 〈N²〉 ≠ 0, function N((n–1)∆t) is used as the noise source as described in more detail in the Experimental section. The two primary outputs from MPM signal analysis are an amplitude matrix $\tilde{V}$ and eigenvalues z_m. In the absence of noise where N((n–1)∆t) = 0 for all choices of n in Eq. (1), the MPM analysis disentangles the linear combination of exponentials to provide an estimate of the signal amplitudes rates from the diagonal elements of the amplitude matrix as A_m = V(m,m), and since λ_m = z_m, as R_m = log(z_m)/∆t.

In the presence of noise, the eigenvalues z_m deviate from their relaxation coefficient, noise-free values, λ_m. Not surprisingly, the closer the λ_m values, the larger the deviation of z_m from λ_m for fixed mean square noise amplitude 〈N²〉. The dependence of the relaxation coefficient values λ_m on the dwell time ∆t suggests that there will be a sampling time ∆t_max that leads to the largest separation between λ_m values. This time can be determined by maximizing the sum of square differences as a function of ∆t

$$SS_{\lambda } = \sum\limits_{m = 1}^{mpts} {(\lambda_{m} - \overline{\lambda })^{2} } \,\,\,\,\,\,\,\,\,\,with\,\,\,\,\,\,\,\,\,\overline{\lambda } = \frac{1}{mpts}\sum\limits_{m = 1}^{mpts} {\lambda_{m} } ,$$

(2)

an effort that will provide the least sensitivity of the z_m values to noise. A useful way to study the resolution limits of the MPM for both theoretical and real experimental data is to compare the sum of square differences obtained in the noise-free case shown in Eq. (2) to the sum of squares based purely on MPM eigenvalues with added noise

$$SS_{z} = \sum\limits_{m = 1}^{mpts} {(z_{m} - \overline{z})^{2} } \,\,\,\,\,\,\,\,\,\,with\,\,\,\,\,\,\,\,\,\overline{z} = \frac{1}{mpts}\sum\limits_{m = 1}^{mpts} {z_{m} }$$

(3)

using the noise tolerance

$$N_{T} = \frac{{SS_{z} - SS_{\lambda } }}{{SS_{z} + SS_{\lambda } }}\, \times \,100$$

(4)

that describes the percent difference between SS_z and SS_λ for many values of ∆t. In practice this would be accomplished by determining SS_λ and SS_z at different dwell times ∆t. This is most easily done by first performing MPM on the raw transient sampled at the dwell time ∆t. Subsequent MPM analyses are then performed on signals generated from the raw transient by taking every other data point at 2∆t, every third data point at 3∆t and so on as described in Fig. 1a. In the noise-free case SS_z = SS_λ and N_T = 0 for all q∆t values. As the noise increases, SS_z deviates from SS_λ, N_T becomes non-zero, and the ability of MPM eigenvalues to reproduce relaxation coefficients degrades. In the case of a purely theoretical comparison, one knows the λ_m values from knowledge of the R_m rates and multiples of the dwell time q∆t enabling calculation of SS_λ. The value of SS_z is obtained from the z_m eigenvalues generated from MPM analysis of a decaying transient signal with added random noise having the same dwell time multiple q∆t. For experimental comparisons, samples with known R_m values in separate containers housed in the detection coil are used. Since the R_m values are known from separate measurements of each container alone, SS_λ can be calculated for each multiple of the dwell time q∆t. The MPM analysis of the bulk signal obtained from all the containers simultaneously placed in the magnet and detection coil provides the z_m eigenvalues at each q∆t needed to calculate SS_z and thus a value for N_T.

Fig. 1

The decay signal shown in (a) as the solid line helps graphically identify the signal points used for the MPM analysis. For example, for q = 1, the points separated by ∆t and indicated with a solid circle, open circle, square, and triangle are used, while for q = 2, the points separated by 2∆t and indicated with a solid circle, open square, open diamond and the s(6∆t) point not shown are used. The plot shown in (b) for two relaxation coefficients shows how the shaded gray region representing noise suggests that the points at ∆t and 4∆t shown as “x” symbols will present large N_T values while the open symbols at 2∆t and 3∆t yield smaller N_T values with the smallest being at ∆t_max. The shaded gray region bounds are described by the f ₊〈N²〉 and f_–〈N²〉 noise terms in the two relaxation coefficient eigenvalues shown in Eq. (12)

To understand the meaning of and how to develop and use SS_λ, SS_z, and N_T, consider the special signal involving just two relaxing components, or equivalently a bi-exponential decay signal where λ₁ = exp(–R₁∆t) and λ₂ = exp(–R₂∆t) for the respective decay rates R₁ and R₂. It should be clear that there are many choices of λ₁ and λ₂ given the dwell time ∆t and that the subscript numbers indicate different relaxation rates not spin–spin versus spin–lattice relaxation rates. Ignoring the shaded gray region of Fig. 1b that will be described below, the solid, curved line in the northwest corner of the plot shows the behavior of λ₁ and λ₂ for fixed R₁ and R₂ values as a function of variable dwell time ∆t. Recall that this line describes the eigenvalues z_+,- in the absence of noise as z_+,- = λ_1,2 in this case. The model transient signal shown in Fig. 1a obtained with a fixed ∆t can be used to explore the solid line curve in Fig. 1b created with a continuous ∆t variable by taking each point at ∆t, every other point at 2∆t, every third point at 3∆t, every fourth point at 4∆t, etc. When operating in this way, the open circle, “x”, and “*” symbols shown on the solid line indicate the λ_1,2 pairs obtained at the chosen q∆t value. Here q ≥ 1 is an index. The dot-dashed line along the diagonal in the plot represents the condition where λ₁ = λ₂. In the absence of noise, λ₁ and λ₂ are not resolved because they are the same. As all the five symbolled points on the solid curve in Fig. 1b for different q∆t are displaced from the λ₁ = λ₂ diagonal, they are in principle able to be resolved in the absence of noise. However, the best resolution will occur for the largest difference between λ₁ and λ₂ or for the ∆t value obtained from maximizing SS_λ shown in Eq. (2) above. For two components, SS_λ reduces to the difference λ₁–λ₂ which is captured by the dotted line along the anti-diagonal in Fig. 1b. When maximized as a function of ∆t, the difference yields the optimized time ∆t_max = log(R₂/R₁)/(R₂–R₁) and the most distinguishable relaxation coefficient $\lambda_{1} {}^{\max}=(\text{R}_{1}/\text{R}_{2})^{\text{R}_{1}/(\text{R}_{2}-\text{R}_{1})}$ and $\lambda_{2} {}^{\max}=(\text{R}_{1}/\text{R}_{2})^{\text{R}_{2}/(\text{R}_{2}-\text{R}_{1})}$. Graphically in the two-component case, the ∆t_max value obtained from maximizing SS_λ corresponds to finding the longest vector perpendicular to the λ₁ =  λ₂ diagonal and between the diagonal and the λ_1,2 curve.

In the absence of noise, any two λ_m values are always completely resolved, although the best resolution occurs at ∆t_max with the associated λ₁^max and λ₂^max values as just described. Moreover, the parameter N_T in Eq. (3) is always zero in the absence of noise as SS_z = SS_λ. To determine the effect of added noise in the case of two relaxing components, MPM can be accomplished analytically, an approach that avoids the uncertainties in addition to the mysteries of the singular valued decomposition included in most MPM algorithms [15, 16]. Here one first uses the measured data array s((n-1)∆t) sampled at the times (n-1)∆t to construct two data matrices from just four points.

$$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{1} = \left[ {\begin{array}{*{20}c} {s(0)} & {s(q\Delta t)} \\ {s(q\Delta t)} & {s(2q\Delta t)} \\ \end{array} } \right]\,\,\,\,\,\,\,and\,\,\,\,\,\,\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{2} = \left[ {\begin{array}{*{20}c} {s(q\Delta t)} & {s(2q\Delta t)} \\ {s(2q\Delta t)} & {s(3q\Delta t)} \\ \end{array} } \right].$$

(5)

Notice that if q = 1, the first four points in the measured data s((n-1)∆t) array are used and λ_m = exp(-R_m∆t). For q = 3, every third data point is used to select the {s(0), s(3∆t), s(6∆t), s(9∆t)} four points for analytical MPM analysis as described in Fig. 1a and λ_m = exp(− 3R_m∆t). The matrix pencil equation

$$(\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{1}^{ - 1} \cdot \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_ 2 - z_{m} \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{E} ) \cdot \vec{p}_{m} = 0$$

(6)

in terms of the unit matrix $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{E}$ and the eigenvectors $\vec{p}_{m}$ is then solved for the eigenvalues z_m. Essentially one finds the inverse matrix

$$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{1}^{ - 1} = \frac{1}{{s(0)s(2q\Delta t) - s(q\Delta t)^{2} }}\,\,\left[ {\begin{array}{*{20}c} {s(2q\Delta t)} & { - s(q\Delta t)} \\ { - s(q\Delta t)} & {s(0)} \\ \end{array} } \right]$$

(7)

and then the eigenvalues of the matrix product

$$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{1}^{ - 1} \cdot \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{2} = \left[ {\begin{array}{*{20}c} 0 & {\frac{{s(2q\Delta t)^{2} - s(q\Delta t)s(3q\Delta t)}}{{s(0)s(2q\Delta t) - s(q\Delta t)^{2} }}} \\ 1 & {\frac{s(0)s(3q\Delta t) - s(q\Delta t)s(2q\Delta t)}{{s(0)s(2q\Delta t) - s(q\Delta t)^{2} }}} \\ \end{array} } \right]$$

(8)

in the usual way as

$$z_{ \pm } = \frac{1}{2}\frac{s(0)s(3q\Delta t) - s(q\Delta t)s(2q\Delta t)}{{s(0)s(2q\Delta t) - s(q\Delta t)^{2} }} \pm \sqrt {\frac{1}{4}\left( {\frac{s(0)s(3q\Delta t) - s(q\Delta t)s(2q\Delta t)}{{s(0)s(2q\Delta t) - s(q\Delta t)^{2} }}} \right)^{2} + \frac{{s(2q\Delta t)^{2} - s(q\Delta t)s(3q\Delta t)}}{{s(0)s(2q\Delta t) - s(q\Delta t)^{2} }}} .$$

(9)

These analytical expressions for MPM analysis can be used to explore the effect of signal noise on the z_± eigenvalues and their deviation from λ₁ and λ₂. For two relaxing components with amplitudes A₁ and A₂, Eq. (1) implies that four consecutive signal points with noise beginning at t = 0 are

$$\begin{gathered} \,\,\,\,\,\,\,\,\,\,s(0) = A_{1} + A_{2} + N(0)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,s(2q\Delta t) = A_{1} \lambda_{1}^{2} + A_{2} \lambda_{2}^{2} + N(2q\Delta t) \hfill \\ s(q\Delta t) = A_{1} \lambda_{1} + A_{2} \lambda_{2} + N(q\Delta t)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,s(3q\Delta t) = A_{1} \lambda_{1}^{3} + A_{2} \lambda_{2}^{3} + N(3q\Delta t). \hfill \\ \end{gathered}$$

(10)

Taylor expansion of Eq. (9) to second order in terms of these four independent noise variables along with knowledge that the average noise and cross-correlation are zero, 〈N(n∆t)〉 = 〈N(n∆t)N(m∆t)〉 = 0, and with non-zero mean square, 〈N(n∆t)²〉 = 〈N²〉 ≠ 0, results in two noise factors

$$\begin{aligned} f_{ + } =\, & [1 + A_{1} A_{2} (\lambda_{1}^{4} + 4\lambda_{1}^{2} + \lambda_{2}^{2} + \lambda_{1}^{4} \lambda_{2}^{2} + 4\lambda_{1}^{3} \lambda_{2} + 4\lambda_{1}^{2} \lambda_{2}^{2} + 4\lambda_{1} \lambda_{2} ) \\ & - A_{2}^{2} (2\lambda_{2}^{4} + 2\lambda_{2}^{2} - \lambda_{1}^{2} + \lambda_{1} \lambda_{2}^{5} - \lambda_{1}^{2} \lambda_{2}^{4} + 2\lambda_{1} \lambda_{2}^{3} - 4\lambda_{1}^{2} \lambda_{2}^{2} - \lambda_{1} \lambda_{2} )]/\left( {\lambda_{1} - \lambda_{2} } \right)^{5} \\ & f_{ - } = [ - 1 - A_{1} A_{2} (\lambda_{2}^{4} + 4\lambda_{2}^{2} + \lambda_{1}^{2} + \lambda_{1}^{2} \lambda_{2}^{4} + 4\lambda_{1} \lambda_{2}^{3} + 4\lambda_{1}^{2} \lambda_{2}^{2} + 4\lambda_{1} \lambda_{2} ) \\ & + A_{1}^{2} (2\lambda_{1}^{4} + 2\lambda_{1}^{2} - \lambda_{2}^{2} + \lambda_{1}^{5} \lambda_{2} - \lambda_{1}^{4} \lambda_{2}^{2} + 2\lambda_{1}^{3} \lambda_{2} - 4\lambda_{1}^{2} \lambda_{2}^{2} - \lambda_{1} \lambda_{2} )]/\left( {\lambda_{1} - \lambda_{2} } \right)^{5} \\ \end{aligned}$$

(11)

that can be used to approximate the two eigenvalues to first order in 〈N²〉 as

$$z_{ + , - } \simeq \lambda_{1,2} + f_{ + , - } \langle N^{2} \rangle .$$

(12)

These two equations are interesting because they suggest that the added signal noise does not directly yield eigenvalues that fluctuate equally about an average z_+,-  = λ_1,2. Rather, the added signal noise induces a shift of the z_+,- eigenvalues away from the noise-free λ_1,2 values. It is this shift that is captured by the shaded gray region in Fig. 1b with bounds indicated by dashed black lines that correspond to the f₊〈N²〉 and f_{_}〈N²〉 factors appended to the λ_1,2 relaxation coefficients in the z_+,- eigenvalues in Eq. (12). A larger average signal noise 〈N²〉 will generate a larger shaded gray region on this plot. Consideration of the eigenvalue difference, z₊ − z_– =  λ₁ − λ₂ + (f₊ − f_–)〈N²〉, indicates that the λ₁–λ₂ noise-free value increases to a larger 〈N²〉-dependent number where the increase calculated from Eq. (11), (f₊ − f_–)〈N²〉, can be used to determine the noise tolerance defined in Eq. (4) as

$$N_{T} = \frac{{(f_{ + } - f_{ - } )\langle N^{2} \rangle }}{{2(\lambda_{1} - \lambda_{2} ) + (f_{ + } - f_{ - } )\langle N^{2} \rangle }}\, \times \,100.$$

(13)

It is important to remember that the N_T value directly depends on q∆t through the definition of the λ_1,2 values and that the best separation of the z_+,- eigenvalues will be when q∆t = ∆t_max and when N_T is minimum. In terms of the plot in Fig. 1b, the N_T value for the points at ∆t and 4∆t labeled with an “x” and firmly embedded in the shaded region will have large, close to 100% N_T values because the noise is so large that the relaxation components are not resolved. The open circles at 2∆t and 3∆t will have smaller N_T values. However, the point for ∆t_max labeled with an “*” will display the minimum N_T value and thus be most easily resolved in the presence of noise. As long as the ∆t_max point on the λ_1,2 curve is greater than the system noise or equivalently lies outside of the gray shaded region in Fig. 1b, the N_T value will approach zero and the two components will likely be resolved. If, however, this same point falls within the shaded gray area, a case not shown here, the N_T value will become much greater and a reliable estimate of the λ_1,2 values is not possible given the inherent signal noise.

A similar graphical interpretation of SS_λ drives the determination of ∆t_max for signals involving more than just two components. In the three-component case, the three λ₁ = exp(-R₁q∆t), λ₂ = exp(-R₂q∆t), and λ₃ = exp(-R₃q∆t) relaxation coefficients for the three R₁, R₂, and R₃ rates label three orthogonal axes and maximization of SS_λ amounts to determining the maximum perpendicular distance from the body diagonal to the λ_1,2,3 curve. Noise in this case broadens the body diagonal line to the surrounding volume and only certain ∆t values will provide λ_1,2,3 points outside of this volume. These ∆t sampling rates produce small N_T values and thus faithful experimental estimates of the true relaxation properties can be recovered. The effect of the added noise can be mathematically considered in the same way it was for two relaxation coefficients above. The additional λ₃ contribution to the signal s(t) requires an additional two signal points s(4q∆t) and s(5q∆t) to accomplish reduced dimensional MPM analysis. In this case, the $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{1}$ and $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{2}$ matrices are constructed as

$$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{1} = \left[ {\begin{array}{*{20}c} {s(0)} & {s(q\Delta t)} & {s(2q\Delta t)} \\ {s(q\Delta t)} & {s(2q\Delta t)} & {s(3q\Delta t)} \\ {s(2q\Delta t)} & {s(3q\Delta t)} & {s(4q\Delta t)} \\ \end{array} } \right]\,\,\,\,\,\,\,\,\,\,\,\,\,\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{2} = \left[ {\begin{array}{*{20}c} {s(q\Delta t)} & {s(2q\Delta t)} & {s(3q\Delta t)} \\ {s(2q\Delta t)} & {s(3q\Delta t)} & {s(4q\Delta t)} \\ {s(3q\Delta t)} & {s(4q\Delta t)} & {s(5q\Delta t)} \\ \end{array} } \right]$$

(14)

and the $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{1}^{ - 1} \cdot \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{2}$ eigenvalues z_1,2,3 are calculated to find λ_1,2,3 and thus R_1,2,3 given q∆t. Although analytical solutions for the z_1,2,3 eigenvalues were generated akin to the two-component case shown in Eq. (12), the overwhelming complexity and number of terms prevent reproduction here. Instead, numerical diagonalization of $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{1}^{ - 1} \cdot \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{2}$ was used to calculate N_T for MPM algorithm performance estimates on purely theoretical and experimental data.

Extension of this approach to more than three relaxing components is straightforward. The MPM analysis of signals with four relaxing components requires eight signal point and five relaxing components requires ten signal points. A signal with m relaxing components will thus require 2m points for treatment with reduced MPM, and of course numerical matrix diagonalization is required.

3 Experimental

Gadopentetic acid (Gd-DTPA) and sodium chloride were obtained from Sigma-Aldrich and used without further purification. In-house deionized water was used to prepare a 0.9% (by mass) sodium chloride stock solution. Six samples in the 0.05 mM < [Gd-DTPA] < 2 mM concentration range were prepared from this stock solution and are listed in Table 1.

Table 1 Summary of water T₂ values as a function of Gd-TPA concentration

All NMR measurements were accomplished using a Tecmag Apollo controlled, Aspect Instruments M100, 1.01 T, magnetic resonance imaging spectrometer operating at a 43.7 MHz ¹H Larmor frequency. All MPM analyses of experimental data were accomplished on transient, exponentially damped signals obtained with the standard Carr–Purcell–Meiboom–Gill (CPMG) pulse sequence [23]. The time between the 1024 consecutive collected spin echoes was 10 ms, the π/2 pulse time with 27 W of applied power was 35 μs and multiple samples easily fit within the 6 cm diameter, 10 cm long NMR detection coil. Signal averaging was used to sum eight separate signal acquisitions and this averaged signal corresponds to one experimental trial. The experimental data shown in all figures correspond to a 512-trial average. Experimental estimates of 〈N²〉 values were obtained by calculating the mean square of the last 200 points of the CPMG transient after scaling the raw signal so that the first point is one.

All data processing and simulations were accomplished using Matlab (Mathworks, Natick, MA). Because the measured transient CPMG signals from the six standard Gd-DTPA-doped samples were single-exponential, regression to exp(-(n-1)∆t/T₂) could be used to accurately determine T₂ values (R² > 0.95). The noise in numerical transient signal simulations was established by scaling random numbers Ɲ(0,1) extracted from a zero-centered, width of one, standard normal distribution by the root mean square signal noise as N((n-1)∆t) = Ɲ(0,1)〈N²〉^1/2 so that the mean square average and cross-correlation at each ∆t reduce to 〈N((n-1)∆t)²〉 = 〈N²〉 and zero. In all calculations, 15,000 noisy CPMG decay signals are used to calculate an average transient decay response that is used to estimate SS_z.

4 Results and Discussion

The goal of this work is to identify the data sampling time ∆t_max that maximizes the resolution of as many relaxation components λ_m as possible in the presence of real signal noise described by 〈N²〉. It is the maximization of SS_λ as a function of ∆t in Eq. (2) that identifies ∆t_max and it is the value of N_T in Eq. (4) as a function of q∆t that reports on the resolution performance. Here N_T values of zero are considered to perform well while N_T values exceeding zero perform less well. The theoretical results in Fig. 2(a–d) for two and three relaxation components reflect these comments.

Fig. 2

Purely theoretical predictions of the behavior of ∆t_max and N_T for two (a, b) and three (c, d) relaxation coefficient cases as a function of added noise. The vertical dashed line in all the plots indicates the value of ∆t_max obtained from maximizing SS_λ in Eq. (2) as a function of ∆t. For easy reference, the appropriate relaxation rate R_m values in units of Hz are shown in the small box in each plot. The solid lines shown in all plots correspond to numerically determined estimates of N_T. The open circles in (a) and (b) represent the two-component analytical result shown in Eq. (13). The two-coefficient comparison shown in (a) at the four 〈N²〉 = 10^–10 (green), 10^–8 (black), 10^–6 (red), and 10^–4 (blue) noise levels used transient signals s(t) with equal amplitude A₁ = A₂ = 0.5 exponential decay functions. One 〈N²〉 = 10^–4 noise value was used to generate (b) for four different exponential decay function amplitude pairs {A₁,A₂} = {0.5,0.5} (green), {0.25,0.75} (black), {0.33,0.67} (red), and {0.1,0.9} (blue). The three-coefficient comparison shown in (c) at the same four 〈N²〉 noise levels used transient signals s(t) with equal amplitude A₁ = A₂ = A₂ = 0.33 exponential decay functions. One 〈N²〉 = 10^–8 noise value was used to generate (d) for four different exponential decay function amplitude triples {A₁,A₂,A₃} = {0.33,0.33,0.33} (green), {0.25,0.25,0.5} (black), {0.2,0.4,0.4} (red), and {0.14,0.29,0.57} (blue). Some of the colors are difficult to see in (d) as the lines overlap (color figure online)

The open circles shown in Fig. 2(a, b) describe the behavior of Eq. (13) as a function of q∆t and 〈N²〉 for two relaxing components with different R₁ and R₂ values. The relevant R_m values in units of Hz are shown in the plot mounted boxes. The solid lines correspond to numerical diagonalization of Eq. (8) followed by averaging over the random number generated noise N(q∆t). The similarity of the analytical behavior of Eq. (13) to the numerical results in the neighborhood of ∆t_max indicated by the vertical dashed line, especially at very low 〈N²〉 value, suggests that it is safe to use numerical diagonalization followed by averaging to determine N_T values. The deviation of the analytical and numerical results displayed in Fig. 2a at large q∆t and 〈N²〉 values is not a limitation of numerical diagonalization. Rather it is due to restricting the eigenvalue expansion to the second-order term in the Taylor expansion used to develop the analytical eigenvalues in Eq. (12). Both the analytical and numerical approaches shown in Fig. 2b demonstrate that as the added noise 〈N²〉 increases the range of q∆t times that yield low to zero N_T values decreases and becomes centered on the ∆t_max value predicted from maximizing SS_λ. It is only when 〈N²〉 exceeds 10^–4 that N_T noticeably increases from zero implying that the added noise is too great to adequately recover meaningful λ_1,2 values from the z_+,- eigenvalues even though operation at ∆t_max is insured. In other words, as signal noise increases, the range of q∆t values providing z_+,- eigenvalues that faithfully report the λ_1,2 relaxation coefficient values decreases. In terms of the plot shown in Fig. 1b, as noise increases, the overall length of the solid line outside of the shaded gray area decreases and when 〈N²〉 exceeds 10^–4, the entire λ_1,2 curve is within the shaded gray area. The similarity of the solid line graphs with minima centered on ∆t_max shown in Fig. 2b demonstrates that relaxation coefficient amplitude has little effect on both ∆t_max and N_T, consistent with the determination of ∆t_max from SS_λ where the relative contribution of the relaxation components to the full signal is not considered.

Numerical diagonalization of the $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{1}^{ - 1} \cdot \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\thicksim}$}}{Y}_{2}$ matrix is used in all cases where more than two components are involved. The three-component case shown in Fig. 2c, like Fig. 2a, considers N_T as a function of q∆t and 〈N²〉 for three equal amplitude components with three different rates R_1,2,3 while Fig. 2d, like Fig. 2b, fixes 〈N²〉 and the R_1,2,3 values while varying the three amplitudes A_1,2,3. The three-component calculated results in Fig. 2(c, d) convey similar information as the two-component results in Fig. 2(a, b). As the noise 〈N²〉 increases, the range of q∆t values where N_T equals zero decreases and centers on the ∆t_max value predicted from maximizing SS_λ as a function of ∆t. There also appears to be very little dependence on the fractional contribution of each relaxation coefficient to the signal at fixed 〈N²〉 value, as judged from the similarity of the numerically generated curves in Fig. 2d. The important difference between the two- and three-component cases shown in Fig. 2 is sensitivity to noise. The three-component results always present N_T values much greater than the corresponding two-component case for the same applied 〈N²〉 value. For example, for 〈N²〉 = 10^-4 in Fig. 2a, N_T is ca. 5% at most, while in Fig. 2c, N_T exceeds 50%. This effect can also be observed in the comparison of Fig. 2(b, d) where similarly shaped curves are produced from drastically different applied 〈N²〉 values, here 10^-6 and 10^-8 respectively. The reason that the three-component case appears to be more sensitive to noise is directly related to the proximity of the chosen relaxation rates. The two-component case uses R₁ = 2.23 Hz and R₂ = 10.2 Hz while the three-component case includes the additional 4.36 Hz rate. It is the smaller difference of 2.13 Hz in the three-component case in comparison to the 7.96 Hz difference in the two-component case that likely leads to this sensitivity. If the added third rate led to a similar 7.96 Hz difference, then a similar noise sensitivity would be expected.

The set of six samples with the Gd-DTPA concentration in Table 1 was prepared, placed in six separate 1 mL plastic containers, and separately analyzed with the CPMG pulse sequence to obtain the listed T₂ time constants. The value of R reported in Table 1 is the inverse of the T₂ value. These samples were prepared to verify the predictions described in Fig. 2 by obtaining transient CPMG decay signals from different pairs, triples and quadruples of these tubes simultaneously placed inside of the NMR detection coil. The relaxation rates in Hz for the samples being simultaneously studied are included in all figure legends.

The two-component predictions shown in Fig. 2(a, b) are reflected in the experimental measurements. Here different pairs of the samples listed in Table 1 were simultaneously placed inside of the NMR detection coil and explored with the CPMG pulse sequence. The left-hand column in Fig. 3 shows how the N_T value changes when the relaxation properties of one of the containers are changed while keeping those for the second container fixed. The relaxation rates obtained from the separate individual containers in Hz are shown in the legends. As the relaxation rates of the separate samples converge, the relaxation coefficients λ_m at all q∆t values become more similar. Experimental results shown as the open circles agree well with the solid line theoretical predictions. The theoretical result is much smoother than the experimental data because the experimental results were generated from an average of 512, 8 scan transient signal trials with noise while the theoretical prediction involved averaging over 15,000 transient signals with the same noise amplitude as the measured data. The right-hand column in Fig. 3 displays an expected small variation in shape as a function of relaxation coefficient amplitude confirming the theoretical predictions summarized in Fig. 2b. The graphs presented in Fig. 3 suggest that the best MPM performance with noise is when q∆t = ∆t_max as the N_T value is minimum. All the N_T values for equal amplitude, similar relaxation coefficients shown in Fig. 3d and for drastically different amplitude relaxation coefficients shown in Fig. 3h are greater than zero. This observation suggests that even operation at ∆t_max will not reliably recover accurate relaxation coefficients and thus relaxation times given the experimental noise level.

Fig. 3

Comparison of experimental (open red circles) to theoretical (solid black line) N_T values for two-component signals. The vertical dashed line in all the plots indicates the value of ∆t_max obtained from maximizing SS_λ in Eq. (2) as a function of ∆t. For easy reference, the appropriate relaxation rate R_m values in units of Hz are shown in the small box in each plot. The results shown in the left-hand column (a–d) explore the effect of different R_m values for equal amplitude components while those shown in the right-hand column (e–h) have fixed R_m values with the respective {0.50,0.50}, {0.59,0.41}, {0.71,0.29}, and {0.91,0.09} relaxation coefficient amplitudes established using variable volume samples. The 〈N²〉 values used in the solid black line simulated results were obtained from the measured raw transient signal (color figure online)

A similar good agreement between experiment and theory is enjoyed in equal amplitude three and four relaxation coefficient studies as shown in the left- and right-hand columns in Fig. 4, respectively. Here three and four separate equal-volume containers, each loaded with one of the samples listed in Table 1, were simultaneously placed in the NMR detection coil and explored with the CPMG pulse sequence. The shape of the curves shown in Fig. 4(b, d) is to be expected by comparison to Fig. 3d. The difference between these results is the addition of the R₁ = 0.58 Hz and 2.23 Hz rate samples in Fig. 4(b, d), respectively. In addition to being more sensitive to signal noise as judged by the larger q∆t range of N_T > 0 values in Fig. 4(b, d), the added relaxation rate sample also shifts ∆t_max to longer values. The N_T values much closer to zero displayed in Fig. 4(a, c) are also to be expected with reference to Fig. 3(b, c), respectively. Here the only difference is the addition of a third R₁ = 0.58 Hz sample container to the NMR detection coil prior to generating the results shown in Fig. 4(a, c). Once again, the ∆t_max values shift to a longer time due to the addition of the third sample. Even though the N_T values for the three-component case suggest that the inherent experimental noise is too great to resolve the relaxation coefficients in Fig. 4d, the overall shape of the N_T curve is easily reconciled with respect to the two relaxation coefficient results shown in Fig. 3. Such a simple understanding does not seem possible when four relaxation coefficients are involved as shown in the right-hand column in Fig. 4. Like Fig. 4d, the N_T value in all the four component cases shown in Fig. 4e–h exceeds zero meaning that the MPM approach does not reliably resolve these four components given the level of experimental noise. Data sampling at the ∆t_max value still produces N_T values greater than zero although in most cases N_T is minimum at q∆t = ∆t_max.

Fig. 4

Comparison of experimental (open red circles) to theoretical (solid black line) N_T values for equal amplitude three (a–d) and four (e–h) component signals. The vertical dashed line in all the plots indicates the value of ∆t_max obtained from maximizing SS_λ in Eq. (2) as a function of ∆t. For easy reference, the appropriate relaxation rate R_m values in units of Hz are shown in the small box in each plot (color figure online)

The agreement between experiment and theory displayed in Figs. 3, 4 validates the reduced MPM approach used here to calculate N_T and its relationship to ∆t_max. This strong agreement allows one to ask how small the experimental noise 〈N²〉 must be to resolve all four relaxation coefficients for the sample quadruples used to develop the results on the right-hand side of Fig. 4. In all cases, it was found that all four components are resolved at ∆t_max if the noise level is decreased by a factor of 100. Since the signal to noise ratio is known to increase with the square of the number of scans [24], increasing the number of signals averaged from 8 per trial to 100^2 × 8 = 80,000 should recover acceptable resolution.

Another useful theoretical exercise considers the dependence of ∆t_max on the density and amplitude of the time constants within a fixed range of values as summarized in Fig. 5. The graph in (a) displays four separate time constant amplitude envelope functions. These functions include a box or equal amplitude distribution across the 0.1 s to 1.0 s time range, a Gaussian distribution and two oppositely skewed, bimodal Gaussian distributions. These amplitude distributions were used to generate the respective plots of ∆t_max as function of the density of relaxation coefficients in (b), or equivalently the number of relaxation coefficients # λ within the displayed fixed 0.1 s < T < 1 s range. The ∆t_max values were obtained by maximizing Eq. (2) as a function of ∆t and increasing number of relaxation coefficients. The ∆t_max value for two relaxation coefficients is identical in all cases as the two time constants T₁ = 0.1 s and T₂ = 1 s yield R₁ = 10 Hz and R₂ = 1 Hz and ∆t_max = log(R₂/R₁)/(R₂–R₁) = log(10 Hz/1 Hz)/(10 Hz–1 Hz) = 0.26 s. As the number of relaxation coefficients increases, the ∆t_max value steadily increases and asymptotically approaches a constant number that reflects the relaxation time distribution mean value. The similarity of the 0.35 s and 0.34 s asymptotic ∆t_max values shown for the solid and dashed curves in Fig. 5b is expected as the corresponding time constant amplitude distributions shown in Fig. 5a are symmetric and have the same mean time constant value. As the shape of the distribution becomes bimodal in Fig. 5a, the asymptotic ∆t_max value changes. The shift of average time constant to shorter T value for the dash-dotted curve in Fig. 5a causes the asymptotic ∆t_max value to shorten to 0.32 s in Fig. 5b. The opposite effect is observed for the dotted curve in Fig. 5b. Here the asymptotic ∆t_max value lengthens to 0.37 s as the average time constant shifts to longer T value in Fig. 5a. It should be clear that operation at ∆t_max does not guarantee that all the relaxation coefficients will be resolved. For example, it is not possible to resolve the ca. 1,000 relaxation coefficients used to reach the ∆t_max asymptotic limit. Rather, the asymptotic ∆t_max values shown in Fig. 5 suggest the best sampling rate needed to extract the maximum number of relaxation coefficients from the measured data provided enough signals have been averaged to keep the experimental noise level low and thus produce near zero N_T values.

Fig. 5

Dependence of ∆t_max on the number of uniformly spaced components in the 0.1 s < T < 1 s range. The asymptotic ∆t_max values shown in (b) for the uniform (solid), Gaussian (dashed), and oppositely skewed, bimodal Gaussian (dash-dotted and dotted) amplitude profiles in (a) are 0.35, 0.34, 0.32, and 0.37 respectively

5 Conclusion

The MPM is an algebraic way to analyze the damped, exponential, and transient signals common to all NMR relaxometry experiments. This work builds upon earlier results [18, 22] and demonstrates that the eigenvalues provided by MPM analysis z_m most faithfully reproduce the relaxation coefficients λ_m when the data sampling time ∆t_max is chosen such that the sum of square λ_m differences SS_λ is maximum as a function of data acquisition sampling time ∆t. The validity of ∆t_max in sample collections contrived to simultaneously produce two, three, and four relaxation coefficients is demonstrated by comparison of experiment to theory via the N_T parameter shown in Eq. (4). The analysis of two relaxation coefficient samples is attractive because a tractable analytical solution for N_T is obtained to first order in the mean square noise 〈N²〉. In the case of Gaussian white noise, higher nth-order analytical corrections in terms of 〈N²〉ⁿ could be developed. But judging from the agreement between the analytical and numerical results in the vicinity of ∆t_max, numerical methods were adopted as they also easily scale to situations presenting more than just two relaxation coefficients.

The obvious way to implement this work in the laboratory begins by obtaining a high signal–to–noise, oversampled transient relaxation decay signal. Oversampling guarantees that ∆t < ∆t_max. The number and the values of the relaxation time constants are estimated from this oversampled signal using conventional MPM, ILT, data fitting, etc. The m estimated time constants are used to construct the λ_m relaxation coefficients and SS_λ, the function that is maximized to yield ∆t_max. The measured transient data is then resampled at ∆t_max and MPM is used to find the inherent relaxation times and their relative contribution to the signal. If the actual number of time constants contributing to the signal is known, then the reduced MPM introduced here can be used in place of applying the MPM to the full resampled signal.