<figure><img src="image_1.png" /></figure>

Published in final edited form as:

NMR Biomed. 2020 December; 33(12): e4366. doi:10.1002/nbm.4366.

# Parsimonious discretization for characterizing multi-exponential decay in magnetic resonance

Jean-Marie Bonny¹,², Amidou Traore¹,², Mustapha Bouhrara³, Richard G. Spencer³, Guilhem Pages¹,²

¹INRAE, UR QuaPA, Saint-Genès-Champanelle, France

²AgroResonance, INRAE, 2018, Nuclear Magnetic Resonance Facility for Agronomy, Food and Health, Saint-Genès-Champanelle, France

³NIH/National Institute on Aging - Baltimore, Maryland, USA

## Abstract

We address the problem of analyzing noise-corrupted magnetic resonance transverse decay signals as a superposition of underlying independently decaying mono-exponentials of positive amplitude. First, we indicate the manner in which this is an ill-conditioned inverse problem, rendering the analysis unstable with respect to noise. Second, we define an approach to this analysis, stabilized solely by the nonnegativity constraint without regularization. This is made possible by appropriate discretization, which is coarser than that often used in practice. Thirdly, we indicate further stabilization by inspecting the plateaus of cumulative distributions. We demonstrate our approach through analysis of simulated myelin water fraction measurements, and compare the accuracy with more conventional approaches. Finally, we apply our method to brain imaging data obtained from a human subject, showing that our approach leads to maps of the myelin water fraction which are much more stable with respect to increasing noise than those obtained with conventional approaches.

## Keywords

cumulative distribution function; inverse Laplace transform; myelin; simulation

# 1 | INTRODUCTION

## 1.1 | Formulation of the inverse problem

Free induction decays,¹⁻⁴ refocused echo trains⁵⁻⁷ and diffusion-weighted signals⁸⁻¹¹ are all examples of MR signals that represent a superposition of monoexponentially-decaying component signals. Multi-exponential analysis to define component amplitudes and decay time constants is then applied with the goal of identifying important structural or functional tissue characteristics.

<footnote>Correspondence: Jean-Marie Bonny, INRAE, UR QuaPA, F-63122 Saint-Genès-Champanelle, France. jean-marie.bonny@inrae.fr.
Richard G. Spencer and Guilhem Pages contributed equally to this article</footnote>
Without lack of generality, we take the observed NMR signal, $s(t)$, to be a decaying function of time. The most general model for this is a continuous distribution of exponentials

$$s(t) = \int_{\check{R}} a(R) \exp(-Rt) dR \quad t \in \check{t}, \qquad (1)$$

where $a$ is the desired distribution of decay rates $R$ within a specified domain $\check{R}$ of permitted values and $\check{t}$ is the domain of sampling times. Here and in the following, we are considering signal behavior within a voxel or region of interest, and ignoring spatial localization. The goal of the experiment is to determine $a(R)$ in order to provide a physical–chemical characterization of the sample. For example, determination of $a(T_2)$ can enable the quantification of macromolecular components by identifying relaxation time components defined by ranges of $T_2$ values, taken in this case to be $T_2 = 1/R$.

Equation 1 is a particular instance of a more general inverse problem¹² arising in many experimental settings¹³:

$$s(t) = (Ka)(t) = \int_{\check{R}} a(R) K(t, R) dR, \qquad (2)$$

known as the Fredholm equation of the first kind, with the kernel function $K$ in effect providing a transformation from the rate domain $\check{R}$ to the sampling time domain $\check{t}$. It is well known that, given $K$, the forward transformation of obtaining $s(t)$ from $a(R)$ is a smoothing operation, while the inversion of the Fredholm equation of the first kind, that is, the problem of determining the distribution $a(R)$ from the observed signal $s(t)$, amplifies noise and is highly unstable in the sense that even small amounts of variability in $s(t)$ can result in large variations in the derived $a(R)$.¹⁴

With $\check{t}$ and $\check{R}$ both restricted to $\mathbb{R}^+$, Equation 1 is a special case of the Laplace transform (LT) with the kernel function restricted to real-valued exponentials, describing a mapping between function spaces of real-valued functions $a \xrightarrow{LT} s$. The more general LT permits complex exponentials, resulting in a transformation from the domain of real-valued functions to the domain of complex-valued functions. This mirrors the manner in which the Fourier transform is a special case of the LT, with the former restricted to purely imaginary exponentials. These special cases can have very different properties from the general case. For the general LT, inversion may be carried out by complex integration through the so-called Bromwich integral; this approach is not available for real-valued functions.¹⁵

There have been an enormous number of MR papers addressing the problems of parameter estimation from multi-exponential decay (to illustrate, a PubMed search on 4 April 2020 for “biexponential” OR “triexponential” OR “multiexponential” and “NMR” OR “MRI” returned 904 hits). In spite of this, there are relatively few reviews comparing the available approaches for obtaining the desired distribution $a(R)$ from measured data $s(t)^{16–18}$; the choice of methods is largely a matter of computing resources and all methods are fundamentally limited by the ill-posed nature of the problem (see Istratov and Vyvenko¹⁶).
## 1.2 | Discretization of the problem

For calculations, a discretized form of Equation 1 is required as

$$s(t) = \sum_{i=1}^{m} a_i \exp(-R_i t) = \sum_{i=1}^{m} a_i \exp\left(-\frac{t}{T_{2,i}}\right) \quad t \in \tilde{t}. \quad (3)$$

The acquisition domain $\tilde{t}$ is defined between two values [$t_1, t_n$] based on experimental limitations and known or assumed sample characteristics related to the expected range of $T_2$ values. The minimum sampling time $t_1 > 0$ is determined by the duration of pulse sequence elements such as excitation and refocusing pulses, and spatial encoding gradients. In addition, sample considerations regarding the smallest values of $T_2$ that are of practical interest also serve to determine $t_1$. The maximum sampling time $t_n$ is similarly determined by sample characteristics, with the additional consideration that as the longest $T_2$ component approaches near-complete decay, the signal is dominated by noise.

