On Using Signal Magnitude in Diffusion Magnetic Resonance Measurements of Restricted Motion

TL;DR

This paper demonstrates that relying solely on signal magnitude in diffusion MRI leads to information loss, impairing accurate inference of tissue microstructure, especially near reflective boundaries.

Contribution

It reveals the limitations of magnitude-based diffusion MRI models and highlights the importance of preserving phase information for microstructural analysis.

Findings

Signal magnitude causes loss of Hermitian symmetry in diffusion MRI data.

Loss of symmetry impairs microstructural inference near reflective boundaries.

Simulation shows magnitude-based models are less accurate in complex geometries.

Abstract

Tissue microstructure has significance as a biomarker, however its accurate inference with diffusion magnetic resonance (MR) is still an open problem. With few exceptions, diffusion weighted (DW) MR models either process diffusion MR data using signal magnitude, whereby microstructural information is forcefully confined to symmetry due to Fourier transform properties, or directly use symmetric basis expansions. Herein, information loss from magnitude utilization is demonstrated by numerically simulating particles undergoing diffusion near a fully reflective infinite wall and an orthogonal corner. Simulation results show that the loss of the Hermitian property when using signal magnitude impedes DW--MR from accurately inferring microstructural information in both of the geometries.

Equations

$$x_i(t)=x_i(t_0)+w_i(t)$$

$$M(t)=\sum _i{m}_i(t)=m_0\sum _i{e}^{-ȷ\gamma {\Omega }_i}.$$

$${\Omega }_i=G_{\mathrm{D}}\cdot \left(\underset{t_{d3}}{\overset{t_{d4}}{\int }}{x}_i(\tau )d\tau -\underset{t_{d1}}{\overset{t_{d2}}{\int }}{x}_i(\tau )d\tau \right).$$

$$W_i^{\mathrm{d}}\doteq \underset{t_{d3}}{\overset{t_{d4}}{\int }}{x}_i(\tau )d\tau -\underset{t_{d1}}{\overset{t_{d2}}{\int }}{x}_i(\tau )d\tau =\underset{t_{d3}}{\overset{t_{d4}}{\int }}{w}_i(\tau )d\tau -\underset{t_{d1}}{\overset{t_{d2}}{\int }}{w}_i(\tau )d\tau \in {\mathbb{R}}^3.$$

$$S_{\mathrm{cfd}}^{\mathrm{nmr}}=m_0\sum _i{e}^{-ȷ\gamma G_{\mathrm{D}}\cdot W_i^{\mathrm{d}}}.$$

$$S_{\mathrm{cfd}}^{\mathrm{nmr}}=\int P_{\mathrm{cfd}}^{\mathrm{nmr}}(W^{\mathrm{d}})e^{-ȷ\gamma G_{\mathrm{D}}\cdot W^{\mathrm{d}}}d{W}^{\mathrm{d}},$$

$$S_{\mathrm{cfd}}^{\mathrm{nmr}}(k_{\mathrm{D}})=\mathcal{F}\{P_{\mathrm{cfd}}^{\mathrm{nmr}}\}(k_{\mathrm{D}}),$$

$$P_{\mathrm{cfd}}^{\mathrm{nmr}}(W^{\mathrm{d}})={\mathcal{F}}^{-1}\{S_{\mathrm{cfd}}^{\mathrm{nmr}}\}(W^{\mathrm{d}}).$$

$$s=\sqrt{6D_{\mathrm{water}}{t}_{\mathrm{step}}}=0.48\mu \text{mbox}m.$$

$$x_i(k+1)=x_i(k)+s\eta (k),$$

$$W_i^{\mathrm{d}}=t_{\mathrm{step}}\left(\left(\sum _{k=K_3}^{K_4}{w}_i(k)\right)-\left(\sum _{k=K_1}^{K_2}{w}_i(k)\right)\right)$$

$$S_{\mathrm{cfd}}\to P_{\mathrm{cfd}}\to \text{mbox}Biologicalproperties.$$

