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

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.


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 1 3 Improving the Resolution of MPM Recovered Relaxometry… 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 twocomponent 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.

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 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) 2 〉 = 〈N 2 〉 ≠ 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 Ṽ 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 2 〉.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 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 using the noise tolerance Improving the Resolution of MPM Recovered Relaxometry… 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 .
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 λ 1 = exp(-R 1 ∆t) and λ 2 = exp(-R 2 ∆t) for the respective decay rates R 1 and R 2 .It should be clear that there are many choices of λ 1 and λ 2 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 λ 1 and λ 2 for fixed R 1 and R 2 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 λ 1 = λ 2 .In the absence of noise, λ 1 and λ 2 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 λ 1 = λ 2 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 λ 1 and λ 2 or for the ∆t value obtained from maximizing SS λ shown in Eq. ( 2) above.For two components, SS λ reduces to the difference λ 1 -λ 2 which is captured by the dotted line along the anti-diagonal in Fig. 1b.When maximized as a ( 4) 1 3 Improving the Resolution of MPM Recovered Relaxometry… function of ∆t, the difference yields the optimized time ∆t max = log(R 2 /R 1 )/(R 2 -R 1 ) and the most distinguishable relaxation coefficient . Graphically in the two-component case, the ∆t max value obtained from maximizing SS λ corresponds to finding the longest vector perpendicular to the λ 1 = λ 2 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 λ 1 max and λ 2 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.
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 in terms of the unit matrix E ∼ and the eigenvectors ⃗ p m is then solved for the eigenval- ues z m .Essentially one finds the inverse matrix and then the eigenvalues of the matrix product in the usual way as These analytical expressions for MPM analysis can be used to explore the effect of signal noise on the z ± eigenvalues and their deviation from λ 1 and λ 2 .For two relaxing components with amplitudes A 1 and A 2 , Eq. (1) implies that four consecutive signal points with noise beginning at t = 0 are ( 5) 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) 2 〉 = 〈N 2 〉 ≠ 0, results in two noise factors that can be used to approximate the two eigenvalues to first order in 〈N 2 〉 as 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 2 〉 and f _ 〈N 2 〉 factors appended to the λ 1,2 relaxation coefficients in the z +,-eigenvalues in Eq. (12).A larger average signal noise 〈N 2 〉 will generate a larger shaded gray region on this plot.Consideration of the eigenvalue difference, z + − z -= λ 1 − λ 2 + (f + − f -)〈N 2 〉, indicates that the λ 1 -λ 2 noise-free value increases to a larger 〈N 2 〉-dependent number where the increase calculated from Eq. ( 11), (f + − f -)〈N 2 〉, can be used to determine the noise tolerance defined in Eq. ( 4) as 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 ( 10) 1 3 Improving the Resolution of MPM Recovered Relaxometry… 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 λ 1 = exp(-R 1 q∆t), λ 2 = exp(-R 2 q∆t), and λ 3 = exp(-R 3 q∆t) relaxation coefficients for the three R 1 , R 2 , and R 3 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 λ 3 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 Y ∼ 1 and Y ∼ 2 matrices are constructed as and the 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 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.

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.
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 1 H Larmor frequency.All MPM analyses of experimental ( 14) 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 2 〉 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 2 ) could be used to accurately determine T 2 values (R 2 > 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 2 〉 1/2 so that the mean square average and cross-correlation at each ∆t reduce to 〈N((n-1)∆t) 2 〉 = 〈N 2 〉 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 .

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 2 〉.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.
The open circles shown in Fig. 2(a, b) describe the behavior of Eq. ( 13) as a function of q∆t and 〈N 2 〉 for two relaxing components with different R 1 and R 2 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 2 〉 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 2 〉 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 2 〉 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 2 〉 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 2 〉 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 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 2 〉 for three equal amplitude components with three different rates R 1,2,3 while Fig. 2d, like Fig. 2b, fixes 〈N 2 〉 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 2 〉 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 2 〉 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 2 〉 value.For example, for 〈N 2 〉 = 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 2 〉 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 1 = 2.23 Hz and R 2 = 10.2Hz 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 2 time constants.The value of R reported in Table 1 is the inverse of the T 2 value.These samples were prepared to 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.
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 righthand 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 1 = 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 1 = 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 .
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 2 〉 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 1 = 0.1 s and T 2 = 1 s yield R 1 = 10 Hz and R 2 = 1 Hz and ∆t max = log(R 2 /R 1 )/ (R 2 -R 1 ) = 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 have been averaged to keep the experimental noise level low and thus produce near zero N T values.

3
Improving the Resolution of MPM Recovered Relaxometry…

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 Fig. 5 Dependence of ∆t max on the number of uniformly 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 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 2 〉.In the case of Gaussian white noise, higher nth-order analytical corrections in terms of 〈N 2 〉 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 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.

Fig. 1
Fig.1The 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 2 〉 and f -〈N 2 〉 noise terms in the two relaxation coefficient eigenvalues shown in Eq.(12)

Fig. 2 3
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 2 〉 = 10 -10 (green), 10 -8 (black), 10 -6 (red), and 10 -4 (blue) noise levels used transient signals s(t) with equal amplitude A 1 = A 2 = 0.5 exponential decay functions.One 〈N 2 〉 = 10 -4 noise value was used to generate (b) for four different exponential decay function amplitude pairs {A 1 ,A 2 } = {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 2 〉 noise levels used transient signals s(t) with equal amplitude A 1 = A 2 = A 2 = 0.33 exponential decay functions.One 〈N 2 〉 = 10 -8 noise value was used to generate (d) for four different exponential decay function amplitude triples {A 1 ,A 2 ,A 3 } = {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)

Fig. 3
Fig. 3 Comparison of experimental (open red circles) to theoretical (solid black line) N T values for twocomponent 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 2 〉 values used in the solid black line simulated results were obtained from the measured raw transient signal (color figure online)

Fig. 4
Fig. 4 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)

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