The parameter estimation problem defined in Equation 3 is generally approached, broadly, in two somewhat different ways in the MR literature. The first is the case in which the number $m$ of exponential components is a priori known or assumed but both the relaxation times $T_{2,i}$ and amplitudes $f_i = \mathcal{I}(T_{2,i})$ must be estimated. The most widely used approach to this is through by minimizing a least squares difference criterion between measured and modeled values of $s(t)$ using various optimization algorithms$^{19,20}$; this is the nonlinear least squares (NLLS) approach. However, this is restricted to two, or at most three, exponential components due to the redundancy of exponential representations. Moreover, it requires initial guesses of the $2m$ unknowns that are sufficiently accurate to ensure convergence towards the global minimum. In fact, such initial guesses may be obtained through application of methods based on linear prediction$^{21}$ or Laplace-Padé approximation,$^{22}$ or through global optimization algorithms that perform one of several types of searches for optimal initial guesses. Note that the quality of least squares solutions depends not only on convergence but also on the properties of the noise, which should be centered, normally distributed, and uncorrelated with parameter values for optimal estimation. These conditions may not be met throughout the sampling domain $\tilde{t}$.

An alternative strategy for analysis of Equation 3 may be implemented when the number of decaying components is unknown, or is thought to be larger than two or three. This consists of discretizing the $\tilde{T}_2$ domain, with an exponentially decaying component potentially occurring at each value of $T_{2,i}$. This approach bypasses the need to explicitly define the number of components and, since the decaying signal is a linear combination of signals with predefined time constants, converts the problem to a linear formulation with the only unknowns being the amplitudes of the components associated with each of these decay components. The linear system may then be defined according to $m$ possible values of $T_2$ in the $\tilde{T}_2$ domain, ranging from $T_{2,1}$ to $T_{2,m}$, with data taken over $n$ time points in the sampling domain $\tilde{t}$. Then the linear formulation appears as

$$Ef = s, \quad (4)$$
where **f** is the column vector of *m* unknown amplitudes, **E** is a matrix of size *n* × *m* defined by $E_{ij} = \exp(-t_i/T_{2j})$ and **s** is the column vector composed of the *n* sampled signal values.

There may be approximately *m* = 100 candidate values of *T*_{2,i} and approximately *n* = 1000 sampling times, so that this linear system is necessarily large. It is well known that the condition number of a LT matrix **E** of this size is extremely large, so that simple inversion of Equation 4 via QR decomposition or the pseudoinverse will yield entirely unreliable results. This corresponds to the fact that the NLLS approach defined above can realistically be implemented only for two or three components.

In four seminal papers,¹²,²³⁻²⁵ Bertero et al considered the problem of LT inversion incorporating limits to the support of the unknown distribution. It was shown that greater resolution of exponential components is achievable through decreasing the ratio $\gamma = T_{2,m}/T_{2,1}$.¹² This motivates a strategy of first performing an inversion of Equation 4 over a wide range to define the support of $f$, and then using this restricted interval for a subsequent inversion.¹⁶ Even over this restricted range, the derived solutions are neither unique nor stable.²⁶ Nonuniqueness means that several solutions for $f$ are able to fit a noisy decay curve with virtually the same level of residual error (see figure 2 in Istratov and Vyvenko¹⁶), while instability means that derived solutions vary widely for different noise realizations; these problems persist even for high SNR signals.

Our goal in the present work is to describe an approach to the optimal interplay of nonnegativity constraints, regularization and discretization to improve the determination of the vector of unknown amplitudes in Equation 4.

1.3 | Strategies to reduce instability

A nonnegativity constraint consists in imposing the positivity on the amplitudes, as negative amplitudes are physically meaningless in the case under consideration. The imposition of this constraint greatly reduces the possible space of solutions, serving to improve the stability of Equation 4. The solution is often obtained through use of the iterative nonnegative least squares (NNLS) algorithm developed by Lawson and Hanson²⁷ or one of its variants²⁸,²⁹; this converges to a unique solution with positive amplitudes at each relaxation time $T_{2,i}$. A price to pay for imposing nonnegativity constraints is the requirement for use of an iterative algorithm for obtaining a solution³⁰ (see below for definition of the parameter $\delta$ which defines this resolution).

A second approach to stabilization is through regularization, such as one of the many variants of Tikhonov regularization.²⁸,³¹ For this, a specific penalty function, or constraint, must be defined based on prior knowledge or assumption regarding the properties of the desired solution; common choices are to penalize solution norm, total variation, or curvature. Having selected this, an optimal degree of regularization, defined by a parameter $\mu$ which scales the constraint, is calculated according to one of the several available methods.³²,³³ One popular strategy consists in iteratively updating $\mu$ to determine the value that results in a misfit of the data to the model by an a priori-defined factor relative to the unregularized fit.³⁴ Overall, the choice involved in a regularization strategy strongly affects the obtained
distributions, defined as the vector **f** in Equation 4, introducing a degree of subjectivity into the derived solution.

Lastly, discretization of $\tilde{T}_2$ has a large effect on stability and should ideally be performed in an optimal manner. Intuitively, fine discretization is desirable insofar as it more closely mimics the ideal continuous form of the problem, Equation 1. However, there are theoretical limits to the degree to which finer discretization can improve the resolution of the components comprising *a*; this limit is determined by the eigenvalue spectrum of the transformation $a \to s$ defined in Equation 2 and comes about physically through the loss of information induced by the integral operator in Equation 2 and corruption of the signal by noise. To illustrate the limits on discretization, let us consider a set of points $T_{2,i}$ in geometric progression$^{25,35}$ starting from $T_{2,1}$ and with $T_{2,i+1} = \delta \cdot T_{2,i}$, so that $T_{2,i+1} = \delta^i \cdot T_{2,i}$. With no further considerations, the maximum resolution available with such a discretization is defined by $\delta$, giving the proportional spacing between two consecutive points. In the unbounded case (ie, $\tilde{T}_2 = \mathbb{R}^+$), the minimum $\delta_i$ (ie, maximum resolution) for obtaining stable solutions depends only on SNR$^{12}$ and can be algebraically expressed as