$${\mathcal{F}}^{-1}\{|S_{\mathrm{cfd}}^{\mathrm{nmr}}|\}\ne P_{\mathrm{cfd}}^{\mathrm{nmr}}={\mathcal{F}}^{-1}\{S_{\mathrm{cfd}}^{\mathrm{nmr}}\},$$

Full text

Herein, information loss from magnitude utilization is demonstrated by numerically simulating particles undergoing diffusion near a fully reflective infinite wall and an orthogonal corner.

Simulation results show that the loss of the Hermitian property when using signal magnitude impedes DW–MR from accurately inferring microstructural information in both of the geometries.

††journal: Magnetic Resonance Imaging

1 Introduction

Tissue microstructure, which can be inferred by molecular motion measurements, has high significance as a biomarker. However, accurate inference of restricted motion with diffusion MR is still an open problem. Several models are in use [1], and with few exceptions [2, 3], DW–MR models process diffusion MR data using signal magnitude thereby confining microstructural information to symmetry due to Fourier transform properties.

Symmetry assumptions might be realistic in isotropic environments such as homogeneous liquids. By contrast, in a complex environment, such as biological tissue, the path traveled by the molecules is determined by the heterogeneous microstructure and thereby displacement asymmetry becomes possible. Two existing basis expansion models, the mean apparent propagator (MAP) model [2] and generalized diffusion tensor [3] model, recognize asymmetry and work on the complex valued DW–MR signal. Likewise, the more recently developed complete Fourier direct magnetic resonance imaging (CFD–MRI) [4] methodology treats the complex valued signal and the distribution of the displacement integrals without any symmetry constraints by establishing a high dimensional Fourier relationship between them.

However, direct utilization of symmetric templates in model matching methods such as symmetric tensors of diffusion tensor imaging (DTI) [5, 6] or spherical harmonics of high angular resolution diffusion imaging (HARDI) [7], and/or symmetry by way of signal magnitude in the Fourier based $q$–space methodologies [8, 9, 10] is now standard. This raises the concern that information is being automatically distorted and/or lost by collapsing the asymmetry within the signal, thereby impeding accurate inference. As analytical solutions of DW–MR signal processing exist only for simple geometries such as infinite parallel plates [11, 12], planar, cylindrical and spherical geometries [13], [8]; a proof of principle is needed for emphasizing and revising the importance of concerns regarding loss of information when using symmetry based methods.

Herein, DW–MR inference inaccuracies when using signal magnitude are demonstrated by simulating MR signal from diffusing molecules near an infinite reflective wall and an orthogonal corner. While these two basic geometries are not necessarily representatives of biological tissue’s microstructure, they provide solid test beds for analyzing and comparing complex valued data versus magnitude usage in the treatment of diffusion MR signal originating from asymmetric environments.

For properly demonstrating the fundamental concept, the simulations were comprised of pulsed gradient spin echo (PGSE) nuclear magnetic resonance (NMR) experiments without any imaging gradients.

2 Material and Methods

2.1 Theory

For modeling coherent or incoherent motion, the time–dependent position of magnetic moments is represented in the most general fashion [4]:

[TABLE]

using non–interacting particles. In Eq. (1), $w_{i}(t)\in\mathbb{R}^{3}$ represents the displacement of the $i^{\mathrm{th}}$ magnetic moment from its initial position $x_{i}(t_{0})$ (i.e., $w_{i}(t_{0})=0$). The initial time $t_{0}$ is chosen as the end time of the $\displaystyle{\pi/2}$ radio frequency (RF) pulse (see Fig. 1). The model covers any type of motion with only the assumption of path continuity since a magnetic moment cannot disappear at a given point and reappear at another.

The NMR signal, on the other hand, is modeled as the sum of the individual transverse magnetization vectors [14], $m_{i},$ which is written in complex number form:

[TABLE]

