Synopsis

Tikhonov regularization and related methods are widely used in recovering relaxation time distributions in magnetic resonance relaxometry. Regularization optimization methods such as the L-curve and generalized cross-validation (GCV) identify a single optimized solution as the best approximation to the underlying distribution. In contrast, we propose a new reconstruction method, Multi-Reg, incorporating a range of regularized solutions. Multi-Reg is based on a dictionary of noise-corrupted regularized reconstructions of distribution basis functions. We demonstrate that Multi-Reg can out-perform L-curve or GCV analyses in simulation analyses of Gaussian distribution components, and with experimental results on mouse spinal cord and human muscle.

Purpose

Multicomponent MR relaxometry defines a distribution of relaxation times from signal amplitudes as a function of echo time. However, this ill-posed inverse problem ¹ exhibits solutions highly unstable to noise. Stabilization is achieved through methods such as Tikhonov regularization ¹. Letting $$$A$$$ be the design matrix, $$$\mathbf{f}$$$ the unknown desired distribution, $$$\mathbf{d}_{\text{obs}}$$$ the acquired data, and $$$\lambda$$$ the regularization parameter, the relaxometry inverse problem is:

$$\mathbf{f}^*=\arg\min_{\mathbf{f}\succeq \mathbf{0}}\left\|\mathbf{d}_{\text{obs}}-A\mathbf{f}\right \|_2^2+\lambda\left\|\mathbf{f}\right\|_2^2\tag{1}$$

Optimal $$$\lambda$$$ is selected through e.g. the L-curve ² or generalized cross-validation (GCV) ³. However, regularization effects differ according to the underlying distribution to be recovered. Thus, additional information is available through incorporating a range of $$$\lambda$$$. We introduce what we believe to be the first analysis of inverse problems incorporating this concept, termed Multi-Reg.

Theory

We define a sequence of distributions, $$$\left\{\mathbf{f}_j\right\}_1^N$$$, obtained from regularization across a range of $$$\left\{\lambda_j\right\}_1^N$$$ and express the unknown distribution as a linear combination:

$$\mathbf{f}_{\text{true}}\approx\sum_{j =1}^N\alpha_j\mathbf{f}_j\tag{2}$$

where $$$\sum_{j =1}^N\alpha_j=1$$$. We also define a Gaussian basis $$$\left\{\mathbf{g}_i\right\}_1^M $$$ for the reconstructed distribution:

$$\mathbf{f}_{\text{true}}\approx\sum_{i =1}^Mc_i\mathbf{g}_i\tag{3}$$

and perturb the time-domain representation of these basis elements by adding the same level of noise as in the measured signal, $$$\mathbf{d}^{(i)}_{\text{obs}}$$$, and reconstruct to obtain $$$\left\{\mathbf{g}_{i, \text{noisy}}\right\}_1^M$$$. We then obtain a sequence of distributions $$$\mathbf{g}_{i,j}$$$ by applying Tikhonov regularization to these perturbed basis elements, across the same range of $$$\lambda_j$$$:

$$\mathbf{g}_{i,j}=\arg\min_{\mathbf{g}_{i,\text{noisy}}\succeq\mathbf{0}}\left\|\mathbf{d}^{(i)}_{\text{obs}}-A\mathbf{g}_{i,\text{noisy}}\right\|_2^2+\lambda_j\left\|\mathbf{g}_{i,\text{noisy}}\right\|_2^2,\quad j=1,2,\cdots,N$$

The noiseless basis elements are represented in terms of the noise-degraded basis elements, defining the effect of noise:

$$\mathbf{g}_{i}\approx\sum_{j=1}^N\beta_{i,j}\mathbf{g}_{i,j}\tag{4}$$

The $$$\mathbf{f}_j$$$ and $$$\mathbf{g}_{i,j}$$$ are obtained with the same $$$\lambda_j$$$, so that we represent each $$$\mathbf{f}_j$$$ as:

$$\mathbf{f}_{j}\approx\sum_{i=1}^N\mathbf{x}_{i,j}\mathbf{g}_{i,j}\tag{5}$$

Combining (2) - (5), we obtain the optimal $$$(\mathbf{c}^*,\alpha^*)$$$ by solving:

$$(\mathbf{c}^*,\mathbf{\alpha}^*)=\arg\min_{\mathbf{c}_i\geq 0,\alpha_j\geq 0}\left\lVert\sum_{j=1}^{N} \alpha_j\sum_{i = 1}^{M}\mathbf{x}_{i,j}\mathbf{g}_{i,j}-\sum_{i=1}^{M}c_i\sum_{j=1}^{N}\beta_{i,j}\mathbf{g}_{i,j}\right\rVert_2^2,\quad\sum_j^N\mathbf{\alpha}_j=1$$

providing the required expansion coefficients, $$$\alpha^*$$$, to determine $$$\mathbf{f}_{\text{true}}$$$. In practice, we perform the above analysis across multiple noise realizations for statistical stability. This approach incorporates the effect of noise on the inverse solution via its impact on the known Gaussian basis through a range of Tikhonov regularizations.

Methods

Simulations: Simulations were tailored to experiments on excised mouse spinal cord and in vivo human thigh muscle in terms of $$$TE$$$, $$$T_2$$$, and noise level.

Mouse Spinal Cord: Multi-echo CPMG (TR/TE = 10s/300μs, 4096 echoes, NEX = 32) spectroscopic data at 9.4T were obtained on formalin-fixed 10mm lengths of cervical and lumbar spinal cord, cross-sectional size ~2x3mm, from a 4-month-old male C57BL/6 mouse using a Bruker Avance III spectrometer, Micro2.5 probe, and 5mm diameter solenoidal coil. Saturation slabs restricted data acquisition to a 2mm slice.

Human Thigh Muscle: After informed consent, data were obtained from a 62-year-old male using a 3T whole-body clinical scanner (Achieva, Philips) with a SENSE Flex-M coil. $$$T_2$$$-weighted scans were collected along the axial plane within the thigh with TE/TR=6ms/5 sec, 72-echo train, in-plane resolution of 3x3mm reconstructed to isotropic 0.98 mm, and 10 mm slice thickness. The data were collected before and after 45-sec intense quadriceps extension exercise.

Results

Figures 1 and 2: Simulations based on spinal cord and muscle experiments. Recovery of known simulated distributions using Multi-Reg, L-curve and GCV are shown. Multi-Reg outperforms the L-curve, and demonstrates increased stability with respect to noise compared with GCV (not shown), which often displays a flat optimization curve.

Figure 3: Distributions from excised spinal cord, with results consistent with previous work ^{4, 5}; the shortest $$$T_2$$$ component is assigned to myelin water, with intermediate and longer components assigned to extra- and intra-axonal water.

Figure 4: Distributions recovered in vivo from the adductor muscle group within the thigh (left: before exercise; right: 12 min after exercise). Two water pools are identified; the more rapidly relaxing pool with $$$T_2$$$~40 ms is assigned to intracellular water and the pool with $$$T_2$$$>200 ms is assigned to extracellular water ⁶. The increase in $$$T_2$$$, from ~260 to ~300 ms, and fraction size of the extracellular pool reflects the known increase in water content after exercise.

Discussion

Acknowledgements

References

1. Sabett, C., Hafftka, A., Sexton, K., & Spencer, R. G. (2017). L1, Lp, L2, and elastic net penalties for regularization of Gaussian component distributions in magnetic resonance relaxometry. Concepts in Magnetic Resonance Part A, 46(2), e21427.

2. Hansen, P. C., & O’Leary, D. P. (1993). The use of the L-curve in the regularization of discrete ill-posed problems. SIAM Journal on Scientific Computing, 14(6), 1487-1503.

3. Wahba, G. (1982). Constrained regularization for ill-posed linear operator equations, with applications in meteorology and medicine. In Statistical decision theory and related topics III (pp. 383-418).

4. Wachowicz, K., & Snyder, R. E. (2002). Assignment of the T2 components of amphibian peripheral nerve to their microanatomical compartments. Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, 47(2), 239-245.

5. Dortch, R. D., Harkins, K. D., Juttukonda, M. R., Gore, J. C., & Does, M. D. (2013). Characterizing inter‐compartmental water exchange in myelinated tissue using relaxation exchange spectroscopy. Magnetic resonance in medicine, 70(5), 1450-1459.

6. Saab, G., Thompson, R. T., Marsh, G. D., Picot, P. A., & Moran, G. R. (2001). Two‐dimensional time correlation relaxometry of skeletal muscle in vivo at 3 Tesla. Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, 46(6), 1093-1098.