$$ \delta_i = \exp\left(\frac{\pi^2}{\operatorname{acosh}(\pi \operatorname{SNR}^2)}\right). \qquad (5) $$

This resolution limit is shown in Figure 1, and indicates the slow increase in possible resolution with increasing SNR. Perhaps most importantly, it indicates a much larger $\delta_i$ value than is used for discretization in typical studies, indicating possible overoptimism on the part of investigators with respect to the value of finer discretization. However, there are at least two caveats to this assertion. First, even if two components of the underlying distribution are readily resolved, there may be value in a more precise knowledge of their centers. This is a somewhat different issue than resolution. Second, the limit defined in Equation 5 is for the worst case scenario, in which the support of the $T_2$ distribution is unknown so that the entire range of (0, infinity) must be considered. In practice, one has a great deal of prior knowledge regarding the plausible $T_2$ values represented in the sample, so that this value of $\delta_i$ can be decreased.

From these considerations, it may be reasonable to use a discretization finer than that defined by Equation 5, but coarser than that conventionally used in practice; we will designate this as a parsimonious (or moderate) discretization. This would serve to improve the condition number of the parameter extraction process as compared with finer discretization, while permitting more accurate peak localization than indicated by Equation 5. Ideally, such an approach, along with nonnegativity constraints, would obviate the need to incorporate regularization into the analysis.

Thus, as an alternative to the conventional regularized NNLS analysis, we study the accuracy and stability of moderately discretized NNLS without regularization (ie, the parsimonious discretization method [PDM]). We further extend this to evaluate cumulative, rather than point, densities, which may have certain advantages of stability when inspecting the plateau(s) of the cumulative distributions.
## 2 | MATERIALS AND METHODS

### 2.1 | Models and simulations

Continuum distributions were simulated as a weighted sum of Gaussians defined as:

$$f(T_2) = \sum_{i=1}^{n_G} \frac{\alpha_i}{\sqrt{2\pi}\sigma_i} \exp\left(-\frac{(T_2 - \mu_i)^2}{2\sigma_i^2}\right), \quad (6)$$

where $\alpha_i$, $\mu_i$ and $\sigma_i$ are normalized component amplitude ($\sum_{i=1}^{n_G} \alpha_i = 1$) and corresponding mean and standard deviation (SD) of the relaxation times within this component, respectively. This defines a continuous density function, for which $f(T_2)$ indicates the relative amplitude of the relaxation time component $T_2$ in the signal.

According to Equation 1, the time domain signal was then generated by

$$s(TE) = \int_{0}^{+\infty} f(T_2) \exp\left(-\frac{TE}{T_2}\right) dT_2 \quad (7)$$

and discretized over a linearly spaced set of TE values,²⁵,³⁵ $TE_j = i \cdot \Delta TE$. Gaussian-distributed noise of SD $\epsilon_n$ was added to these simulated decay curves to generate $p = 10 000$ noisy datasets for each defined SNR, defined as the ratio of signal intensity at TE = 0 divided by the noise SD. SNR was varied from 50 to 1000 in the simulations quantifying myelin water fraction (MWF). The decaying signal was sampled until its value decreased to $0.1 \times \epsilon_n$ for all simulations. This strategy provides a balance between the advantages of working with a large set of sampled data while limiting the influence of data that is primarily noise.

### 2.2 | NNLS inversion

NNLS was performed with the algorithm proposed by Lawson and Hanson²⁷ using MATLAB software (Natick, MA, USA). A geometric progression of $T_2$ values from a starting value of $T_{2,1}$ and defined by $T_{2,i+1} = \delta \cdot T_{2,i}$ was used in all cases. The ratio of the upper to lower limit of the defined support of $T_2$ was defined as $\gamma = \frac{T_{2,m}}{T_{2,1}}$. As shown by Bertero et al,¹² $\gamma$ is of central importance in the analysis, indicating in effect the degree of prior knowledge of the support of $T_2$. The degree of discretization was obtained by changing the number $m$ of $T_2$ values within a given support. Thus, each discretization is fully described by indicating $\gamma$ and $m$, and the resulting density $\delta = \frac{T_{2,i+1}}{T_{2,i}}$.

The result of NNLS analysis consists of a real-valued vector $\mathbf{x}$ of length $m$, with each element $i$ defining the value of the probability density function at the $i$-th $T_2$ value in the discretization, $f(T_{2,i})$. This vector minimizes the least-squares misfit $\|\mathbf{Ex} - \mathbf{S}\|^2$, where $\mathbf{x} \ge 0$, with the notation indicating that all elements of $\mathbf{x}$ are greater than or equal to zero, and $\mathbf{E}$ is the kernel matrix of size $n \times m$ described below (see Equation 4).
Tikhonov regularization was implemented by minimizing the function ‖Ex − S‖²+μ‖x‖², with the regularization parameter μ set through the discrepancy principle³⁴ as described above. μ was iteratively selected so that the chi-squared (χ²) value obtained from the regularized fit was between 102% and 102.5% of the χ² from the nonregularized fit.

Each discrete cumulative value cxᵢ at T₂ᵢ was calculated from all NNLS outputs according to

$$
cx_i = \frac{\sum_{j=1}^{i} x_j}{\sum_{j=1}^{m} x_j}. \qquad (8)
$$