Here, $\gamma$ is the gyromagnetic ratio and $m_{0}$ is the initial magnetization tipped to the transverse plane at the end of the $\displaystyle{\pi/2}$ RF pulse at $t_{0}$ (see Fig. 1). The phase is obtained by multiplying $\gamma$ and $\Omega_{i},$ which is ideally a function of the magnetic field gradients and the position of the magnetic moment $x_{i}\in\mathbb{R}^{3}.$ When ideal rectangular gradient pulses with amplitude vector $G_{\mathrm{D}}\in\mathbb{R}^{3}$ as in Fig. 1 are applied, $\Omega_{i}$ becomes:

[TABLE]

where $t_{d1},$ $t_{d3}$ represent on– and $t_{d2},$ $t_{d4}$ off–times of the pulsed gradients.

Following the derivations in [4] that use the time parameters of Fig. 1, the DW–MR signal model’s central theme is the displacement integral. For the PGSE sequence of Fig. 1 the displacement integral of the $i^{\mathrm{th}}$ magnetic moment is defined as [4]:

[TABLE]

Additionally, the PGSE parameters $\delta$ (pulse length) and $\Delta$ (pulse separation) are defined as: $\delta\doteq t_{d2}-t_{d1}=t_{d4}-t_{d3}$ and $\Delta\doteq t_{d3}-t_{d1}=t_{d4}-t_{d2}$ (see also Fig. 1). This is basically the difference of the integrals of the paths traveled by magnetic moments during the pulsed gradients.

According to Eq. (2) and applying Eq. (3) for the pulsed gradients of the PGSE sequence of Fig. 1, the total magnetization at the echo time $T_{E}$ is a function of the pulsed gradient vector (i.e. the motion sensitizing magnetic field gradient vector), $G_{\mathrm{D}},$ and the displacement integral $W^{\mathrm{d}}_{i}$

[TABLE]

Rather than computing Eq. (5) over the (finite) count of magnetic moments, a more convenient expression is obtained by defining $P_{\mathrm{cfd}}^{\mathrm{nmr}}(W^{\mathrm{d}})$ as the fraction of magnetic moments with displacement integral value equal to $W^{\mathrm{d}}$ (per unit interval, $dW^{d}$). Using the distribution of displacement integrals, $P_{\mathrm{cfd}}^{\mathrm{nmr}}(W^{\mathrm{d}}),$ $S_{\mathrm{cfd}}^{\mathrm{nmr}}$ at the echo time $T_{E}$ in Fig. 1 is calculated in the (continuous) displacement integral space for obtaining the signal $S_{\mathrm{cfd}}^{\mathrm{nmr}}:$

[TABLE]

where both sides of the equation are divided by $m_{0}$ for ease of notation. This formulation is general and existing models can in fact be derived starting from Eq. (6) [15, 1].

The expression of Eq. (6) is exactly the Fourier transform of $P_{\mathrm{cfd}}^{\mathrm{nmr}}$ evaluated at the three dimensional frequency vector $k_{\mathrm{D}}:$

[TABLE]

where $k_{\mathrm{D}}=[k_{\mathrm{D}1},\,k_{\mathrm{D}2},\,k_{\mathrm{D}3}]=\gamma\,G_{\mathrm{D}}\in\mathbb{R}^{3}$ represents the pulsed (or diffusion/motion sensitizing) gradient vector with its entries equal to the amplitudes of the pulsed magnetic field gradients of PGSE [4] scaled by $\gamma.$

In other words, for each point $k_{\mathrm{D}}\in\mathbb{R}^{3},$ the Fourier transform of $P_{\mathrm{cfd}}^{\mathrm{nmr}}$ is sampled at the corresponding acquisition. Consequently, $P_{\mathrm{cfd}}^{\mathrm{nmr}},$ which is the physical quantity of interest, is recovered with the inverse Fourier transform of the signal obtained from the MR scanner:

[TABLE]

