Magnetic resonance (MR) relaxometry and related experiments provide useful information about tissues, including brain, cartilage, and muscle^{1,2,3,4,5,6,7}. One-dimensional MR relaxometry experiments reveal the distribution of a parameter in a sample, such as the decay constant of transverse relaxation, $$$T_2$$$. The process required to recover the distribution, $$$m$$$, describing underlying macromolecular tissue components, from the observed data, $$$d$$$, is a Tikhonov-regularized discrete inverse Laplace transform (ILT) with non-negativity constraints. The ILT is highly ill-conditioned; low noise levels can result in large changes in the solution. Regularization reduces noise-sensitivity but decreases component resolution and distorts in component sizes. However, previous work has shown by simulations that substantially better stability is possible in 2D than in 1D, providing the potential for greatly improved characterization of tissue by relaxometry^8,9.

Here, new results are presented indicating the theoretical basis for improved stability in 2D. We present a new theorem bounding the condition number of the solution with respect to perturbations in data. For a series of simulated data consisting of two peaks at various positions, representing tissue components, we analyze stability in 1D and 2D by computing this condition number. These theoretical results are illustrated by extensive numerical experiments focusing on the ability to resolve two tissue components.

The problem to be solved is

$$$\min_{m\geq 0} ||Km-d||_2^2 + \alpha^2 ||m||_2^2. $$$ (1)

The problem can be solved efficiently using the Venkataramanan-Song-Hurlimann algorithm^10,11,12.

The data $$$d$$$ is of size $$$M \times 1$$$, the solution $$$m$$$ is of size $$$N\times 1$$$, the kernel $$$K$$$ is of size $$$M\times N$$$, and $$$\alpha$$$ controls regularization. The solution $$$m$$$ describes the distribution of relaxation components. In 2D the problem takes the same form, with $$$K:=K_1\otimes K_2$$$ and with $$$m$$$ and $$$d$$$ describing 2D data arrays flattened row-wise.

If $$$d$$$ is perturbed to $$$\tilde{d}$$$, the solution $$$m^{*}$$$ of (1) is perturbed to $$$\tilde{m}^{*}$$$. Without nonnegativity constraints, known results bound the perturbation in the solution⁹. We present a new theorem that extends the known condition number to (1).

Theorem: There exists a neighborhood $$$N(d)$$$ containing $$$d$$$ such that if $$$\tilde{d}\in N(d)$$$, the solution $$$\tilde{m}^{*}$$$ of (1) with $$$d$$$ replaced by $$$\tilde{d}$$$ satisfies $$\frac{||\tilde{m}^{*}-\tilde{m}^{*}||}{||m^{*}||} \leq\mathcal{K}\frac{||\tilde{d}-d||}{||d||},$$ where $$\mathcal{K}:=\frac{||d||}{||m^{*}||\sqrt{\sigma_{\min}(K)^2+\alpha^2}}\approx\frac{||d||}{||m^{*}||\alpha}=:\mathcal{K}_{\text{appx}}.$$

For various peak positions and regularization levels, we quantify $$$\mathcal{K}_{\text{appx}}$$$ and resolution, as measured by the height $$$\rho$$$ of the saddle point between peaks, with the height of the larger peak normalized to 1. Here $$$\rho=1$$$ indicates that two peaks are overlapping, and hence not well-resolved, while $$$\rho=0$$$ indicates that the peaks are perfectly resolved.

Define $$$\text{linspace}(a,b,n)$$$ to be $$$n$$$ points linearly-spaced from $$$a$$$ to $$$b$$$ and $$$\text{logspace}(a,b,n)$$$ to be $$$n$$$ points log-spaced from $$$10^a$$$ to $$$10^b$$$, as in Matlab. Define $$$T_2$$$ values $$$t[i]=\text{logspace}(-2.5,0.5,128)$$$ and acquisition times $$$\tau[i]=\text{logspace}(-3,1,128)$$$. The 2D kernel $$$K = K_1\otimes K_2$$$, where $$$K_l[i,j] = \exp(-\tau[i]/t[j])$$$, for $$$l = 1,2.$$$. Each $$$K_l$$$ is normalized so that $$$\sigma_{\max}(K_l)$$$ = 1.

For $$$\alpha =10^{i/2}$$$, $$$i=5,6\ldots,10$$$, data is simulated from initial distributions $$$m$$$ of size $$$128 \times 128$$$ with two nonzero entries, corresponding to peaks at positions $$$(0.1,0.1)$$$ and $$$(p_1,p_2)$$$. The relative peak positions $$$\lambda_i := \log_{10}(p_i / 0.1)$$$ are drawn from $$$\text{linspace}(-1.5,1.5,64).$$$ Gaussian noise is added at $$$\text{SNR}=1024$$$ and results are averaged over 4 realizations. Analogous simulations are performed in 1D for comparison.

The condition number from the simulations is substantially smaller in 2D than in 1D, indicating greatly improved stability. At $$$\alpha = 10^{-5}$$$, $$$\mathcal{K}_{\text{appx}}\sim5000$$$ for 2D and $$$\sim20000$$$ for 1D, indicating $$$4\times$$$ better stability. The contour plots of $$$\rho$$$ indicate peaks are well-resolved outside of a rounded-square shaped region. If two macromolecular components are close in 1D, introducing a second dimension in which the peaks are much more separated can greatly improve resolution and component quantification.We have presented a new theorem guaranteeing stability of the Tikhonov-regularized non-negative least squares problem arising in multidimensional MR relaxometry and demonstrated that the stability bound is substantially better in 2D than 1D. In extensive numerical simulations, we have measured stability and resolution, as quantified by the condition number $$$\mathcal{K}_{\text{appx}}$$$ and saddle height $$$\rho$$$, showing directly that the introduction of a second dimension can greatly improve signal component quantification and resolution if parameter separation is greater in the newly introduced dimension. These results provide a theoretical basis for the previously-described improvement in stability of the 2D ILT as compared to the 1D ILT. These results apply equally well to other 2D relaxometry and related experiments, including T1-T2 and ADC-T2 experiments. The stability of the 2D ILT may form the basis for a wide range of accurate experimental determinations of biophysical properties of tissue.