Datasets and corresponding T₂ probability and cumulative probability distributions were obtained over 10 000 noise realizations of the underlying model with NNLS analysis.

2.3 | Single Gaussian simulations

A single Gaussian distribution was first simulated for comparing cumulative and density representations without incorporation of a regularization term. The discretization was fixed according to T<sub>2,1</sub> = 1 ms and T<sub>2,m</sub> = 50 ms (ie, γ = 50), with m = 30 (δ = 1.14). The decaying signal was sampled with ΔTE = 1 ms.

2.4 | Simulation analysis of MWF

Simulations mimicking the relaxation of myelin water offer a special case of the multi-Gaussian model given in Equation 6. In this framework, the T₂ distribution is obtained, with all values below a fixed value assigned to MWF and all greater values assigned to intra/extracellular water (IEW); we selected a cut-off value of 40 ms.³⁶

Simulation parameters are based on measurements performed on the genu of the corpus callosum in human subjects,³⁹ so that n<sub>G</sub> = 2. The component fraction α<sub>1</sub> of the myelin pool was assigned values of α<sub>1</sub> = 5%, 10% or 15%, centered on μ<sub>1</sub> = 15 ms with SD σ<sub>1</sub> = 3 ms, while IEW was assigned a relative fraction of α<sub>2</sub> = 1 − α<sub>1</sub>, centered on μ<sub>2</sub> = 70 ms with SD σ<sub>2</sub> = 5 ms. The decaying signal was sampled with ΔTE = 10 ms across a range of SNR to reflect a conventional acquisition scheme. The values T<sub>2,1</sub> = 10 ms and T<sub>2,m</sub> = 320 ms were fixed, encompassing the expected range of relaxation times (ie, γ = 32), with m varying and defining the density of discretization. We used linear interpolation to define F(T<sub>2</sub>) between the values of T<sub>2</sub> spanning the MWF cut-off of 40 ms.

MWF simulations were first performed for studying discretization only. Hence, the objective function for minimization did not include a regularization term. Next, the utility of additional regularization was evaluated by comparing two approaches which differed in terms of discretization and regularization: a moderately discretized NNLS (m = 10, γ = 32) without regularization (PDM) and a densely discretized NNLS (m = 80, γ = 32) including L2 regularization.
## 2.5 | Quantification of MWF by quantitative $T_2$ MRI/real data

These approaches were also tested on $T_2$ relaxation data obtained on a cognitively unimpaired 24-year-old male volunteer. Three-dimensional gradient and spin echo (GRASE) images were acquired with $n = 32$ echoes with $TE_i = i \cdot \Delta TE$ with $\Delta TE = 11.3$ ms, repetition time = 1000 ms, echo planar imaging factor = 3, field of view = 278 mm × 200 mm × 30 mm, acquisition matrix size = 185 × 133 × 10, acquisition voxel size = 1.5 mm × 1.5 mm × 3 mm, reconstructed to 1 mm × 1 mm × 3 mm using zero filling in k-space. The total acquisition time was 10 minutes. All experiments were performed on a 3 T whole body Philips MRI system (Achieva, Best, the Netherlands) using the internal quadrature body coil for transmission and an eight-channel phased-array head coil for reception.

The 32 multi-echo images were then denoised using the nonlocal estimation of multispectral magnitudes (NESMA) filter, which has been shown to improve MWF estimates from such GRASE imaging data.⁴⁰ Gaussian noise was added to these filtered images and was the dominant determinant of SNR, allowing us to vary SNR across a defined range. The SD $\epsilon_n$ of the added noise ranged from 0.1 to 10, with, for example, $\epsilon_n = 10$ leading to an overall SNR of ~110.

In line with the previous simulations, $p = 100$ MWF maps were first reconstructed in WM voxels by taking the value for cumulative distribution at $T_2 = 40$ ms. Two inversion methods were investigated with the decomposition basis being fixed to $T_{2,1} = 9$ ms and $T_{2,m} = 300$ ms ($\gamma = 33.3$): $m = 10$ without regularization (PDM) and $m = 80$ with $L_2$ regularization. Mean MWF maps were then computed on a pixel-by-pixel basis, and compared with the ones obtained without added noise.

To study the effect of noise on the MWF estimation in more detail, inversions were also performed on several pixels from different parts of the WM, with the analysis as described above performed with a larger $p = 10 000$.

# 3 | RESULTS

Figure 2 illustrates the markedly reduced variability of $cx_i$, a cumulative sum which is an estimation of $F$, as compared with $x_i$ in describing $f(T_2)$. The example shown is for a single Gaussian pool with $\mu_1 = 10$ ms and $\sigma_1 = 3$ ms. Similar results were obtained for a wide range of conditions exhibiting varying degrees of instability (data not shown).

Furthermore, these simulations permit inspection of how well the discrete cumulative distribution obtained by NNLS without regularization (ie, computed by Equation 8) approximates a continuous cumulative distribution function. This is especially important in the case of a coarse discretization, to confirm adequate sampling. We observed a good correspondence between the estimated discrete distributions and the underlying continuous ones, taking into account the fact that the $T_2$ axes cannot coincide perfectly simply due to the mechanics of discretization.

Simulation analyses allowed us to heuristically determine that $cx_i \approx F(\dot{T}_{2,i})$, where $\dot{T}_{2,i} = \frac{T_{2,i} + T_{2,i+1}}{2}$. This result also serves to demonstrate how the NNLS discrete outputs
(shown in Figure 2B) are related to the continuous density function *f* (see Appendix). The SD plots show that both cumulative and density discrete outputs exhibit their greatest variability near the peak of the density curve *f*. In other words, the NNLS amplitude estimates are very noisy due to instability, and therefore difficult to use on a pixel-by-pixel basis. In contrast, the calculated cdf values exhibit greater stability. On the other hand, the skewness assesses the extent to which the distribution of these cumulative values is symmetrical from either side of the mean. Here, there is only relatively weak asymmetry at the plateau as illustrated by the skewness calculated from the *p* cumulative values obtained at each *T*<sub>2,i</sub> by simulation. This figure also shows that the density values are seen to be distributed more asymmetrically compared with the cumulative ones.