Using the properties of the Fourier transform, Eq. (8) provides crucial information on the nature of $S_{\mathrm{cfd}}^{\mathrm{nmr}}.$ The Fourier transform of a real valued function is a complex valued Hermitian function [16]. By the reciprocity property of the Fourier transform [16] only symmetric real valued functions’ Fourier transforms are real valued. Therefore, $S_{\mathrm{cfd}}^{\mathrm{nmr}}(k_{\mathrm{D}})$ in Eq. (7) is complex valued Hermitian since $P_{\mathrm{cfd}}^{\mathrm{nmr}}(W^{\mathrm{d}})$ is real valued. Furthermore, when there are coherent displacements such as bulk motion, Eq. (4) incorporates a bias term which in turn reflects as a phase shift in $S_{\mathrm{cfd}}^{\mathrm{nmr}}(k_{\mathrm{D}})$ [4].

By contrast, existing $q$–space models [8, 9, 10] use the magnitude of the signal in practice when evaluating the Fourier transform with the hope of filtering out bulk motion [9, pp.1378]. However, experimental DW–MR signal is not necessarily Hermitian symmetric as various factors including but not limited to susceptibility distort it, whereby the signal magnitude becomes asymmetric. When Fourier transformed, the transform of the asymmetric signal magnitude is a complex valued Hermitian function and therefore cannot describe a physically meaningful distribution. The solution offered by the $q$–space methods is to take the magnitude one more time creating a real symmetric function due to Hermitian symmetry in the complex domain.

In short, by taking the magnitude before and after the Fourier transform, $q$–space methods constrain the outcome to be real valued and symmetric. This raises the concern that more than what was initially intended (e.g., bulk motion removal) is in fact being eliminated unnecessarily.

As it is impossible to fully analytically describe the distribution functions of the displacement integrals, herein, simulations were used to obtain the functions numerically for showing the detrimental effects of using signal magnitude in the calculations.

2.2 Simulations

Scenarios imitating water protons near reflective walls were constructed with $n_{\mathrm{particles}}=120000$ particles. With a fixed simulation step size (the standard deviation of the increments), non–interacting particle motion was implemented using a hindered random walk with normally distributed increments ($w_{i}(t)$ in Eq. (1)) with zero mean. This procedure simulated molecular motion propelled by thermal energy (see [17] for an early simulation example). For modeling hinderance caused by the structure, the displacements were modified when a path crossed a wall. Effectively, the path was readjusted by computing the outcome of the elastic collision(s).

There were $3000$ simulation steps executed for each particle. For matching the experimental values of PGSE parameters used in an earlier biological phantom work [18] $(\delta=15\mbox{ms}$ and $\Delta=30\mbox{ms}),$ simulation values of the same parameters were chosen as $\delta_{\mathrm{sim}}=1000$ steps and $\Delta_{\mathrm{sim}}=2000$ steps with a step time of $\displaystyle{t_{\mathrm{step}}=15\times 10^{-6}\,\mbox{s}},$ e.g. $\delta=t_{\mathrm{step}}\,\delta_{\mathrm{sim}}.$

The step size, $s,$ of the random walk was calculated based on the diffusion coefficient of water, $D_{\mathrm{water}}=2.6\times 10^{-3}\,\mbox{mm}^{2}/\mbox{s}:$

[TABLE]

Numerical simulations consisted of:

DW–MR Simulations

At the initial time step, $k=0,$ uniformly distributing initial positions of the particles near the walls, 2. 2.

Calculating the position of the $i^{\mathrm{th}}$ magnetic moment in discrete time propelled by thermal energy:

[TABLE]

where $\eta(\cdot)\in\mathbb{R}^{3}$ is (pseudorandom) normally distributed. 3. 3.

If the line through $(x_{i}(k+1),x_{i}(k))$ crosses a wall, i.e. when hitting a wall, correcting the path by elastic collision computations, 4. 4.

Evaluating numerically the displacement integral differences of Eq. (4) using the displacements, $w_{i}(k):$