Figure 3 shows bias and uncertainty in MWF estimation as a function of SNR, performed as above; note that these results involve only discretization with different values of *m* without regularization. The left-hand panel (A) shows that for each level of discretization, bias is fairly constant across a wide range of SNR, but increases substantially with increasing *m*. Note in the right-hand panel (B) that the finer discretization, as expected, exhibit greater SD than the coarser and parsimonious discretizations. Because the variance of derived solutions is markedly worse at low than at high SNR, the nonnegativity constraint is clearly insufficient. These results illustrate (i) why an additional regularization was introduced for reducing bias when fine discretization is used, and (ii) the need for high SNR data when attempting to reduce or eliminate the need for regularization.

Results incorporating Tikhonov regularization implemented as described above are shown in Figure 4. Here, the goal is to compare our proposed scheme PDM with the more conventional approach, which involves fine discretization with regularization. Comparison of the red curve in the right-hand panel (B) with the magenta curve in Figure 3B shows that instability is reduced by regularization, as expected. Of greatest interest is the observation of the left-hand panel (A) that bias is lower with a parsimonious basis (*m* = 10) without regularization, and across a wider SNR range, than with conventional discretization and regularization.

Further simulations of the type shown in Figure 4 were performed for a range of MWF values (Figure 5). The results for all three values of MWF investigated were qualitatively similar to those of Figure 4, with bias being lower with a parsimonious basis (*m* = 10) without regularization, and across a wider SNR range, than with conventional discretization and regularization. One point of difference is for the regime of very low SNR for the lowest MWF; this produces highly inaccurate parameter estimates. Note that the SD is virtually the same for both of the low- and high-resolution cases and for all three of the MWF values, since the red curve has finer discretization but is regularized, while the blue one has coarser regularization without regularization (PDM).

Figure 6 shows MWF maps obtained on a human subject. Panel (A) shows PDM, without regularization, while panel (B) shows the more conventional dense discretization with regularization. These results are very similar, showing that both methods are successful at high SNR. Panel (C), for comparison with panel (A), shows the results of MWF determination using PMD, averaged over 100 noise realizations with SNR severely
decreased by the addition of noise of SD $\epsilon_n = 1$. The results are virtually indistinguishable from panel (A), indicating the lack of bias in this approach. Panel (D), for comparison with panel (B), shows the results of MWF determination using the conventional fine discretization with regularization analysis, averaged over noise realizations in the same way as described for panel (C). Clear bias is seen in comparisons with panel (B), particularly in the frontal white matter.

Mean and variance can be assessed more quantitatively by investigating results over a larger ensemble ($p = 10\,000$) of noise realizations (Figure 7). This figure shows the mean and SD of MWF estimates for six representative pixels in the brain image. The mean of derived MWF values as a function of the SD of added noise is shown, starting from the best available SNR, and variance is shown as the SD of MWF estimates across that same range. While comparable in many voxels, MWF estimates obtained by the two methods with added noise of SD of 0.1 also differ in several voxels (eg, voxels e, d, f in the occipital white matter), indicating bias. In addition, the SD of MWF values resulting from addition of noise was generally greater for the conventional analysis. Thus, parsimonious discretization displays superior stability properties compared with the conventional analysis.

## 4 | DISCUSSION

Multi-exponential inversion based on the NNLS algorithm is widely used in MR relaxometry and related experiments. As expressed in Equation 3, the desired component coefficients $f_i = \int T_{2,i} d\mathbf{r}$ can be defined on only a finite set of points,²⁵ representing an approximation of the continuum of potential $T_2$ components. Therefore, intuitively, it would appear that increased discretization in $T_2$ would improve the accuracy of the inversion, and permit a more faithful representation of this continuum. However, in contrast to other settings in numerical analysis, finer discretization in fact produces a problem more similar to the underlying ill-posed Fredholm equation of the first kind, resulting in greater instability in the inversion. Although this effect has been discussed extensively in the mathematical literature on inverse problems, it does not seem to be well recognized in the MR literature. Rather, fine discretization is performed, potentially beyond the useful limit as defined in Equation 5 with the resulting instability managed through use of regularization schemes.

Here we propose an alternative procedure, in which we implement a coarser, that is parsimonious discretization to eliminate the need for regularization. This avoids the difficult and, to some extent, subjective task of selecting and optimizing a regularization approach. Indeed, we find that across a large range of SNR values, this regularization-free approach can lead to higher quality $T_2$ distributions (Figures 4 and 5). These conclusions were confirmed using human brain imaging data. All of these results are in line with the seminal works of Bertero et al,¹² who remarked that the problem with continuous data can be very well approximated by a problem with a small number of geometrically spaced $T_{2,i}$.

Having decided upon a parsimonious discretization, there is still a choice to be made regarding $T_{2,1}$ and $T_{2,m}$. These can generally be chosen based on physical a priori knowledge of the system under study, with additional consideration of experimental limitations in recovering signal from very short $T_2$ species. However, an appropriate
choice of $T_{2,1}$ is dependent not only upon the shortest echo time available, but also on SNR considerations. Studies of the selection of initial sampling times and in fact optimal sampling times throughout the signal decay have appeared in the literature,⁴¹⁻⁴³ but further discussion of this important topic is beyond the scope of the present work. We emphasize that while the limits defined by $T_{2,1}$ and $T_{2,m}$ should be generous enough to capture the range of $T_2$ values in the sample, in accord with previous work,¹² the ratio $\gamma = \frac{T_{2,m}}{T_{2,1}}$ should be kept as small as possible to maintain maximum accuracy.

Figures 4, 5 and 7 highlight that the stability of this inverse problem is highly SNR-dependent as uncertainty generally decreases faster than linearly as a function of SNR. This indicates the need for optimal acquisition conditions as well as careful implementation of noise reduction in postprocessing. There are a number of effective filters for this purpose, including nonlinear multispectral filtering, which has been shown to improve MWF estimation without losing fine anatomical detail⁴⁰,⁴⁴

In this work, we have emphasized the stabilizing effect of using a moderate discretization in the $T_2$ reconstruction domain. The density of discretization as defined by $m$ is directly linked with the limit defined by Bertero et al (equation 5¹²) and the value conventionally used with the densely discretized approach including regularization. The ultimate choice should be driven by a compromise between achieving sufficient resolution to adequately define the desired $T_2$ distribution, and the desire to avoid regularization, or to at least minimize the bias associated with it.

Finally, we note that in practice, a NNLS inversion of decay data is sometimes first performed to define the number of components in a sample and to select reasonable starting values for subsequent NLLS analysis. To the extent that the PDM improves stability and decreases error in NNLS results, it stands to reason that it would improve the accuracy of this sequential procedure.

## 5 | CONCLUSION

The conventional use of a dense decomposition basis in NNLS analysis has resulted in nearly universal application of some type of regularization for analysis of multi-exponential decay in MR. However, we show through simulations of MWF estimation by $T_2$ relaxometry that nonnegativity constraints alone, combined with parsimonious discretization of the reconstruction domain, may be more effective over a large range of SNR values than use of a dense discretization with regularization. Additional investigations are warranted to evaluate this parsimonious discretization approach for other distributions and experimental settings.

## ACKNOWLEDGEMENTS

This work was supported in part by an AgroResonance intramural funding (JMB, AT and GP) and by the Intramural Research Program of the National Institute on Aging of the NIH (RGS and MB).

## APPENDIX A

By definition:
$$ \frac{dF(T_2)}{dT_2} = f(T_2). \qquad (A1) $$

Experimentally, it has been observed that $cx_i \approx F(\dot{T}_{2,i})$. We can then derive the cdf using finite differences:

$$ \begin{aligned} \frac{cx_{i+1} - cx_i}{\dot{T}_{2,i+1} - \dot{T}_{2,i}} &= \frac{2x_{i+1}}{T_{2,i+2} - T_{2,i}} \approx f\left(\frac{\dot{T}_{2,i+1} + \dot{T}_{2,i}}{2}\right) \\ &= f\left(\frac{T_{2,i} + 2T_{2,i+1} + T_{2,i+2}}{4}\right). \end{aligned} \qquad (A2) $$

This expression indicates that the discrete value $\frac{2x_i + 1}{T_{2,i+2} - T_{2,i}}$ is an estimate of the continuous density $f$ for $T_2 = \frac{T_{2,i} + 2T_{2,i+1} + T_{2,i+2}}{4}$, and was used to represent in concert the point-wise density function $f$ derived from NNLS outputs and the underlying continuous density function (Figure 2).

Abbreviations used:

<table><tr><td>cdf</td><td>cumulative density function</td></tr><tr><td>df</td><td>density function</td></tr><tr><td>ILT</td><td>inverse Laplace transform</td></tr><tr><td>LT</td><td>Laplace transform</td></tr><tr><td>MRI</td><td>magnetic resonance imaging</td></tr><tr><td>MWF</td><td>myelin water fraction</td></tr><tr><td>NLLS</td><td>nonlinear least squares</td></tr><tr><td>NMR</td><td>nuclear magnetic resonance</td></tr><tr><td>NNLS</td><td>nonnegative least squares</td></tr><tr><td>PDM</td><td>parsimonious discretization method</td></tr><tr><td>SNR</td><td>signal-to-noise ratio</td></tr><tr><td>TE</td><td>echo time</td></tr></table>

REFERENCES

1. Yablonskiy DA, Haacke EM. Theory of NMR signal behavior in magnetically inhomogeneous tissues: the static dephasing regime. Magn Reson Med. 1994;32(6):749–763. [PubMed: 7869897]
2. Yablonskiy DA, Sukstanskii AL, He X. Blood oxygenation level-dependent (BOLD)-based techniques for the quantification of brain hemodynamic and metabolic properties - theoretical models and experimental approaches. NMR Biomed. 2013;26(8):963-986. [PubMed: 22927123]

3. Sati P, van Gelderen P, Silva AC, et al. Micro-compartment specific T2* relaxation in the brain. Neuroimage. 2013;77:268-278. [PubMed: 23528924]

4. Hwang D, Kim DH, Du YP. In vivo multi-slice mapping of myelin water content using T2* decay. Neuroimage. 2010;52(1):198-204. [PubMed: 20398770]

5. MacKay A, Whittall K, Adler J, Li D, Paty D, Graeb D. In vivo visualization of myelin water in brain by magnetic resonance. Magn Reson Med. 1994;31 (6):673-677. [PubMed: 8057820]

6. Does MD, Snyder RE. T2 relaxation of peripheral nerve measured in vivo. Magn Reson Imaging. 1995;13(4):575-580. [PubMed: 7674853]

7. Drenthen GS, Backes WH, Aldenkamp AP, Op't Veld GJ, Jansen JFA. A new analysis approach for T2 relaxometry myelin water quantification: Orthogonal Matching Pursuit. Magn Reson Med. 2019;81(5):3292-3303. [PubMed: 30444019]

8. Milne ML, Conradi MS. Multi-exponential signal decay from diffusion in a single compartment. J Magn Reson. 2009;197(1):87-90. [PubMed: 19121965]