[TABLE]

where t_{\mathrm{step}}\times\left[\matrix{K_{1}&K_{2}&K_{3}&K_{4}}\right]=\left[\matrix{t_{d1}&t_{d2}&t_{d3}&t_{d4}}\right]. 5. 5.

Repeating Step 2 through Step 4 for all of the magnetic moments and 6. 6.

Computing the distribution, $P_{\mathrm{cfd}}^{\mathrm{nmr}},$ of the displacement integral values.

In–house Matlab®, (ver. 2014b, Mathworks, Natick, MA USA) programs were used for all of the computations and visualization. Distributions with $128$ bins were calculated using hist and hist3 functions for one and two dimensional simulations respectively. For implementing discrete Fourier transform (DFT) Matlab®’s fast Fourier transform (FFT) routines were used.

3 Results and Discussion

Although there are solutions in simple microstructures with planar, spherical and cylindrical geometries [11, 12, 13, 8], a full description of information loss from signal magnitude usage in a general geometry is difficult to obtain due to the lack of analytical characterizations of the distribution functions. Herein, Brownian motion near reflecting walls in one and two dimensions, similar to the work in [19, 20], was simulated with the algorithm of Section 2.2 for numerically calculating the distributions.

For one dimensional scenario, a wall was placed on the origin and particles were initially uniformly distributed in an interval starting from the origin ending at $5$ different positions: $(0\;6\mu m],\,(0\;13\mu m],\,(0\;20\mu m],\,(0\;27\mu m],\,(0\;34\mu m]$ (see Fig. 2). Simulations were executed and $P_{\mathrm{cfd}}^{\mathrm{nmr}}(W^{\mathrm{d}})$ was obtained from the numerically calculated $W^{\mathrm{d}}_{i},\>i=1,\ldots,n_{\mathrm{particles}}.$

Figure 2 compares distribution functions obtained from different initial conditions. Initial condition intervals extending farther from the wall resulted in more symmetric and spread $P_{\mathrm{cfd}}^{\mathrm{nmr}},$ resembling more isotropic motion. In contrast, $P_{\mathrm{cfd}}^{\mathrm{nmr}}$’s asymmetry is more prominent when particles start near the wall.

For all of the distributions, $\mathcal{F}\{P_{\mathrm{cfd}}^{\mathrm{nmr}}\}$ was calculated for simulating DW–MR signal coming out of the MR scanner. As all of the distributions were asymmetric, their Fourier transforms were complex valued Hermitian. Therefore, the magnitude, $|\mathcal{F}\{P_{\mathrm{cfd}}^{\mathrm{nmr}}\}|,$ was a real valued symmetric function. As the (inverse) Fourier transform of a real symmetric function is also real symmetric [16], $\mathcal{F}^{-1}\{|\mathcal{F}\{P_{\mathrm{cfd}}^{\mathrm{nmr}}\}|\}$ was real and symmetric. This is demonstrated on the right of Fig. 2 by showing the mismatch between $P_{\mathrm{cfd}}^{\mathrm{nmr}}$ (asymmetric) and $\mathcal{F}^{-1}\{|\mathcal{F}\{P_{\mathrm{cfd}}^{\mathrm{nmr}}\}|\}$ (symmetric).

In the one dimensional scenario which did not incorporate any bulk motion, using signal magnitude, which was claimed to be filtering out bulk motion [9, pp.1378], should have left the distribution intact. Figure 2 shows that distribution asymmetry is lost and the symmetric distribution obtained from signal magnitude implies wrongly an isotropic medium. In reality, for particles near an elastic wall the motion is anisotropic as reported by the asymmetry of the original distribution.

In two dimensions, two scenarios were run using a vertical wall and a corner, both reflective and placed at the origin (see Fig. 3). For the vertical wall, the initial positions were uniformly distributed on the $(0,\;20\mu m]\times(0,\;5.76mm]$ strip ($5.76mm=0.1\times s\times n_{\mathrm{particles}}$) and for the corner, the initial distribution was uniform on the square: $(0,\;20\mu m]\times(0,\;20\mu m]$ (Fig. 3, left panel).