9. Niendorf T, Dijkhuizen RM, Norris DG, van Lookeren Campagne M, Nicolay K. Biexponential diffusion attenuation in various states of brain tissue: implications for diffusion-weighted imaging. Magn Reson Med. 1996;36(6):847-857. [PubMed: 8946350]

10. Maier SE, Vajapeyam S, Mamata H, Westin CF, Jolesz FA, Mulkern RV. Biexponential diffusion tensor analysis of human brain diffusion data. Magn Reson Med. 2004;51(2):321-330. [PubMed: 14755658]

11. Marchand AJ, Hitti E, Monge F, et al. MRI quantification of diffusion and perfusion in bone marrow by intravoxel incoherent motion (IVIM) and nonnegative least square (NNLS) analysis. Magn Reson Imaging. 2014;32(9):1091-1096. [PubMed: 25093628]

12. Bertero M, Boccacci P, Pike ER. On the recovery and resolution of exponential relaxation rates from experimental data: a singular-value analysis of the Laplace transform inversion in the presence of noise. Proc R Soc A. 1982;383:15-29.

13. Shtrauss V Decomposition filters for multi-exponential and related signals. Int J Math Mod Meth Appl Sci. 2007;1(3):137-142.

14. Wing GM, Zahrt JD. A Primer on Integral Equations of the First Kind: The Problem of Deconvolution and Unfolding. Philadelphia, PA: Society for Industrial and Applied Mathematics; 1991.

15. Fordham EJ, Venkataramanan L, Mitchell J, Valori A. What are, and what are not, inverse Laplace transforms. Diff Fundamentals. 2017;29(2):1-8. https://diffusion.uni-leipzig.de/pdf/volume29/ diff_fund_29(2017)2.pdf

16. Istratov AA, Vyvenko OF. Exponential analysis in physical phenomena. Rev Sci Instrum. 1999;70(2):1233-1257.

17. Whittall KP. Inversion of magnetic resonance data: How to T2. 2006. http://www3.telus.net/ InversionsInc/t2.html

18. Borgia CG, Bortolotti V, Brown RJS, Castaldi P, Fantazzini P, Soverini U. A comparison among different inversion methods for multi-exponential NMR relaxation data. Magn Reson Imaging. 1994;12(2):209-212. [PubMed: 8170301]

19. Bjorck A Numerical methods for least squares problems. 51 Philadelphia, PA: Siam; 1996.

20. Hansen PC, Pereyra V, Scherer G. Least squares data fitting with applications. Baltimore, MD: JHU Press; 2013.

21. Lin YY, Ge NH, Hwang LP. Multiexponential analysis of relaxation decays based on linear prediction and singular-value decomposition. J Magn Reson. 1993;105(1):65-71.

22. Yeramian E, Claverie P. Analysis of multiexponential functions without a hypothesis as to the number of components. Nature. 1987;326(6109): 169-174.

23. Bertero M, Brianzi P, Parker P, Pike ER. Resolution in diffraction-limited imaging, a singular value analysis. III The effect of sampling. Opti Acta Int J Opt. 1984;31(2):181-201.
24. Bertero M, Brianzi P, Pike ER. On the recovery and resolution of exponential relaxational rates from experimental data: Laplace transform inversions in weighted spaces. Inverse Problems. 1985;1(1):1–15.

25. Bertero M, Boccacci P, Pike ER. On the recovery and resolution of exponential relaxation rates from experimental data II. The optimum choice of experimental sampling points for Laplace transform inversion. Proc Roy Soc Lond A-Math Phy Sci. 1984;393:51–65.

26. Kabanikhin SI. Definitions and examples of inverse and ill-posed problems. J Inverse Ill-Pose P. 2008;16:317–357.

27. Lawson CL, Hanson RJ. Solving least squares problems. Philadelphia, PA: Society for Industrial and Applied Mathematics; 1995.

28. Whittall KP, MacKay AL. Quantitative interpretation of NMR relaxation data. J Magn Reson. 1989;84:134–152.

29. Bro R, De Jong S. A fast non-negativity-constrained least squares algorithm. J Chemometr. 1997;11(5):393–401.

30. Bertero M, Brianzi P, Pike ER, Rebolia L. Linear regularizing algorithms for positive solutions of linear inverse problems. Proc Roy Soc Lond A-Math Phy Sci. 1988;415:257–275.

31. Bjarnason TA, Mitchell JR. Analyze NNLS: magnetic resonance multiexponential decay image analysis. J Magn Reson. 2010;206(2):200–204. [PubMed: 20688549]

32. Zou Y, Xie RH, Arad A. Numerical estimation of choice of the regularization parameter for NMR T2 inversion. Petr Science. 2016;13(2):237–246.

33. Sabett C, Haffftka A, Sexton K, Spencer RG. L1, Lp, L2, and elastic net penalties for regularization of Gaussian component distributions in magnetic resonance relaxometry. Concept Magn Reson A. 2017;46A(2):1–20. e21427.

34. Graham SJ, Stanchev PL, Bronskil MJ. Criteria for analysis of multicomponent tissue T2 relaxation data. Magn Reson Med. 1996;35(3):370–378. [PubMed: 8699949]

35. Ostrowsky N, Sornette D, Parker P, Pike ER. Exponential sampling method for light scattering polydispersity analysis. Opt Acta. 1981;28(8): 1059–1070.

36. Kolind SH, Mädler B, Fischer S, Li DKB, MacKay AL. Myelin water imaging: Implementation and development at 3.0T and comparison to 1.5T measurements. Magn Reson Med. 2009;62(1):106–115. [PubMed: 19353659]

37. Dortch RD, Yankeelev TE, Yue Z, Quarles CC, Gore JC, Does MD. Evidence of multiexponential T2 in rat glioblastoma. NMR Biomed. 2009;22(6): 609–618. [PubMed: 19267385]