In Fig. 3’s middle panel, accurate inference of the wall geometry was attained as original distributions were analyzed. By contrast, on the right panel, the Fourier transformed magnitude produces symmetric distributions with the appearance of a tube (top) and a box (bottom) rather than a vertical wall and a corner respectively.

Also, taking the magnitude transfers the imaginary portion to the real axis, thereby adding a positive constant to the signal which in turn becomes a spike at $[0,\,0]$ in Fig. 3 and at [math] in Fig. 2 when calculating the Fourier transform. This, in turn, inaccurately indicates a large number of magnetic moments with zero displacement integral values, which might mislead inference of mathematical models that include terms for stationary magnetic moments, e.g., of water trapped in small compartments [21, 22], likewise the positive constant term in models using cumulant expansion [23].

4 Conclusions

Biological interpretation of DW–MR information is the subject of active research with new developments using existing model attributes, such as the decrease in the apparent diffusion coefficient caused by neurite beading after ischemic stroke [24] and cerebral cortical gray matter development [25]. New theories describing novel attributes are also being developed for the same purpose, e.g., the MAP model [2] where the focus is the microstructure in brain tissue.

The purpose of the DW–MR, like any other imaging modality, is the inference of biological properties with the aim of producing informative biomarkers. This is accomplished in two steps starting from the MR scanner’s output signal that is Fourier transformed for reconstructing the physically meaningful displacement integral distribution:

[TABLE]

In the process, intermediate steps must be carried out carefully while respecting physically meaningful information content.

This work demonstrated that, as shown in Fig. 2 and Fig. 3, the usage of magnitude may result in incorrect information:

[TABLE]

potentially hampering accurate inference of microstructure.

The proof was achieved in an ideal simulated environment with minimal assumptions, e.g. without including loss of magnetization at the walls [26, 27], wall permeability [28] and particle interactions. While no mathematical model is perfect, falsification when validating under ideal conditions is a legitimate reason for abandoning a model or at least using prudence in interpreting the model’s outcomes. For future in and ex vivo diffusion imaging research, rather than assuming that diffusion profiles in biological tissue are symmetric for any reason, it would be prudent to adopt DW–MR models incorporating asymmetry assumptions and let the data analysis indicate symmetry and/or lack thereof.

Accordingly, instead of using signal magnitude, the solution offered by CFD–MR [4] is based on re–establishing Hermitian symmetry of the complex valued signal collected from the MR scanner by using systemic phase correction algorithms. For example, bulk motion appears as a linear phase which is rigorously handled with a linear phase correction algorithm in the Fourier domain [4]. However, factors which are beyond the scope of this manuscript, including but not limited to susceptibility effects, distort the Hermitian symmetry of the DW–MR signal in a nonlinear fashion. Development of algorithms for analyzing and subsequently correcting those effects remains an open question for future research. Furthermore, using complex valued signal in CFD–MR automatically removes nonlinear Rician restrictions [29] on noise modeling in DW–MR.

In conclusion, this paper mainly concentrated on the first step that constitutes the foundation for biological property inference, namely removing signal magnitude processing from DW–MR methods with the aim of improving accurate reconstruction of physical quantities. Simulations in fundamental geometries demonstrated that the usage of signal magnitude in diffusion MR can significantly distort the estimation of motion describing distributions, and hence may impede accurate inference of microstructure.

Acknowledgements

This work was supported by the Scientific and Technological Research Council of Turkey’s (TUBITAK) $3001$ funding program, titled “Simulating Microstructure Inferred from Diffusion Magnetic Resonance Imaging and Modelling Effects of the MR Scanner on the Diffusion MR Signal”, grant number $118$E$304.$

The author would also like to extend his special thanks to John Gore and the reviewers for their valuable insights and progressive feedback.