38. Reiter DA, Irrechukwu O, Lin PC, et al. Improved MR-based characterization of engineered cartilage using multiexponential T2 relaxation and multivariate analysis. NMR Biomed. 2012;25(3):476–488. [PubMed: 22287335]

39. Whittall KP, MacKay AL, Graeb DA, Nugent RA, Li DK, Paty DW. In vivo measurement of T2 distributions and water contents in normal human brain. Magn Reson Med. 1997;37(1):34–43. [PubMed: 8978630]

40. Bouhrara M, Reiter DA, Maring MC, Bonny JM, Spencer RG. Use of the NESMA filter to improve myelin water fraction mapping with brain MRI. J Neuroimaging. 2018;28(6):640–649. [PubMed: 29999204]

41. Shrager RI, Weiss GH, Spencer RG. Optimal time spacings for T2 measurements: monoexponential and biexponential systems. NMR Biomed. 1998;11 (6):297–305. [PubMed: 9802472]

42. Moody JB, Xia Y. Analysis of multi-exponential relaxation data with very short components using linear regularization. J Magn Reson. 2004;167(1): 36–41. [PubMed: 14987596]

43. Jones JA, Hodgkinson P, Barker AL, Hore PJ. Optimal sampling strategies for the measurement of spin-spin relaxation times. J Magn Reson. 1996;113 (1):25–34.

44. Bouhrara M, Bonny JM, Ashinsky BG, Maring MC, Spencer RG. Noise estimation and reduction in magnetic resonance imaging using a new multispectral nonlocal maximum-likelihood Filter. IEEE Trans Med Imaging. 2017;36(1):181–193. [PubMed: 27552743]
<figure><img src="image_2.png" /><figcaption>FIGURE 1.</figcaption></figure>

Theoretical resolution limit $\delta$ as a function of SNR for an infinite support, $T_2 = \mathbb{R}^+$. Black, purple and red lines show the resolutions used in references,³⁶⁻³⁸ respectively
<figure><img src="image_3.png" /><figcaption>FIGURE 2.</figcaption></figure>

Comparison of (A) cumulative and (B) density representations obtained by averaging 10 000 solutions provided by NNLS without regularization (red crosses). Blue lines correspond to the theoretical single Gaussian distribution (μ₁ = 10, σ₁ = 3). The inversions were performed with SNR = 1000 and midrange resolution γ = 50, m = 30 (δ = 1.14).
<figure><img src="image_4.png" /><figcaption>FIGURE 3.</figcaption></figure>

Illustration of the influence of the decomposition basis density on the accuracy of the estimates as a function of SNR (without regularization); 10 000 myelin-like model (MWF = 10%) simulations were performed and the averaged (A) bias (ie, the difference between the estimated and input values) and (B) SD for (blue) low resolution $\gamma = 32$, $m = 10$ ($\delta = 1.47$), (green) midrange resolution $\gamma = 32$, $m = 30$ ($\delta = 1.13$) and (magenta) high resolution $\gamma = 32$, $m = 80$ ($\delta = 1.04$) are plotted. Hence, the blue curve corresponds to PDM
<figure><img src="image_5.png" /><figcaption>FIGURE 4.</figcaption></figure>

Comparison of the (blue) low resolution PDM, with $\gamma = 32$ and $m = 10$ ($\delta = 1.47$), with the (red) dense discretization approach with regularization, with $\gamma = 32$ and $m = 80$ ($\delta = 1.04$). The average (A) bias and (B) SD over 10 000 simulations of the myelin-like model are plotted as a function of the data SNR
<figure><img src="image_6.png" /><figcaption>(B)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(C)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(D)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(E)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(F)</figcaption></figure>

<figure><img src="image_11.png" /></figure>

FIGURE 5.

Sensitivity of MWF estimates. The input MWF was set to (top) 5%, (middle) 10% and (bottom) 15%. The averaged bias (A, C, E) and SD (B, D, F) of the 10 000 simulations of the myelin-like model are plotted as a function of the data SNR. Blue: PDM, ie, low resolution discretization with $\gamma$ = 32 and $m$ = 10 ($\delta$ = 1.47); red: dense discretization with regularization, with $\gamma$ = 32 and $m$ = 80 ($\delta$ = 1.04)
<figure><img src="image_12.png" /><figcaption>(B)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>MWF (%)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(C)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(D)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>MWF (%)</figcaption></figure>

<figure><img src="image_17.png" /></figure>

FIGURE 6.

White matter MWF mapping of a human subject using NNLS analysis with (A, C) PDM, ie, low resolution $\gamma$ = 33.3 and $m$ = 10 ($\delta$ = 1.47), and (B, D) the dense discretization approach with regularization, ie, high resolution $\gamma$ = 33.3 and $m$ = 80 ($\delta$ = 1.04). For (A) and (B), no noise was added to the original high SNR (~790) filtered images, while for (C) and (D), additional noise ($\epsilon_n$ = 10) was added and the mean MWF was estimated over $p$ = 100 noise realizations
<figure><img src="image_18.png" /><figcaption>FIGURE 7.</figcaption></figure>

Analysis of individual voxels of the human brain imaging dataset. The average and SD of MWF estimates over 10 000 noise realizations are plotted as a function of the SD of added noise. MWF was calculated as described in the text, using a cut-off value of 40 ms. Voxels were chosen either in the (A-C) frontal or (D-F) occipital region within the (A, E) left, (C, D) middle and (B, F) right hemispheres. Analysis was performed with (blue) PDM, ie, low resolution discretization with $\gamma$ = 33.3 and $m$ = 10 ($\delta$ = 1.47), and (red) dense discretization with regularization, with $\gamma$ = 33.3 and $m$ = 80 ($\delta$ = 1.04).