A diffusion MRI model for random walks confined on cylindrical surfaces: towards non-invasive quantification of myelin sheath radius

Introduction: Quantifying the myelin sheath radius of myelinated axons in vivo is important for understanding, diagnosing, and monitoring various neurological disorders. Despite advancements in diffusion MRI (dMRI) microstructure techniques, there are currently no models specifically designed to estimate myelin sheath radii.

Methods: This proof-of-concept theoretical study presents two novel dMRI models that characterize the signal from water diffusion confined to cylindrical surfaces, approximating myelin water diffusion. We derive their spherical mean signals, eliminating fiber orientation and dispersion effects for convenience. These models are further extended to account for multiple concentric cylinders, mimicking the layered structure of myelin. Additionally, we introduce a method to convert histological distributions of axonal inner radii from the literature into myelin sheath radius distributions. We also derive analytical expressions to estimate the effective myelin sheath radius expected from these distributions.

Results and Discussion: Monte Carlo (MC) simulations conducted in cylindrical and spiral geometries validate the models. These simulations demonstrate agreement with analytical predictions. Furthermore, we observe significant correlations between the effective radii derived from histological distributions and those obtained by fitting the dMRI signal to a single-cylinder model. These models may be integrated with existing multi-compartment dMRI techniques, opening the door to non-invasive in vivo assessments of myelin sheath radii. Such assessments would require MRI scanners equipped with strong diffusion gradients, allowing measurements with short echo times. Further work is required to validate the technique with real dMRI data and histological measurements.

1 Introduction

White matter (WM) primarily consists of axons [1], which are often enveloped by myelin produced by oligodendrocytes [2]. Myelin serves as an insulating sheath that enables nerve signals to propagate faster along the axon [3, 4]. The axon-myelin unit interacts through complex molecular signaling and cellular processes, regulating the development and maintenance of myelin and the overall axon radius. Disruptions in the axon-myelin unit, such as demyelination or axon damage, are associated with neurological disorders such as multiple sclerosis [5], severe psychiatric conditions [6, 7], and Alzheimer’s disease [8]. These disorders are known to impair diverse cognitive functions [9]. Quantifying the microstructural properties of myelinated axons in vivo is crucial for enhancing our understanding of neurological diseases. This will ultimately improve diagnosis, early disease detection, and treatment of neurological disorders that affect millions worldwide.

Magnetic Resonance Imaging (MRI) is the primary technique for in vivo, non-invasive imaging of WM in the human brain. Many MRI techniques have been developed to characterize distinct WM properties [10–12]. For example, diffusion-weighted MRI (dMRI) measures the random motion of water molecules within and around axons. This sensitivity enables the estimation of spatial maps for various WM characteristics, such as axon orientations [13–30], dispersion [31, 32], axon volume fraction [33–35], inner axon radii [10–12, 36–46], intra- and extra-axonal water diffusivities [47, 48], and T2 relaxation times [49, 50]. In contrast, multi-echo T2 relaxometry [51–62] provides estimates closely correlated with myelin volume.

Despite considerable progress, challenges and research gaps remain in estimating the full range of WM microstructural features. One of them is the absence of specialized dMRI models explicitly designed for in vivo estimation of myelin sheath radii. Understanding water diffusion dynamics within myelin bilayers is essential, as the “apparent” radial diffusivity of myelin water likely depends on the myelin sheath radius. This connection is promising, as it could enable myelin sheath radius estimation using dMRI data.

Accurately estimating myelin water diffusivities is challenging. This is because myelin water contributes minimally to the dMRI signal due to its short T2 time (i.e., 15 ms [52]), compared to the longer echo times (TE∼80 ms) used in standard dMRI sequences. Nevertheless, various ex-vivo studies attempted to estimate myelin water diffusivities using T2 and T1 relaxation selective measurements. A diffusion-relaxation hybrid experiment proposed by [63], using a Carr-Purcell-Meiboom-Gill sequence, surprisingly revealed minor diffusional anisotropy and large parallel and radial diffusivities for the short T2 component associated with myelin water in the bovine optic nerve. Another approach employed T2 relaxation time to characterize myelin water selectively in the frog’s peripheral nerve [64]. However, this ex-vivo study did not report myelin water diffusivities. On the other hand, T1 and T2 relaxation times have been utilized to observe myelin water in the excised frog sciatic nerve [65]. The T1-based method employed a double inversion recovery (DIR) sequence to nullify non-myelin water components, resulting in signals predominantly (>90%) derived from myelin water. This study found that myelin water diffusivities were lower when selected based on T1 characteristics with DIR-T1 measures (yielding parallel and radial diffusivities of D_∥ = 0.37–0.43 μm²/s and D_⟂ = 0.13–0.17 μm²/ms, respectively) compared to T2 characteristics (D_∥ = 0.8 μm²/s and D_⟂ = 0.19 μm²/ms).

Conversely, various in-vivo human brain studies have attempted to make the dMRI signal sensitive to the microstructure of myelin tissue. For instance [66], implemented a magnetization transfer (MT) prepared stimulated-echo diffusion tensor imaging technique. The short TE = 34 ms enabled by the stimulated-echo acquisition preserved a significant signal from the myelin water component with short T2, while the MT preparation further provided differentiating sensitization to this signal. Compared to the diffusion tensor derived from the conventional dMRI sequence acquired without MT preparation, the myelin water weighted tensor exhibited a significant increase in fractional anisotropy, most likely explained by the lower radial diffusivity of myelin water. In recent years, the diffusion-T2 relaxation approach has gained momentum thanks to the emergence of human scanners with strong diffusion gradients G [67–69], allowing the use of diffusion sequences with shorter TEs. TE can be further reduced by using dMRI sequences with spiral readouts; for example, in the work by [70, 71], TEs of 21.7 and 30 ms were achieved for b = 1,000 and 6,000 s/mm² respectively, with G = 300 mT/m, whereas [72] reduced the TE to 19 ms for b = 1,000 s/mm² with G = 200 mT/m.

These recent studies suggest that it is possible to acquire dMRI data significantly weighted by myelin water. Therefore, this is an opportune time to develop new dMRI models for this often-overlooked WM compartment. In this theoretical and numerical proof of concept study, we propose a novel dMRI model for the diffusion process within a series of impermeable concentric cylinders separated by infinitesimal gaps filled with water, which could be employed as a first approximation to estimate myelin sheath radius. We derive the analytical dMRI signal and a Gaussian approximation with time-dependent radial diffusivity for this geometrical model and used Monte Carlo (MC) diffusion simulations to validate the proposed models.

This article is organized as follows. Section 2 presents our study’s mathematical derivations, beginning with the geometrical model for the diffusion process in multiple concentric cylinders separated by infinitesimal distances (Section 2.1). We then model the dMRI signal as the product of signals generated by displacements parallel and perpendicular to the main cylinder’s axis (Section 2.2) and introduce the diffusion propagator formalism to derive the analytical dMRI signal under the narrow-pulse approximation for pulsed-gradient spin-echo (PGSE) acquisitions (Section 2.3). A Gaussian approximation is presented in Section 2.4, followed by a refinement of these models in Section 2.5 to account for PGSE sequences with rectangular or trapezoidal diffusion gradients with non-narrow pulses. In Section 2.6, we derive the spherical mean signals, simplifying the modeling by eliminating fiber orientation and dispersion effects. In Section 2.7, we explore theoretical approximations to clarify how the estimated cylinder radius should be interpreted when fitting these models to measured data. The Methods section (Section 3) details the dMRI MC simulations designed to validate the proposed models. The results are presented in Section 4, followed by a discussion of their significance and the study’s limitations in Section 5.

2 Theory

2.1 General description – geometrical model

Oligodendrocytes extend their cell membranes to wrap around axons in WM, creating multiple concentric layers of myelin. Each turn of wrapping adds another bilayer of myelin with a thickness of approximately d_m = 4–5 nm. This process results in a multilayer spiral structure, with gaps of about d_w = 3 nm thick [73] between the layers, filled by myelin water. Figure 1A shows a schematic transverse section of a myelinated axon.

Figure 1

Figure 1. Schematic representation of an axon and its myelin sheath. (A) Cross-sectional view of a myelinated axon showing the spiral trajectory of compact myelin bilayers (in yellow-orange). Each myelin bilayer has a thickness of approximately 4–5 nm and is separated by myelin water gaps (i.e., cytoplasmic and extracellular water) (in blue) with a thickness of approximately 3 nm [73]. (B) Cross-section of multiple concentric alternating cylinders representing the myelin bilayers and myelin water. This simplified geometrical model is used to study the diffusion process. (C) Example of myelin water molecules (represented by blue dots) diffusing on a cylindrical surface, where denotes the polar angle, quantifying the displacement of a water molecule along the 2D surface in the x-y plane. This plane is assumed to be perpendicular to the main axis of the cylinder, which is oriented along the z-axis.

In this study, we approximate the diffusion process along this spiral trajectory as diffusion within a series of impermeable concentric solid cylinders separated by infinitesimal water-filled gaps (see Figure 1B). The rationale for this approximation is as follows: For a given diffusion time, a diffusing water molecule traveling a total displacement of 2πaN (where a is the myelin radius at the starting position and N is an arbitrary number) along the spiral trajectory experiences a net radial displacement of about N (d_w + d_m) (see cross-sectional plane shown in Figure 1A. This displacement remains negligible, even for molecules traveling long distances. For example, for a = 0.5 µm and N = 10, the path length along the spiral is 31.4 µm, and the net radial displacement is approximately 0.08 µm, hence significantly smaller than the minimum displacement required to attenuate the dMRI signal in state-of-the-art scanners [39, 74]. Moreover, since spin echo dMRI sequences designed to be sensitive to myelin water employ short TEs (equivalently short diffusion times), most molecules will travel relatively short distances along the spiral trajectory, minimizing the net radial displacement.

For this reason, we propose to simplify the spiral trajectory by using concentric cylinders of similar size. As infinitesimal distances separate the cylinders, we assume that the underlying diffusion process is equivalent to random walks confined to the cylinder surfaces. Therefore, we will first derive the diffusion propagator for Brownian motion on the cylinder surface, see Figure 1C, and then extend this model to multiple cylinders. Moreover, to eliminate fiber orientation and dispersion effects (confounding factors), we will derive the spherical mean dMRI signal for this model. This approach will help us to interpret the mean radius estimated by fitting a single-cylinder-surface model to the dMRI signal arising from multiple cylindrical surfaces.

2.2 Decoupling diffusive motions

To simplify our model, we will consider an infinitely long cylinder whose main axis is oriented along the z-axis, with its transverse section lying in the x-y plane. An important aspect of this model is that the dMRI signal can be decomposed into contributions from spin particles diffusing parallel and perpendicular to the cylinder’s main axis. In this coordinate frame of reference, these diffusion processes are statistically independent. Therefore, the displacement probability distribution $P\left(\mathbf{r},t\right)=P\left({\mathbf{r}}_{xy},t\right)P\left({\mathbf{r}}_z,t\right)$ can be expressed as the product of the distributions for motion in the perpendicular $P\left({\mathbf{r}}_{xy},t\right)$ and parallel $P\left({\mathbf{r}}_z,t\right)$ directions. The net displacement vector $\mathbf{r}={\mathbf{r}}_{xy}+{\mathbf{r}}_z$ at diffusion time $t$ can be decomposed into the displacement vectors perpendicularly and parallel to the cylinder’s axis. Note that ${\mathbf{r}}_{xy}=r_x\widehat{\mathbf{i}}+r_y\widehat{\mathbf{j}}$ and ${\mathbf{r}}_z=r_z\widehat{\mathbf{k}}$, where $r_x$, $r_y$ and $r_z$ are the vector’s lengths along the unit vectors $\widehat{\mathbf{i}}$, $\widehat{\mathbf{j}}$, $\widehat{\mathbf{k}}$ associated with the x-, y-, and z-axes, respectively.

For this type of decoupled diffusive motion [75], showed that the dMRI signal can be expressed as the product of the dMRI signals arising from displacement parallel and perpendicular to the cylinder’s axis: $E\left(\mathbf{q},t\right)=E_{\perp }\left({\mathbf{q}}_{xy},t\right)E_{\parallel }\left({\mathbf{q}}_z,t\right)$, where ${\mathbf{q}}_{xy}=q_x\widehat{\mathbf{i}}+q_y\widehat{\mathbf{j}}$ and ${\mathbf{q}}_z=q_z\widehat{\mathbf{k}}$, $\mathbf{q}={\mathbf{q}}_{xy}+{\mathbf{q}}_z=\gamma \mathbf{g}\delta$, $\gamma$ is the gyromagnetic ratio of the diffusing spin particles (e.g., hydrogen nuclei), $\mathbf{g}=G\widehat{\mathbf{g}}$ denotes the applied diffusion gradient with magnitude $G$ and unit orientation vector $\widehat{\mathbf{g}}$, and $\delta$ is the duration of the diffusion gradient pulses. Note that $t$ should be expressed in terms of the dMRI sequence time parameters. A general detailed derivation of this decoupled signal model is provided in [75].

2.3 Diffusion MRI signal and diffusion propagator: narrow-delta approximation

In this section, we will derive the analytical expressions for $E_{\parallel }\left({\mathbf{q}}_z,t\right)$ and $E_{\perp }\left({\mathbf{q}}_{xy},t\right)$ necessary to provide the full dMRI signal model. This derivation follows the diffusion propagator formalism under the narrow-pulse (narrow-delta) approximation, which assumes that the duration of the diffusion gradient is very short ($\delta \to 0$). Thus, under this formalism and for a PGSE sequence [76], the diffusion time is equal to the time difference between the onset of the two diffusion gradients $t=\Delta$.

The dMRI signal $E_{\parallel }\left({\mathbf{q}}_z,t\right)$ arising from displacements parallel to the cylinder’s main axis ${\mathbf{r}}_z$ is related to the 1D displacement probability distribution by the following Fourier-relationship:

$\frac{{\mathit{E}}_{\parallel }\left({\mathbf{q}}_{\mathrm{z}},\mathrm{t}\right)}{{\mathit{E}}_{\parallel }\left({\mathbf{q}}_{\mathrm{z}}=0,\mathrm{t}\right)}={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}\mathit{P}\left(\left.\mathbf{z}\right|,{\mathbf{z}}^{\prime },\mathit{t}\right)\mathrm{P}\left({\mathbf{z}}^{\prime }\right){\mathit{e}}^{\mathit{i}{\mathit{q}}_{\mathit{z}}\left(\mathit{z}-{\mathit{z}}^{\prime }\right)}\mathit{d}{\mathit{z}}^{\prime }\mathit{dz},(1)$

where $P\left({\mathbf{z}}^{\prime }\right)$ is the probability for a particle to be at position ${\mathbf{z}}^{\prime }=z^{\prime }\widehat{\mathbf{k}}$ at initial time $t=0$, and $P\left(\left.\mathbf{z}\right|,{\mathbf{z}}^{\prime },t\right)$ is the probability that a particle initially located at position ${\mathbf{z}}^{\prime }$ migrate to position $\mathbf{z}=z\widehat{\mathbf{k}}$ in time $t$. Assuming that at $t=0$ all particles are uniformly distributed along the cylinder’s axis (i.e., $P\left({\mathbf{z}}^{\prime }\right)$ is constant) and using the change of variables ${\mathbf{r}}_z=\mathbf{z}-{\mathbf{z}}^{\prime }$ to quantify displacements, Equation 1 can be rewritten as

$\frac{E_{\parallel }\left({\mathbf{q}}_z,t\right)}{E_{\parallel }\left({\mathbf{q}}_z=0,t\right)}={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}}P\left({\mathbf{r}}_z,t\right)e^{i{q}_z{r}_z}d{r}_{z,}(2)$

Since the motion of particles along the cylinder’s main axis is unrestricted (assuming an infinitely long cylinder), we assume 1D Gaussian diffusion with a characteristic myelin water diffusivity $D_{\parallel }$ on the cylinder’s surface:

$\mathit{P}\left({\mathbf{r}}_{\mathrm{z}},\mathit{t}\right)=\frac{1}{\sqrt{4\mathrm{\pi }{\mathit{D}}_{\parallel }\mathit{t}}}{\mathit{e}}^{-\frac{{\mathit{r}}_{\mathit{z}}^2}{4{\mathit{D}}_{\parallel }\mathit{t}}}(3)$

By inserting Equation 3 into Equation 2, we obtain the familiar dMRI signal expression for Gaussian diffusion,

$E_{\parallel }\left({\mathbf{q}}_z,t\right)=E_{\parallel }\left({\mathbf{q}}_z=0,t\right)e^{-q_z^2{D}_{\parallel }t}.(4)$

Likewise, the dMRI signal arising from displacements perpendicular to the cylinder’s axis $E_{\perp }\left({\mathbf{q}}_{xy},t\right)$ depends on the 2D displacement probability distribution by the following Fourier-relationship:

$\frac{E_{\perp }\left({\mathbf{q}}_{xy},t\right)}{E_{\perp }\left({\mathbf{q}}_{xy}=0,t\right)}={\displaystyle \underset{{\mathbb{R}}^2}{\int }}{\displaystyle \underset{{\mathbb{R}}^2}{\int }}P\left(\left.{\mathbf{r}}_{xy}\right|,{\mathbf{r}}_{xy}^{\prime },t\right)P\left({\mathbf{r}}_{xy}^{\prime }\right)e^{i{\mathbf{q}}_{xy}\left({\mathbf{r}}_{xy}-{\mathbf{r}}_{xy}^{\prime }\right)}d{\mathbf{r}}_{xy}^{\prime }d{\mathbf{r}}_{xy,}(5)$

where $P\left({\mathbf{r}}_{xy}^{\prime }\right)$ and $P\left(\left.{\mathbf{r}}_{xy}\right|,{\mathbf{r}}_{xy}^{\prime },t\right)$ are the probability of finding a particle at position ${\mathbf{r}}_{xy}^{\prime }$ in the x-y plane at $t=0$, and the probability of moving from ${\mathbf{r}}_{xy}^{\prime }$ to ${\mathbf{r}}_{xy}$ in time $t$.

As the particle displacements in the plane perpendicular to the cylinder’s axis are confined on a circle, it is convenient to rewrite the integrals in Equation 5 in polar coordinates due to the polar symmetry of this system,

$\begin{array}{l}\frac{E_{\perp }\left({\mathbf{q}}_{xy},t\right)}{E_{\perp }\left({\mathbf{q}}_{xy}=0,t\right)}={\displaystyle \underset{0}{\overset{\infty }{\int }}}{\displaystyle \underset{0}{\overset{\infty }{\int }}}{\displaystyle \underset{0}{\overset{2\pi }{\int }}}{\displaystyle \underset{0}{\overset{2\pi }{\int }}}P\left(\left.\rho ,\theta \right|,{\rho }^{\prime },{\theta }^{\prime },t\right)P\left({\rho }^{\prime },{\theta }^{\prime }\right)\\ \times e^{i{q}_{xy}\rho \cos \left(\phi +\theta \right)}{e}^{-i{q}_{xy}{\rho }^{\prime }\cos \left(\phi +{\theta }^{\prime }\right)}\rho {\rho }^{\prime }d{\rho }^{\prime }d\rho d\theta d{\theta }^{\prime },\end{array}(6)$

where the cartesian components of the 2D vectors, ${\mathbf{r}}_{xy}$, ${\mathbf{r}}_{xy}^{\prime }$, ${\mathbf{q}}_{xy}$, are rewritten in terms of their magnitudes, $\rho$, ${\rho }^{\prime }$, $q_{xy}$, and angles of orientation, $\theta$, ${\theta }^{\prime }$, $\phi$, respectively: ${\mathbf{r}}_{xy}=\left(\rho \cos \left(\theta \right),\rho \sin \left(\theta \right)\right)$, ${\mathbf{r}}_{xy}^{\prime }=\left({\rho }^{\prime }\cos \left({\theta }^{\prime }\right),{\rho }^{\prime }\sin \left({\theta }^{\prime }\right)\right)$, and ${\mathbf{q}}_{xy}=\left(q_{xy}\cos \left(\phi \right),q_{xy}\sin \left(\phi \right)\right)$.

In Supplementary Appendix A, we show that Equation 6 can be simplified to

$\frac{E_{\perp }\left({\mathbf{q}}_{xy},t\right)}{E_{\perp }\left({\mathbf{q}}_{xy}=0,t\right)}=\frac{1}{2\pi }{\displaystyle \underset{0}{\overset{2\pi }{\int }}}{\displaystyle \underset{0}{\overset{2\pi }{\int }}}P\left(\left.\mathrm{\Phi }\right|,a,t\right)e^{i{q}_{xy}a\cos \left(\psi \right)}{e}^{-i{q}_{xy}a\cos \left(\psi -\mathrm{\Phi }\right)}d\psi d\mathrm{\Phi },(7)$

where we used the change of variables $\psi =\phi +\theta$ and $\mathrm{\Phi }=\theta -{\theta }^{\prime }$, and $P\left(\left.\mathrm{\Phi }\right|,a,t\right)$ is the probability that the particles’ motion on the circle with radius $a$ covers a polar angle $\mathrm{\Phi }$ at a time $t$, see Figure 1C.

We model $P\left(\left.\mathrm{\Phi }\right|,a,t\right)$ as a wrapped Gaussian distribution [77] with diffusivity $D$:

$\begin{array}{l}P\left(\left.\mathrm{\Phi }\right|,a,t\right)=\frac{1}{\sqrt{4\pi Dt}}{\displaystyle \sum _{p=-\infty }^{\infty }}{e}^{-\frac{a^2}{4Dt}{\left(\mathrm{\Phi }+2\pi p\right)}^2}\\ =\frac{1}{2\pi }\left[1+2{\displaystyle \sum _{p=1}^{\infty }}{e}^{-p^2\frac{Dt}{a^2}}\cos \left(p\mathrm{\Phi }\right)\right].\end{array}(8)$

This distribution results from wrapping the 1D Gaussian distribution (on the infinite line) around the circle’s circumference. It takes into account that during the diffusion process, a population of particles could travel distances larger than $2\pi ap$, where $2\pi a$ is the perimeter of the circle, and $p=1,2,\dots ,\infty$. The second expression in Equation 8 provides a helpful alternative representation of this function [77–79]. It is the solution of the diffusion equation of Brownian particles confined in a circle $S^1$ [80–82]. However, note that in [81], the function was normalized with the circle’s perimeter, whereas our distribution is normalized with the angle, i.e., ${\int }_0^{2\pi }P\left(\left.\mathrm{\Phi }\right|,a,t\right)d\mathrm{\Phi }=1$.

Assuming that the translational diffusion parallel to the cylinder’s main axis and along the “unwrapped” circle are equal, then $D=D_{\parallel }$. After substituting Equation 8 into Equation 7 we obtain

$\frac{{\mathit{E}}_{\perp }\left({\mathbf{q}}_{\mathit{xy}},\mathit{t}\right)}{{\mathit{E}}_{\perp }\left({\mathbf{q}}_{\mathit{xy}}=0,\mathit{t}\right)}={\mathit{J}}_0^2\left(\mathit{a}{\mathit{q}}_{\mathit{xy}}\right)+2{\displaystyle \sum _{\mathit{p}=1}^{\infty }}{\mathit{J}}_{\mathit{p}}^2\left(\mathit{a}{\mathit{q}}_{\mathit{xy}}\right){\mathit{e}}^{-{\mathit{p}}^2\frac{{\mathit{D}}_{\parallel }\mathit{t}}{{\mathit{a}}^2}},(9)$

where $J_p$ is the p-th Bessel function of the first kind. The complete derivation is shown in Supplementary Appendix A. This expression does not depend on the orientation $\phi$ of vector ${\mathbf{q}}_{xy}$ in the plane perpendicular to the cylinder’s axis due to transverse symmetry, as expected. In the limit $D_{\parallel }t\gg a^2$, Equation 8 becomes a uniform distribution and Equation 9 tends to

$\frac{E_{\perp }\left({\mathbf{q}}_{xy},t\gg a^2/{\mathit{D}}_{\parallel }\right)}{E_{\perp }\left({\mathbf{q}}_{xy}=0,t\gg a^2/{\mathit{D}}_{\parallel }\right)}\approx J_0^2\left(a{q}_{xy}\right).(10)$

Note that Equation 10 does not depend on $t$.

An independent derivation of Equation 9 was reported in [64, 83]. However, the result reported by [64] was obtained by assuming a Gaussian distribution for the angular motion instead of a Wrapped Gaussian, which solution only tends to Equation 9 in the limit case when $a^2\gg D_{\parallel }t$.

By merging results from Equations 4, 9, we obtain the final signal model for a single cylinder:

$\begin{array}{l}E\left(\mathbf{q},t\right)=E_{\parallel }\left({\mathbf{q}}_z,t\right)E_{\perp }\left({\mathbf{q}}_{xy},t\right)\\ =E\left(\mathbf{q}=0,t\right)e^{-q^2\cos {\left(\beta \right)}^2{D}_{\parallel }t}\left[J_0^2\left(aq\sin \left(\beta \right)\right)+2{\displaystyle \sum _{p=1}^{\infty }}{J}_p^2\left(aq\sin \left(\beta \right)\right)e^{-p^2\frac{D_{\parallel }t}{a^2}}\right],\end{array}(11)$

where $\beta$ is the angle between the diffusion gradient orientation and the cylinder’s axis, $q_{xy}=q\sin \left(\beta \right)$ and $q_z=q\cos \left(\beta \right)$. For practical purposes, the signal can be adequately approximated by the first $p=1,\dots ,P$ terms in the series.

2.4 Gaussian approximation

When the displacement probability distribution in the x-y plane (perpendicular to the cylinder’s axis) is approximated by an isotropic bivariate Gaussian distribution, the mean-squared displacement of particles $〈{\left|{\mathbf{r}}_{xy}\right|}^2〉$ is related to the ‘apparent’ radial diffusivity in the 2D plane according to $D_{\perp }^{app}=〈{\left|{\mathbf{r}}_{xy}\right|}^2〉/4t$. For such an isotropic Gaussian diffusion process, the corresponding dMRI signal $E_{\perp }\left({\mathbf{q}}_{xy},t\right)$ is given by

$E_{\perp }\left({\mathbf{q}}_{xy},t\right)=E_{\perp }\left({\mathbf{q}}_{xy}=0,t\right)e^{-q_{xy}^2{D}_{\perp }^{app}t}(12)$

The expression for $D_{\perp }^{app}$ in Equation 12 depends on the diffusion time and circle radius $a$ as

$D_{\perp }^{app}=\frac{a^2}{2t}\left[1-{\mathrm{e}}^{-\frac{D_{\parallel }t}{a^2}}\right],(13)$

where we assumed $D=D_{\parallel }$, like in Equation 9. The full derivation is presented in Supplementary Appendix B. For very short diffusion times, $t\to 0$, the apparent radial diffusivity does not depend on the circle’s radius, $D_{\perp }^{app}=D_{\parallel }/2$, because no structural features are probed at such a small time-scale. Conversely, for very long diffusion times, $D_{\perp }^{app}\to a^2/2t$.

The final dMRI signal, considering both the parallel and radial diffusion components, is given by

$E\left(\mathbf{q},t\right)=E\left(\mathbf{q}=0,t\right)e^{-q^2\cos {\left(\beta \right)}^2{D}_{\parallel }t}{e}^{-q^2\sin {\left(\beta \right)}^2{D}_{\perp }^{app}t}(14)$

This analytical form is equivalent to an axially symmetric diffusion tensor signal, as described in Equation 5 in [84]. However, note that the radial diffusivity depends on the diffusion time and the size of the confining geometry, i.e., the cylinder radius.

2.5 Correction for non-narrow deltas

Our previous derivations are based on the q-space formalism (see Equations 1, 5). This approach is valid for PGSE sequences [76] using diffusion-encoding gradients with infinitesimal duration $\delta$. Consequently, the proposed signal models are not valid for sequences that do not fulfill this requirement. In this section, we will use the q-space correction approach presented by [85] to provide more general signal approximations beyond this acquisition protocol.

Under the narrow pulse approximation, the dephasing of the spins due to their motion during the application of the diffusion gradients is neglected. Thus, the diffusion time is equal to the time difference between the onset of the two diffusion gradients. However, for finite $\delta$ it is unclear what diffusion time derived from the PGSE sequence must be used in the diffusion propagator to evaluate the dMRI model. This problem was tackled by [85], who proposed a general relationship between the signal attenuation ${〈e^{i\varphi }〉}_{\Delta ,\delta ,\mathbf{g}}$ for the PGSE sequence and the displacement probability

$\frac{{〈e^{i\varphi }〉}_{\Delta ,\delta ,\mathbf{g}}}{{〈e^{i\varphi }〉}_{\Delta ,\delta ,\mathbf{g}=0}}={\displaystyle \underset{{\mathbb{R}}^3}{\int }}P\left(\mathbf{r},t_{\exp }\right){〈\left.e^{i\varphi }\right|\mathbf{r}〉}_{\Delta ,\delta ,\mathbf{g}}d\mathbf{r},(15)$

where the integral in Equation 15 is over the infinite three-dimensional space, $t_{\exp }$ is the total diffusion time of the experiment between the onset of the first gradient and the termination of the second gradient, and ${〈\left.e^{i\varphi }\right|\mathbf{r}〉}_{\Delta ,\delta ,\mathbf{g}}$ denotes the average signal attenuation (dephasing) of the population of spins experiencing a net displacement $\mathbf{r}$ in time $t_{\exp }$. Note that $t_{\exp }=\Delta +\delta$ for PGSE sequences with rectangular diffusion gradients and $t_{\exp }=\Delta +\delta +\xi$ for trapezoidal diffusion gradients, where $\xi$ is the rise time of the trapezoidal ramp [86].

In Supplementary Appendix C, we provide a compact re-derivation of Lori’s approach, which found the following approximation:

$\frac{{〈e^{i\varphi }〉}_{\Delta ,\delta ,\mathbf{g}}}{{〈e^{i\varphi }〉}_{\Delta ,\delta ,\mathbf{g}=0}}\approx {\displaystyle \underset{{\mathbb{R}}^3}{\int }}P\left(\mathbf{r},t_{\exp }\right)e^{i{\mathbf{q}}^{\prime }\mathbf{r}}d\mathbf{r},(16)$

where ${\mathbf{q}}^{\prime }=\mathbf{q}\sqrt{t_{eff}/t_{\exp }}$ is a scaled q-space vector; $t_{eff}$ denotes the ‘effective’ diffusion time that appears in the b-value definition, i.e., $b=q^2{t}_{eff}$, which is equal to $t_{eff}=\Delta -\delta /3$ and $t_{eff}=\Delta -\delta /3+{\xi }^3/30{\delta }^2-{\xi }^2/6\delta$ for rectangular and trapezoidal diffusion gradients, respectively [86]. According to this result, the q-space formalism can still be employed to relate the diffusion propagator and the dMRI signal attenuation produced by a PGSE sequence with finite $\delta$. However, it must be corrected by evaluating the diffusion propagator at the total diffusion encoding time $t_{\exp }$ and using a modified q-space vector ${\mathbf{q}}^{\prime }$. Note that for narrow pulses, the correction converges to the classical q-space formalism with $t_{\exp }=t_{eff}=\Delta$, and ${〈\left.e^{i\varphi }\right|\mathbf{r}〉}_{\Delta ,\delta ,\mathbf{g}}=e^{i\mathbf{\text{qr}}}$, as expected.

The theoretical result in Equation 16 was confirmed in [85] by numerical simulations for homogeneous Gaussian diffusion, heterogeneous diffusion in permeable microscopic Gaussian domains, and diffusion inside restricted spherical reflecting domains. In all the analyses, this correction produced better results than using the original q-vector and the relationship $t_{\exp }=\Delta -\delta /3$, for rectangular diffusion gradients. It is important to notice that this approach may only provide a precise correction for displacement distributions that do not deviate significantly from a Gaussian distribution.

In this study, we will use this correction to evaluate our signal models in Equations 11, 14.

2.6 Spherical mean signals

The previous signal models, see Equations 11, 14, are based on the assumption of a single cylindrical surface. In the case of a distribution of cylinders with equal radius but multiple orientations, the orientation effect can be removed from Equation 11 by computing the orientation-averaged spherical mean signal $〈E〉$. Following the approach of [33, 87–89], we obtain,

$\begin{array}{l}\frac{〈E\left(q,\Delta ,\delta ,\xi ,a\right)〉}{E\left(q=0\right)}=\frac{1}{2}[{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{j=0}^k}\frac{{\left(-1\right)}^k}{{\left(k!\right)}^2}\left(\begin{array}{l}2k\\ k\end{array}\right){\left(\frac{aq}{2}\sqrt{\frac{t_{eff}}{t_{\exp }}}\right)}^{2k}\\ \times \left(\begin{array}{l}k\\ j\end{array}\right){\left(-1\right)}^j\frac{\mathrm{\Gamma }\left(j+\frac{1}{2}\right)-\mathrm{\Gamma }\left(j+\frac{1}{2},q^2{D}_{\parallel }{t}_{eff}\right)}{{\left(q^2{D}_{\parallel }{t}_{eff}\right)}^{j+1/2}}\\ +{\displaystyle \sum _{p=1}^{\infty }}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{j=0}^{p+k}}{e}^{-p^2\frac{D_{\parallel }{t}_{\exp }}{a^2}}\frac{{\left(-1\right)}^k}{k!\left(2p+k\right)!}\left(\begin{array}{c}2\left(p+k\right)\\ p+k\end{array}\right){\left(\frac{aq}{2}\sqrt{\frac{t_{eff}}{t_{\exp }}}\right)}^{2\left(p+k\right)}\\ \times \left(\begin{array}{l}p+k\\ j\end{array}\right){\left(-1\right)}^j\frac{\mathrm{\Gamma }\left(j+\frac{1}{2}\right)-\mathrm{\Gamma }\left(j+\frac{1}{2},q^2{D}_{\parallel }{t}_{eff}\right)}{{\left(q^2{D}_{\parallel }{t}_{eff}\right)}^{j+1/2}}],\end{array}(17)$

where $t_{eff}$ and $t_{\mathit{exp}}$ depend on the experimental parameters $\left\{\Delta ,\delta ,\xi \right\}$.

A detailed derivation of this expression is presented in Supplementary Appendix D, which also includes Lori’s q-space correction described in the previous section.

On the other hand, for the Gaussian diffusion model in Equation 14 with time-dependent radial diffusivity, the spherical mean signal is equivalent to that from an axis-symmetric diffusion tensor [33, 37, 43, 84, 87, 90]:

$\frac{〈E\left(q,\Delta ,\delta ,\xi ,a\right)〉}{E\left(q=0\right)}=\sqrt{\frac{\pi }{4}}{e}^{-b{D}_{\perp }^{app}}\frac{\mathit{erf}\left(\sqrt{b\left(D_{\parallel }-D_{\perp }^{app}\right)}\right)}{\sqrt{b\left(D_{\parallel }-D_{\perp }^{app}\right)}},(18)$

where $\mathit{erf}$ denotes the error function. In our model, the radial diffusivity $D_{\perp }^{app}\left(a,t=\Delta +\delta \right)$ depends on the cylinder radius $a$ and the total diffusion time according to the model defined in Equation 13 and incorporating Lori’s correction. Note that this correction does not affect the b-value since $b=q^2{t}_{eff}={q^{\prime }}^2{t}_{\exp }$ for rectangular and trapezoidal diffusion gradients.

2.7 Estimating the mean myelin sheath radius: what do we measure?

In this section, we will derive the spherical mean dMRI signal for a distribution of cylinders with different radii. Specifically, we will consider two cases:

1. Multiple concentric cylinders: This model represents the diffusion process of myelin water within a single axon. Each cylinder corresponds to a layer of the myelin sheath. The diffusion process is confined to these cylindrical surfaces, and the overall dMRI signal is the sum of contributions from each cylindrical layer; see Figure 1B.

2. Distribution of multiple concentric cylinders with different radii: This model represents a voxel with multiple axons, where the inner axon radius follows a Gamma distribution. The Gamma distribution is a flexible choice that can model a wide range of axon radius distributions observed in neural tissues [40]; see Figure 2.

We aim to define and estimate the ‘effective’ myelin sheath radius by approximating the signal from multiple cylindrical surfaces with the signal from a single cylindrical surface. The effective myelin sheath radius simplifies the complex distribution into a single representative value. This approach is analogous to axon diameter mapping techniques, which estimate an effective radius from an underlying distribution of inner axon radii [10, 11, 36–41, 45].

2.7.1 Effective myelin sheath radius for a single axon

The spherical mean dMRI signal $〈E_{axon}〉$ arising from N concentric cylindrical surfaces is given by

$\begin{array}{l}\frac{〈E_{axon}\left(q,\Delta ,\delta ,\xi \right)〉}{E_{axon}\left(q=0\right)}=\frac{1}{E_{axon}\left(q=0\right)}{\displaystyle \sum _{i=1}^N}〈E\left(q,\Delta ,\delta ,\xi ,a_i\right)〉\\ {\displaystyle \sum _{i=1}^N}\frac{E\left(q=0,a_i\right)}{E_{axon}\left(q=0\right)}〈S\left(q,\Delta ,\delta ,\xi ,a_i\right)〉,\\ ={\displaystyle \sum _{i=1}^N}\frac{a_i}{\left({\displaystyle \sum _{j=1}^N}{a}_j\right)}〈S\left(q,\Delta ,\delta ,\xi ,a_i\right)〉.\end{array}(19)$

where the summation is over all cylinder’s radii from the inner radius $a_1$ to the outer radius $a_N$, and $〈S\left(q\right)〉=〈E\left(q\right)〉/E\left(q=0\right)$ denotes the spherical mean dMRI signal produced by each cylinder normalized by its baseline signal without diffusion weighting (q = 0 image); see Equations 17, 18. The term $E_{axon}\left(q=0\right)$ was included on both sides of the equation on purpose. Since $E\left(q=0\right)$ is proportional to the number of diffusing spin particles, $E\left(q=0,a_i\right)/E_{axon}\left(q=0\right)$ is the ratio of the number of those particles on the cylinder with the radius $a_i$ and the total number on all cylinders. Assuming the same proton density (i.e., number of particles per unit surface area) and cylinder length, this ratio is the surface area of the i-th cylinder divided by the total surface area of all cylinders, or equivalently, the radius of the i-th cylinder divided by the sum of all radii.

We can substitute the normalized spherical mean signal obtained for the general model (Equation 17) or the Gaussian approximation with time-dependent radial diffusivity (Equation 18) in Equation 19. When the resulting signal is approximated by the signal from a single cylindrical surface, then

$\frac{〈E\left(q,\Delta ,\delta ,\xi ,a_{eff}\right)〉}{E\left(q=0\right)}\approx {\displaystyle \sum _{i=1}^N}\frac{a_i}{\left({\displaystyle \sum _{j=1}^N}{a}_j\right)}〈S\left(q,\Delta ,\delta ,\xi ,a_i\right)〉,(20)$

where $a_{eff}$ is the effective radius. This effective radius defined by Equation 20 represents the MRI-visible radius that considers that the measured signal is weighted by the radius, such that the outer cylinder contributes more than the inner cylinder to the measured data. Assuming that all cylinders have the same distance between them, then $a_{eff}$ will be more biased towards $a_N$ than towards $a_1$ from the arithmetic mean $a_{eff}\ge 〈a〉$, defined by

$\begin{array}{l}〈a〉=\frac{1}{N}{\displaystyle \sum _{i=1}^N}{a}_i,\\ =\frac{1}{N}{\displaystyle \sum _{i=1}^N}\left(a_1+\left(i-1\right)\Delta a\right),\\ =a_1+\left(\frac{N-1}{2}\right)\Delta a,\\ =\frac{a_i+a_o}{2},\end{array}(21)$

where $\Delta a=a_{i+1}-a_i$ is the distance between two consecutive cylinders, and the outer cylinder’s radius is $a_N=a_1+\left(N-1\right)\Delta a$. In Equation 21, we replaced $a_1$ and $a_N$ with the inner and outer axon radii, $a_i$, $a_o$, respectively.

2.7.2 Effective myelin sheath radius for a distribution of axon radii

For a sample of myelinated axons with the same g-ratio, $g=a_i/a_o$, and the distribution of inner axon radius parameterized by $P\left(a_i\right)$, the marginal distribution of myelin sheath (cylinder) radii is given by

$P\left(a\right)=\eta {\displaystyle \underset{0}{\overset{\infty }{\int }}}U\left(a\left|a_i,\right.a_o\right)P\left(a_i\right)d{a}_i,(22)$

where $\eta$ is the normalization constant ensuring that ${\int }_0^{\infty }P\left(a\right)da=1$, and $U\left(a\left|a_i,\right.a_o\right)$ is a uniform distribution modeling the myelin layers of each axon as uniformly distributed cylinders in the interval $\left[a_i,a_o\right]$,

$\begin{array}{l}U\left(a\left|a_i,\right.a_o\right)=\frac{1}{a_o-a_i}{\mathbf{1}}_{\left[a_i,a_o\right]}\left(a\right),\\ =\frac{g}{a_i\left(1-g\right)}{\mathbf{1}}_{\left[a_i,a_i/g\right]}\left(a\right),\end{array}(23)$

which is written in terms of $a_i$ and $g$. The indicator function ${\mathbf{1}}_{\left[a_i,a_o\right]}\left(a\right)$ is equal to 1 if $a_i\le a\le a_o$ and 0 otherwise.

We assume a Gamma distribution for the inner radius as in [40]:

$P\left(a_i\right)=\frac{{\kappa }^{\mu }}{\mathrm{\Gamma }\left(\mu \right)}{a}_i^{\mu -1}{e}^{-\kappa a_i},\text{for}{a}_i>0,\mu ,\kappa >0,(24)$

where $\mathrm{\Gamma }\left(\mu \right)$ is the Gamma function, and $\mu$ and $\kappa$ are the shape and inverse scale parameters, respectively, such that the mean radius and variance are $〈a〉=\mu /\kappa$ and ${\sigma }^2=\mu /{\kappa }^2$.

Inserting Equations 23, 24 into Equation 22, and considering that at a given radius $a$ only those cylinders in the range from $\left[a\cdot g,a/g\right]$ contribute to the integral (i.e., the population of cylinders from axons with inner and outer radii ranging from [$a_i=a\cdot g$, $a_o=a$] to [$a_i=a$, $a_o=a/g$]), we obtain

$\begin{array}{l}P\left(a\right)=\eta \frac{{\kappa }^{\mu }}{\mathrm{\Gamma }\left(\mu \right)}\frac{g}{\left(1-g\right)}{\displaystyle \underset{a\cdot g}{\overset{a/g}{\int }}}{a}_i^{\mu -2}{e}^{-\kappa a_i}d{a}_i,\\ =\frac{\kappa }{\mathrm{\Gamma }\left(\mu \right)}\frac{g}{\left(1-g^2\right)}\left[\mathrm{\Gamma }\left(\mu -1,a\cdot g\cdot \kappa \right)-\mathrm{\Gamma }\left(\mu -1,\frac{a\cdot \kappa }{g}\right)\right],\end{array}(25)$

where $\mathrm{\Gamma }\left(s,x\right)$ denotes the upper incomplete Gamma function. The complete derivation is developed in Supplementary Appendix E. Note that for axons with a very small number of myelin layers, $g\to 1$ and $P\left(a\right)\approx P\left(a_i\right)$. Figure 2 shows an example of a distribution of inner axon radius sampled from the splenium of the corpus callosum of a human brain reported by [45, 91] and the corresponding marginal distribution of myelin sheath radii assuming $g=0.6$.

Figure 2

Figure 2. Distribution of radius. Left Panel: The diagram illustrates a population of axons within a voxel, displaying varying inner radii while maintaining a constant g-ratio. Right Panel: This graph presents the distribution of inner axon radii sampled from the splenium of the Corpus Callosum of an ex-vivo human brain (data from [91]). The Gamma distribution fitting the measured inner radii is depicted in blue, and the corresponding marginal distribution of the myelin sheath radius calculated using Equation 25 and assuming a constant g-ratio of 0.6, is shown in yellow-orange. The Gamma distribution was fitted to the data using a Maximum Likelihood approach, as implemented in the gamfit function in @Matlab. This visualization highlights the relationship between the inner axon radius distribution (mean = 0.68 µm, variance = 0.11 µm²) and the myelin sheath radius distribution (mean = 0.77 µm, variance = 0.195 µm²).

The spherical mean dMRI signal produced by such a distribution of cylinders is

$\frac{〈E_{dist}\left(q,\Delta ,\delta ,\xi \right)〉}{E_{dist}\left(q=0\right)}=\frac{{\displaystyle \underset{0}{\overset{\infty }{\int }}}aP\left(a\right)〈S\left(q,\Delta ,\delta ,\xi ,a\right)〉da}{{\displaystyle \underset{0}{\overset{\infty }{\int }}}aP\left(a\right)da},(26)$

If the distribution of the inner radius $P\left(a_i\right)$ and the g-ratio are known from histological measurements, we can estimate $P\left(a\right)$ from Equation 25. The dMRI signal in Equation 26 can be computed numerically using Equation 17 or Equation 18 for a given set of PGSE acquisition parameters, and the effective radius $a_{eff}$ can then be estimated by fitting the single-cylinder model to the resulting signal.

Following the approach described by [41], the effective radius can be approximated by the weighted-mean radius

$\begin{array}{l}{a}_{eff}\approx \frac{{\displaystyle \underset{0}{\overset{\infty }{\int }}}aP\left(a\right)N\left(a\right)da}{{\displaystyle \underset{0}{\overset{\infty }{\int }}}P\left(a\right)N\left(a\right)da},\\ =\frac{{\displaystyle \underset{0}{\overset{\infty }{\int }}}{a}^{2}P\left(a\right)da}{{\displaystyle \underset{0}{\overset{\infty }{\int }}}aP\left(a\right)da}\end{array}(27)$

where $N\left(a\right)$ is the number of diffusing particles as a function of the radius $a$. In our case, $N\left(a\right)$ is proportional to the surface area of the cylinder and, therefore, to its radius. Consequently, the signal contribution from each cylinder is approximately proportional to its radius. Thus, we expect $a_{eff}$ to correlate with the ratio $〈a^2〉/〈a〉$ determined from the underlying distribution $P\left(a\right)$.

Alternatively, another approximation can be obtained by following the approach presented by [37] using the Gaussian approximation with time-dependent radial diffusivity. When assuming small myelin sheath radii such that $D_{\parallel }t\gg a^2$ and $D_{\parallel }\gg D_{\perp }^{app}$, for low and moderate b-values, the normalized spherical mean dMRI signal can be approximated by:

$〈S\left(q,\Delta ,\delta ,\xi ,a\right)〉\approx \sqrt{\frac{\pi }{4}}\frac{\mathit{erf}\left(\sqrt{b{D}_{\parallel }}\right)}{\sqrt{b{D}_{\parallel }}}\left(1-b\frac{a^2}{2t_{\exp }}\right),(28)$

where we used Equations 13, 18 and Lori’s correction. Inserting this equation into the right-hand side of Equation 26 and equating this expression to the signal arising from a single cylindrical surface with radius $a_{eff}$ we obtain:

$\begin{array}{l}〈S\left(q,\Delta ,\delta ,\xi ,a_{eff}\right)〉\approx \sqrt{\frac{\pi }{4}}\frac{\mathit{erf}\left(\sqrt{b{D}_{\parallel }}\right)}{\sqrt{b{D}_{\parallel }}}\frac{{\displaystyle \underset{0}{\overset{\infty }{\int }}}aP\left(a\right)\left(1-b\frac{a^2}{2t_{\exp }}\right)da}{{\displaystyle \underset{0}{\overset{\infty }{\int }}}aP\left(a\right)da},\\ =\sqrt{\frac{\pi }{4}}\frac{\mathit{erf}\left(\sqrt{b{D}_{\parallel }}\right)}{\sqrt{b{D}_{\parallel }}}\left(1-\frac{b}{2t_{\exp }}\frac{{\displaystyle \underset{0}{\overset{\infty }{\int }}}{a}^{3}P\left(a\right)da}{{\displaystyle \underset{0}{\overset{\infty }{\int }}}aP\left(a\right)da}\right).\end{array}(29)$

Comparing Equations 29, 28 we obtain

$a_{eff}^2=\frac{{\displaystyle \underset{0}{\overset{\infty }{\int }}}{a}^{3}P\left(a\right)da}{{\displaystyle \underset{0}{\overset{\infty }{\int }}}aP\left(a\right)da},(30)$

Thus, we might also expect $a_{eff}$ to correlate with the expression ${\left(〈a^3〉/〈a〉\right)}^{1/2}$.

In the Results section, we will compare these two effective radius definitions with the numerical effective radius determined by fitting Equation 26 to the theoretical model corresponding to a single cylinder. This evaluation will use histological measurements of inner axon radii sampled from four regions of the Corpus Callosum in a human brain [91], which will be converted into distributions of myelin sheath radii according to Equation 25.

3 Methods

3.1 Monte Carlo simulations

Monte Carlo Diffusion Simulations (MCDS) were employed as a benchmark to validate the proposed models. We used an MC simulator developed by our group, available at https://github.com/jonhrafe/Robust-Monte-Carlo-Simulations [92]. This simulator has been validated against analytical models across multiple geometries, including impermeable planes, cylinders, and spheres [92]. For this study, we extended its capabilities to incorporate new myelin water diffusion models, implementing two geometrical structures: 3D infinite, impermeable cylinders and spiral surfaces.

The analytical models were validated by comparing their predicted dMRI signals to those generated by the MC simulations for identical impermeable cylindrical surfaces. Additionally, the dMRI signals from concentric cylinders were compared with those from spiral surfaces to assess the assumption presented in Section 2.1 (Figure 1). This assumption suggests that net radial displacements along the spiral trajectory are negligible, which allows the diffusion process in the more complex spiral geometry to be approximated as that in concentric cylinders.

3.2 Geometrical models

3.2.1 Cylindrical surfaces

We simulated diffusion on infinite, impermeable cylindrical surfaces. The diffusion process was simulated using a fixed step size along both the z-axis (aligned with the main axis of the cylinder) and the curved trajectory in the x-y plane, given by $l=\sqrt{2D_{\parallel }t/N_t}$​, where $N_t$ is the number of Monte Carlo steps and $t$ is the total diffusion time. At each step, the particle’s z-coordinate was updated as $z\leftarrow z\pm l$, with the direction randomly selected to simulate upward and downward motion. In the x-y plane, the angular displacement was selected to maintain a constant arc length $l$, i.e., $\theta \leftarrow \theta \pm l/a$, allowing particles to move in either rotational direction. The radius $a$ was constant, reflecting the cylindrical surface’s geometry.

For each b-value, dMRI signals were generated from 50 independent cylinders with radii uniformly spaced from 0.1 µm to 5.0 µm in increments of 0.1 µm. To simulate the myelin water dMRI signal from a single axon with specific inner and outer radii, we calculated the radius-weighted sum of the signals from all cylindrical surfaces in this range, following Equation 19.

To replicate the myelin water dMRI signal based on voxelwise realistic distributions of myelin radii, we performed the following steps:

1. Converted histological distributions of inner axon radii from [91] into myelin sheath radii using Equation 25, assuming a constant g-ratio of 0.7.

2. Computed the spherical mean dMRI signal for each resulting distribution by evaluating the integral in Equation 26, discretized using the same grid of 50 radii ranging from 0.1 to 5.0 µm as used in the MC simulations.

3.2.2 Spiral surfaces

The diffusion process was similarly simulated for the spiral surfaces using a fixed step size $l$ along the z-axis and the x-y plane. The curved trajectory in the x-y plane was determined by the particle’s position on the spiral. The radius $a\left(\theta \right)$ of the spiral varies with the polar angle $\theta$ in the x-y plane, according to $a\left(\theta \right)=a_i+\left(d_s/2\pi \right)\theta$, where $a_i$ is the inner radius and $d_s$ is the distance between successive layers of the spiral. The inter-layer distance was fixed to $d_s$ = $d_m+d_w$ = 7.5nm, based on histological data reported by [73]. In this context, $d_m$ and $d_w$ represent the thickness of the myelin layer and the spacing filled by myelin water, respectively. Therefore, $d_s$ corresponds to the distance between the centers of the gaps filled by myelin water in an axon. The polar angle $\theta$ ranged from 0 to the maximum value for which $\left(\theta \right)=a_o$.

To assess whether the dMRI signals from water molecules confined to spiral surfaces can be approximated by those from concentric cylindrical surfaces, we generated spiral geometries with g-ratios of 0.6, 0.7, and 0.8, consistent with values reported in histological studies [93, 94]. Since the results across different g-ratios were comparable, we present findings for g-ratio = 0.7, using three geometries with inner and outer radii of 0.5/0.7 µm, 0.7/1.0 µm, and 1.0/1.4 µm, respectively.

The resulting signals were compared to those from corresponding cylindrical surfaces using the same PGSE sequence parameters described in the next section.

3.3 Simulation protocol

The diffusion process was simulated for both geometrical models using a total diffusion time of $t$ = 20 ms and $N_t$ = 15,000 steps per particle. We conducted a bootstrap-based analysis to ensure convergence of the simulations, as outlined in [92]. A total of 75,000 particles were uniformly distributed on each cylindrical or spiral surface. Three values of parallel diffusivity (D_∥ = 0.3, 0.5, 0.8 μm²/ms) were used to cover the range of myelin water diffusivities reported by [65].

A PGSE sequence with trapezoidal diffusion gradients was used to generate dMRI signals. The sequence was based on the specifications of a Connectome 2.0 scanner, employing a maximum gradient strength of G = 500 mT/m and a maximum slew rate of SR = 600 T/m/s [68], yielding to a trapezoidal ramp rise time $\xi$ = G/SR = 0.833 ms. The protocol included 90° and 180° pulse durations of 2 ms and 4 ms, respectively. Six b-values were selected using the shortest possible TE for each case while maintaining maximum G and SR, following the implementation described in [70, 71]. Table 1 details the experimental parameters.

Table 1

Table 1. Experimental parameters for Monte Carlo simulations using a PGSE sequence with trapezoidal diffusion gradients.

For each b-value, dMRI signals were generated for 92 gradient orientations uniformly distributed on the unit sphere, along with the signal for b = 0. The subsequent analyses focused on the spherical mean signal normalized by the b = 0 signal.

4 Results

4.1 Diffusion diffraction pattern: single cylinder

Figure 3 illustrates the theoretical spherical mean dMRI signal from a cylindrical surface, as generated by the general model presented in Equation 17 using a PGSE sequence with trapezoidal diffusion gradients. The signal is shown for b-values ranging from 0 to 100 ms/μm² and for three cylinders with radii of 0.3 µm, 1.0 µm, and 3.0 µm.

Figure 3

Figure 3. Theoretical spherical mean signal attenuation for cylindrical surfaces. The signal was generated using the general model presented in Equation 17 for b-values ranging from 0 to 100 ms/μm², with diffusion time parameters of $\Delta =9.446$ ms, $\delta =4.612$ ms, and $D_{\parallel }=0.8$ µm²/ms. The signal decay for b-values from 0 to 3 ms/μm² is displayed in a separate zoomed-in axis, as indicated by the rectangular box. The signal attenuation is plotted on a logarithmic scale for three cylinders with radii of 0.3 µm (blue), 1.0 µm (green), and 3.0 µm (orange) as a function of the b-value.

For relatively low b-values (approximately below 3 ms/μm²), the logarithm of the signal approximates a linear relationship. This linearity suggests that a Gaussian model could be valid in this regime. However, as the b-value increases, deviations from Gaussianity become apparent, and signal oscillations, known as diffraction patterns, emerge. These diffraction-like patterns have been reported in other geometries where diffusion is confined, such as planar, cylindrical, and spherical domains [95–97].

4.2 Single cylinder dMRI signal using “realistic” acquisition parameters vs. MC simulations

To assess the accuracy of the new analytical models proposed in this study, we compared the predicted dMRI signals with those generated by MC simulations. Figure 4 shows the theoretical spherical mean dMRI signals from cylindrical surfaces as a function of the radius, as predicted by both the general analytical model and the Gaussian approximation with time-dependent radial diffusivity (Equations 17, 18, respectively) using a PGSE sequence with trapezoidal diffusion gradients. Additionally, the figure includes the dMRI signals obtained from the MC simulations for validation purposes. This comparison was conducted over a range of parallel diffusivities (D_∥ = 0.3, 0.5, 0.8 μm²/ms) and practical b-values from 0.8 to 3.0 ms/μm², achievable in preclinical and human scanners equipped with strong diffusion gradients.

Figure 4

Figure 4. Sensitivity of the spherical mean dMRI signal as a function of myelin sheath radii for different diffusivities. The signals were generated using the general model (Equation 17, continuous lines), the Gaussian approximation (Equation 18), dashed lines), and Monte Carlo (MC) numerical simulations (dots) for the following b-values: [0.8, 1.0, 1.5, 2.0, 2.5, 3.0] ms/µm², using a PGSE sequence with parameters listed in Table 1. (A–C) show results corresponding to diffusivities of $D_{\parallel }=0.3$, $D_{\parallel }=0.5$, and $D_{\parallel }=0.8$ µm²/ms, respectively. The normalized signal amplitudes from the analytical models are displayed for myelin sheath radii ranging from 0 to 5 μm, and the MC signals were generated for 50 discrete radii ranging from 0.1 to 5 µm.

Increasing the b-value results in greater attenuation of the dMRI signal as a function of the radius across all three diffusivity values. At a b-value of 3.0 ms/μm², the signal exhibits maximum sensitivity to myelin sheath radii in the 0.5–3.0 µm range. However, at this higher b-value, we observe the largest, albeit still minor, deviations between the signals predicted by the analytical models and those generated by the MC simulations. Notably, the agreement between the models and simulations is strongest for the lowest diffusivity (D_∥ = 0.3 μm²/ms, panel A). It diminishes as diffusivity increases, with the largest discrepancy observed at D_∥ = 0.8 μm²/ms (panel C).

For this acquisition protocol, the signal shows minimal sensitivity to myelin radii smaller than 0.5 µm and larger than 3.5–4.0 µm. This result indicates that the method is best suited for detecting myelin sheath sizes in the 0.5–3.5 µm range. Across all b-values, the Gaussian approximation closely follows the analytical model, particularly for radii below 4.0 µm, further confirming the accuracy of the approximation in this parameter range.

4.3 Spiral surfaces vs. concentric cylinders: MC simulations and analytical models

The results from the experiment comparing the spherical mean dMRI signals generated by MC simulations for spiral geometries and multiple concentric cylinders are presented in Figure 5. Specifically, Figure 5 shows the dMRI signals as a function of the six b-values employed. The signal from a spiral geometry with inner and outer radii of 0.7 µm and 1.0 µm is compared with the radius-weighted signal from multiple concentric cylinders within the same radius range, calculated using Equation 19. Additionally, we display the signals from individual cylindrical surfaces with radii ranging between 0.7 µm and 1.0 µm, obtained from both MC simulations and the analytical models. Panels A and B correspond to results for diffusivities of D_∥ = 0.3 μm²/ms and D_∥ = 0.8 μm²/ms, respectively.

Figure 5

Figure 5. Comparison of dMRI signals from spiral surfaces and concentric cylinders. Spherical mean dMRI signals as a function of six employed b-values and results from Monte Carlo (MC) simulations for spiral geometries and multiple concentric cylinders. The signals are generated for a spiral with inner and outer radii of 0.7 µm and 1.0 µm, respectively, alongside radius-weighted signals from concentric cylinders within the same radius range. Signals from individual cylindrical surfaces with radii between 0.7 µm and 1.0 µm are plotted using both MC simulations and analytical models. (A, B) show results for $D_{\parallel }=0.3$ µm²/ms and $D_{\parallel }=0.8$ µm²/ms, respectively. In (B), we highlight a region where the most significant discrepancies were observed between the signals computed using the analytical models and those obtained from the MC simulations.

For both diffusivity values, we observe a strong agreement between the MC-generated signals for the spiral geometry and the radius-weighted aggregation of signals from concentric cylinders with the same range of radii. This result suggests that the spiral geometry can be accurately approximated by multiple concentric cylinders. Notably, for the lower diffusivity (D_∥ = 0.3 μm²/ms, panel A), the analytical model’s predictions for individual cylinders closely match the signals generated by MC simulations. Furthermore, the signal produced by the spiral geometry is very similar to that of a single cylinder with a radius intermediate to the inner and outer radii. This implies that when fitting these signals with a single-radius model, the estimated effective radius would likely correspond to a value close to the average radius of the spiral.

However, for simulations at the higher diffusivity (D_∥ = 0.8 μm²/ms, panel B), the signal decay predicted by the analytical models as a function of the b-value is more pronounced than the decay observed in the MC simulations. This result indicates potential inaccuracies in the analytical model at higher diffusivities and larger b-values. Consequently, the effective radius predicted by the analytical models will likely be biased towards a smaller value than the actual radius.

The results for spirals with other inner and outer radii were consistent with these findings. Specifically, the observed discrepancy for D_∥ = 0.8 μm²/ms was reduced for the spiral with a larger inner radius of 1.0 µm. Conversely, the disagreement increased for the smaller spiral with an inner radius of 0.5 µm (results not shown).

4.4 Effective radius from histological measurements for distributions of cylinders

Figure 6 compares the effective radii estimated from simulated dMRI data against three different metrics derived from the distribution of myelin sheath radii in four regions of interest within the Corpus Callosum: axons connecting the prefrontal, motor, parietal, and visual cortices. The inner axon radii for these regions, as reported by [91], were modeled using Gamma distributions. These distributions were subsequently transformed into myelin sheath radii distributions using Equation 25 and a constant g-ratio of 0.7.

Figure 6

Figure 6. Distributions of inner axon and myelin sheath radii and estimated effective radius. Four subplots are presented, each corresponding to a different region of interest in the Corpus Callosum of a human brain. Each subplot includes a histogram of the measured inner axon radius (data from [91]), along with the best-fitting Gamma distribution (in blue) and the derived myelin sheath radius distribution estimated using Equation 25 (in yellow-orange). The effective radius $a_{eff}$, estimated as the radius from the single-cylinder model (see Equation 17) that best fits the signal generated from the whole distribution of myelin sheath radius (see Equation 26), is plotted, along with three representative metrics of the distribution, including the mean value $〈a〉$ and the second- and third-moment based metrics $〈a^2〉/〈a〉$ and ${\left(〈a^3〉/〈a〉\right)}^{1/2}$ derived in Equations 27, 30, respectively. These results correspond to simulations using D_∥ = 0.5 μm²/ms.

We then generated the spherical mean dMRI signals corresponding to these distributions by discretizing Equation 26 and employing the MC simulated signals. We assumed a parallel diffusivity of D_∥ = 0.5 μm²/ms. The generated signals were fitted to the general single-cylinder model in Equation 17 to estimate the effective radius. Figure 6 presents the effective radii $a_{eff}$, the mean radii $〈a〉$ obtained from the distributions, and the second- and third-moment-based radii metrics $〈a^2〉/〈a〉$ and ${\left(〈a^3〉/〈a〉\right)}^{1/2}$​, as defined in Equations 27, 30.

The myelin sheath radii distributions in Figure 6 exhibit slightly longer right-hand tails and lower frequency values for small radii compared to the inner axon radii distributions, as expected. This difference arises because the myelin sheath radii represent all possible layer radii within the range defined by the inner and outer radii for all axons. Hence, it includes contributions from myelin layers near the inner and outer boundaries. These two distributions converge further as the g-ratio increases, as described by Equation 25. This trend is noticeable when comparing the distributions in Figure 2 for a g-ratio of 0.6 with those in Figure 6 employing a g-ratio of 0.7.

The results show that for the distributions with smaller radii (Prefrontal and Parietal regions), the estimated effective radius $a_{eff}$ closely matches the mean radius $〈a〉$. However, for the Motor and Visual regions, with larger radii distributions, the effective radius aligns more closely with the second-moment-based metric $〈a^2〉/〈a〉$, followed by the third-moment-based metric ${\left(〈a^3〉/〈a〉\right)}^{1/2}$. These findings suggest that the appropriate descriptor of the distribution may depend on the range of radii in each region.

To further investigate the relationships between the effective radius and the derived metrics from the myelin sheath radii distributions, we present a correlation analysis in Figure 7. This figure illustrates the correlations between the effective radius and the three descriptive metrics across experiments conducted with three distinct diffusivities.

Figure 7

Figure 7. Correlation between the effective radius and descriptive metrics. This figure shows the correlations between the effective radius $a_{eff}$ (y-axis), estimated from dMRI signals, and three derived metrics (x-axis) from myelin sheath radii distributions: mean radius $〈a〉$ (blue), second-moment-based radius $〈a^2〉/〈a〉$ (orange), and third-moment-based radius ${\left(〈a^3〉/〈a〉\right)}^{1/2}$ (green). Panels A, B, and C depict the results for three distinct diffusivity values: D_∥ = 0.3, 0.5, and 0.8 μm²/ms. The Pearson’s Correlation Coefficient (PCC) and the corresponding p-value are reported for each analysis. Each set of points represents the values estimated from the four distributions shown in Figure 6.

As shown in Figure 7, although these metrics reflect different aspects of the myelin sheath radii distributions, they exhibit significant correlations with the effective radius. Notably, the second-moment-based radius $〈a^2〉/〈a〉$ demonstrated the strongest linear correlation (and smallest p-value) with $a_{eff}$​ across all diffusivity values, indicating its potential as a reliable descriptor of effective radii. The third-moment-based radius ${\left(〈a^3〉/〈a〉\right)}^{1/2}$ closely followed this trend, while the average radius showed less strong correlations. Interestingly, the analysis reveals a trend where the estimated effective radius tends to decrease with increasing diffusivity, particularly pronounced in distributions characterized by smaller axon radii.

5 Discussion

In this proof-of-concept study, we developed two models for the dMRI signal arising due to water molecular displacements on cylindrical surfaces. We focused on potential applications for modeling the dMRI signal associated with myelin water in brain tissues. In the first, more general model, we derived an exact analytical expression for the dMRI signal using the diffusion propagator formalism based on the narrow pulse approximation. The second model employs a Gaussian approximation with time-dependent radial diffusivity, offering a simpler analytical relationship. We also developed approximate signal expressions for PGSE protocols with trapezoidal and rectangular diffusion gradients, extending beyond the narrow pulse assumption.

We derived the spherical mean signal expressions for both models, which are theoretically independent of axonal orientation effects. The spherical mean signal remains invariant to orientation dispersion, as it is approximately equivalent to whether the axons within a voxel have varying orientations or are aligned to the same orientation [33, 98]. While it is theoretically feasible to estimate both fiber orientations and the effective radius of the myelin sheath, such fitting procedures may be unstable. To address this challenge, we adopted a strategy inspired by previous studies on axon diameter mapping. These studies also employ the spherical mean approach to minimize the influence of orientation effects [37, 43, 46, 99], a well-known confounding factor that can bias axon diameter estimates. Indeed, when the dispersion is not accurately incorporated into the model, it could alter the estimated radial and parallel myelin water diffusivities. Conversely, when the spherical mean signal is used, the estimated diffusivities more accurately reflect the intrinsic diffusivities of myelin water.

We also derived expressions for the dMRI signal from multiple concentric cylinders as the radius-weighted sum of signals. This was done to account for the dependence of signal intensities on the cylinders’ surface areas and, thus, their radii. We further generalized this approach to consider a distribution of myelin sheath radii. Various approximations were introduced to enhance our understanding of the effective radius—the radius estimated by fitting the signal from a radius distribution to a single-radius model. Finally, we extended our MC diffusion simulation toolbox to simulate the diffusion process confined on cylindrical and spiral surfaces to compare the analytical and numerical dMRI signals.

Validating the proposed models would require comparing the effective radii estimated from dMRI data and the corresponding values measured from histology on the same brain regions. However, since histological studies typically report only the inner radius distribution, we introduced a new analytical approach to convert this distribution into a distribution of myelin sheath radii based on the assumption of a constant g-ratio across all axons in the sample. It is important to emphasize that this analytical relationship is primarily a practical tool for leveraging existing histological data. If new histological studies provide direct measurements of myelin sheath radii, we would no longer need to rely on this approximation for validation.

The proposed models can potentially estimate the effective myelin sheath radius from real dMRI data. For example, our models could be directly applied in diffusion-T1 experiments using inversion recovery sequences that effectively isolate signals from myelin water, as outlined in [65]. Similarly, for acquisition sequences where signals from other compartments are not entirely suppressed—such as in diffusion-T2 hybrid sequences proposed by [63, 64] or the magnetization-prepared dMRI sequence described by [66]—our models could be integrated with existing multi-compartment dMRI frameworks, e.g. [31, 33, 100, 101], to concurrently fit the myelin water component along with parameters for other compartments. Additional investigations are needed to identify the optimal acquisition protocols for these multi-compartment fittings, focused on mitigating model fitting degeneracies [102]. These approaches could be applied to both ex vivo and in vivo data using scanners with strong diffusion gradients, leveraging recent advances [70–72] that enhance the myelin water dMRI signal by reducing echo times.

Our MC simulations employed parallel diffusivity values as reported by [65], specifically D_∥ = 0.37 μm²/ms in excised frog sciatic nerve for the double-inversion-recovery sequence. Since their experiments were conducted within 1 hour post-euthanasia and lasted approximately 90 min, this relatively short post-mortem interval likely helped preserve some of the tissue’s original diffusion properties compared to in vivo studies, thereby minimizing significant alterations due to dehydration or tissue degradation. However, the reduced temperature (20°C) relative to the typical in vivo temperature (around 37°C) may have contributed to decreased diffusivity. Hence, we expect the diffusivity values they reported to be lower than those observed in vivo. On the other hand, we anticipate that myelin water diffusivity will be lower than in other WM compartments due to its higher bound water content, which results in shorter relaxation times and reduced mobility. Therefore, we employed myelin water parallel diffusivities in the 0.3–0.8 μm²/ms range.

This study is not the first to simulate the dMRI signal from myelin water. To our knowledge, two previous works have specifically addressed the multi-wrapping nature of myelin [103, 104]. In the first study [103], this aspect was modeled implicitly by assuming a higher myelin water diffusivity in the tangential direction than the radial one. MC simulations were employed to assess the sensitivity of dMRI models to the diffusive properties of myelin water. Their findings indicate that myelin water could influence the apparent diffusion coefficient and kurtosis measured transverse to the orientation of WM tracts. In contrast, the second study [104] conducted MC simulations to examine water exchange through myelin sheaths by explicitly creating a spiraling myelin structure. They observed sub-second exchange times for thin axons with fewer wraps, highlighting the importance of modeling water exchange across WM compartments, especially in clinical studies on demyelinating diseases and the developing infant brain. Conversely, a slow exchange rate was observed in axons with more than eight myelin sheaths, typical of healthy WM in humans, supporting the assumption of impermeable membranes.

While other methods exist for quantifying WM microstructure parameters, including the inner axon radius and myelin content, each has inherent limitations. Myelin volume, often combined with the fiber volume fraction estimated from dMRI data to calculate the mean g-ratio, is typically determined using Magnetization Transfer (MT) or Multi-echo T2 (MET2) relaxometry techniques. However, although MT and MET2 techniques are known for their sensitivity to changes in myelin content, they are not exclusively specific to myelin, as other tissue compartments can also influence the measured signal [53, 105, 106]. Similarly, inner axon radius mapping techniques based on dMRI data face a resolution limit below which the radii of smaller axons cannot be reliably estimated [39, 74]. As such, the estimated effective inner radius typically represents the right-hand tail of the inner axon radius distribution rather than the entire distribution [37]. As myelin imaging techniques (i.e., MT and MET2) are not affected by the same resolution limit, care should be taken when combining estimates from these techniques to predict total myelin thickness (i.e., the difference between the outer and inner axon radii).

To the best of our knowledge, we present the first models for estimating myelin sheath radii exclusively using dMRI data, offering a novel imaging biomarker for detecting changes in myelin thickness. Although the method does not directly estimate the distance between the inner and outer layers of the myelin, it provides an integrated measure representing the entire distribution of myelin layer radii. The effective myelin sheath radius is derived by fitting a single-cylinder-surface model to the dMRI signal. In a hypothetical sample of identical axons with the same g-ratio, the effective radius closely approximates the mean of the inner and outer axon radii. In more realistic scenarios, where axon radii vary, and each axon has a distinct g-ratio, it reflects a population-weighted average with larger myelin layers contributing more substantially to the overall value.

Although our results are promising, several limitations need to be addressed in future work:

i. While the analytical models closely match MC simulations under various experimental conditions, discrepancies emerge at high b-values and large diffusivities. These inaccuracies arise from the approximations introduced to facilitate modeling. We initially derived our models using the narrow pulse approximation and later applied a correction framework to extend their applicability beyond this scheme. However, it is important to note that this correction framework primarily provides a valid approximation for Gaussian diffusion. The diffusion process deviates from Gaussian behavior in scenarios involving small cylinder radii, high diffusivities, and high b-values. One potential approach to address this limitation is to adapt the multiple propagator method introduced by [107] and refined by [108] to our specific models. Additionally, exploring a data-fitting approach based on a dictionary of precomputed MC signals may allow us to circumvent the limitations imposed by the theoretical approximations.

ii. The myelin sheath radius estimations are constrained by a resolution limit, influenced by both the strength of the diffusion gradient and the signal-to-noise ratio (SNR). For the employed acquisition parameters (i.e., G_max = 500 mT/m), our results indicate that signals for myelin sheath radii smaller than 0.5 μm and higher than 3.5 μm are indistinguishable (Figure 4). This range shifts with the diffusion gradient strength: weaker gradients make it harder to detect smaller myelin sheaths, whereas stronger gradients, like those in preclinical scanners (e.g., G_max = 1,500 mT/m), improve sensitivity to thinner layers. We did not conduct a formal resolution analysis akin to [39, 74] for estimating inner axon diameters, which would involve determining the exact resolution limit and its dependence on G_max and SNR. However, combining measurements acquired with different diffusion gradient strengths could extend the sensitivity range, although this approach is more feasible in preclinical settings where stronger diffusion gradients are available. In practice, the myelin water dMRI signal attenuation is primarily influenced by myelin layers with radii within the detectable range, with greater sensitivity to the right-hand tail of the radii distribution. Therefore, clinical applications should target pathologies involving larger axons, as smaller myelin layers may fall below the resolution limit. This limitation is not unique to our method. Similar constraints affect other dMRI-based techniques, such as those used to estimate inner axon diameters [37, 109].

iii. This study does not include a numerical evaluation of the model’s robustness to noise and artifacts in dMRI data. The numerical stability depends on the specific dMRI sequence and experimental parameters, such as diffusion gradient strength, diffusion times, and TE. For example, combining diffusion-weighted and double-inversion recovery sequences optimized to suppress non-myelin water signals would enable direct fitting of the proposed models to the measured data. In contrast, diffusion-T2 acquisitions require a multi-compartment model incorporating the proposed methodology. In future work, we plan to address these issues, employing Cramér-Rao bound analyses to optimize acquisition parameters for different sequences and evaluate the fitting stability under varying noise levels.

iv. All results presented in this study are based on synthetic signals derived from the proposed analytical models or MC simulations. Validation with real dMRI data, including histological analyses of various brain regions, is crucial for future work. Additionally, the diffusivity values used in this study are based on those reported by [65]. Still, variations in reported myelin water diffusivities in other experimental [63] and numerical studies [110–112] suggest the need for further work to reconcile these discrepancies and identify more accurate ex vivo and in vivo myelin water diffusivities.

v. Our MC simulations and proposed models assume straight cylinders, thus neglecting axonal undulations and beading, which are known to influence diffusion in WM [99, 113–115]. Incorporating more realistic axonal geometries constitutes a critical direction for future research, as modeling these effects could enhance the generalizability of our approach. To address these limitations, we plan to conduct numerical evaluations to assess their impact on the estimated effective myelin radius and adapt the models to include geometrical variations informed by histological data. Axonal undulations and beading are expected to reduce the apparent parallel diffusivity and increase the radial diffusivity of myelin water relative to values observed for straight cylinders. Based on the relationship between the radius and myelin water diffusivities provided by the Gaussian approximation in Equation (13), these effects would likely lead to overestimating the effective myelin radius compared to the actual value.

vi. In severe pathological conditions, such as certain multiple sclerosis lesions, where the myelin sheath breaks down and undergoes vacuolization, leading to the separation of adjacent spirals as well as axonal dissociation and degeneration [116, 117], the assumptions underlying the proposed model are no longer valid. In such cases, increased water permeability and alterations in myelin water layer thickness would compromise the applicability of the proposed formalism. Therefore, this model is likely more suited for studying healthy brains and pathological conditions at earlier stages with milder alterations.

vii. All data were generated based on an acquisition protocol potentially feasible with a Connectome 2.0-like human scanner equipped with a diffusion gradient of 500 mT/m, where the TE can be further reduced by employing an image readout technique starting at the center of k-space (e.g., spiral). Future studies should investigate a range of acquisition protocols, including stronger diffusion gradients available in preclinical scanners [37], as well as the 300 mT/m diffusion gradients utilized in the Connectome 1.0 [67, 109] and GE SIGNA MAGNUS scanners. The recently introduced MAGNETOM Cima. X clinical scanner, with a diffusion gradient strength of 200 mT/m, should also be considered. Determining the optimal acquisition parameters for each scenario is crucial for improving sensitivity to myelin sheath radii.

In summary, this work introduces dMRI models capable of characterizing myelin water diffusion, enabling the estimation of the effective myelin sheath radius per voxel. This water pool has been largely overlooked in previous dMRI studies due to the strong signal suppression it experiences when long TEs are used in clinical applications due to its short T2 relaxation time. However, recent advancements in dMRI sequences and the advent of MRI scanners equipped with stronger diffusion gradients make it possible to acquire dMRI signals significantly weighted by myelin water. This progress underscores the importance of having available models for this specific tissue compartment.

Nevertheless, the applicability of the proposed methodology is limited by hardware availability. Its use is restricted to a few human scanners with strong diffusion gradients and preclinical animal scanners with higher gradient strengths (e.g., G = 300–1,500 mT/m). This limitation highlights the need for broader access to such advanced MRI systems to fully exploit the potential of these models for both research and clinical applications. Additionally, pathologies involving vacuolization of myelin sheaths or significant separation of adjacent spirals result in altered myelin water layer thickness and increased permeability, which could compromise the validity of the proposed formalism. Consequently, the model is best suited for studies of healthy brains and pathological conditions at earlier stages, where tissue alterations are less severe.

By addressing the discussed limitations and validating the models with real dMRI data and histological measurements, future research may enhance the accuracy and applicability of the proposed models, contributing to the development of novel MRI biomarkers of WM tissue microstructure.

Data availability statement

The datasets and code presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/ejcanalesr/myelin-water-diffusion-models.

Author contributions

EC-R: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing. CT: Investigation, Methodology, Resources, Software, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing. EF-G: Investigation, Resources, Writing–original draft, Writing–review and editing. DJ: Investigation, Resources, Writing–original draft, Writing–review and editing. J-PT: Investigation, Resources, Writing–original draft, Writing–review and editing. JR-P: Data curation, Formal Analysis, Investigation, Methodology, Resources, Software, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. EC-R was supported by the Swiss National Science Foundation (SNSF), Ambizione fellowship PZ00P2_185814 and SNSF grant number 10000706. CT is supported by a Sir Henry Wellcome Fellowship (215944/Z/19/Z), and EF-G is supported by the SNSF, grant number 10000706.

Acknowledgments

To facilitate open access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript arising from this submission.

Conflict of interest

The authors declare that the research was conducted without any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The authors declare that Generative AI was used in the creation of this manuscript. We used ChatGPT (GPT-3.5, free version) and Grammarly (premium) to assist in identifying grammatical errors and typos in this manuscript. All intellectual contributions, including the development of ideas, analysis, and interpretation, were made solely by the authors.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2025.1516630/full#supplementary-material

References

1. Kandel ER, Koester JD, Mack SH, Siegelbaum SA. Principles of neural science. 6th ed. (2021).

Google Scholar

2. Aboitiz F, Scheibel AB, Fisher RS, Zaidel E. Fiber composition of the human corpus callosum. Brain Res (1992) 598:143–53. doi:10.1016/0006-8993(92)90178-C

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Goldstein SS, Rall W. Changes of action potential shape and velocity for changing core conductor geometry. Biophys J (1974) 14:731–57. doi:10.1016/S0006-3495(74)85947-3

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Hursh JB. Conduction velocity and diameter of nerve fibers. Am J Physiol Content (1939) 127:131–9. doi:10.1152/ajplegacy.1939.127.1.131

CrossRef Full Text | Google Scholar

5. Laule C, Leung E, Li DKB, Traboulsee AL, Paty DW, MacKay AL, et al. Myelin water imaging in multiple sclerosis: quantitative correlations with histopathology. Mult Scler (2006) 12:747–53. doi:10.1177/1352458506070928

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Chen J, Patel Z, Liu S, Bock NA, Frey BN, Suh JS. A systematic review of abnormalities in intracortical myelin across psychiatric illnesses. J Affect Disord Rep (2024) 15:100689. doi:10.1016/J.JADR.2023.100689

CrossRef Full Text | Google Scholar

7. Valdés-Tovar M, Rodríguez-Ramírez AM, Rodríguez-Cárdenas L, Sotelo-Ramírez CE, Camarena B, Sanabrais-Jiménez MA, et al. Insights into myelin dysfunction in schizophrenia and bipolar disorder. World J Psychiatry (2022) 12:264–85. doi:10.5498/wjp.v12.i2.264

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Papuc E, Rejdak K. The role of myelin damage in Alzheimer’s disease pathology. Arch Med Sci (2020) 16:345–1. doi:10.5114/AOMS.2018.76863

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Gong Z, Bilgel M, Kiely M, Triebswetter C, Ferrucci L, Resnick SM, et al. Lower myelin content is associated with more rapid cognitive decline among cognitively unimpaired individuals. Alzheimers Dement (2023) 19:3098–107. doi:10.1002/ALZ.12968

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Alexander DC, Dyrby TB, Nilsson M, Zhang H. Imaging brain microstructure with diffusion MRI: practicality and applications. NMR Biomed (2019) 32:e3841. doi:10.1002/nbm.3841

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Dyrby TB, Innocenti GM, Bech M, Lundell H. Validation strategies for the interpretation of microstructure imaging using diffusion MRI. Neuroimage (2018) 182:62–79. doi:10.1016/j.neuroimage.2018.06.049

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Novikov DS, Fieremans E, Jespersen SN, Kiselev VG. Quantifying brain microstructure with diffusion MRI: theory and parameter estimation. NMR in Biomedicine (2019) 32:e3998. doi:10.1002/nbm.3998

CrossRef Full Text | Google Scholar

13. Jeurissen B, Tournier JD, Dhollander T, Connelly A, Sijbers J. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. Neuroimage (2014) 103:411–26. doi:10.1016/j.neuroimage.2014.07.061

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Tournier JD, Yeh CH, Calamante F, Cho KH, Connelly A, Lin CP. Resolving crossing fibres using constrained spherical deconvolution: validation using diffusion-weighted imaging phantom data. Neuroimage (2008) 42:617–25. doi:10.1016/j.neuroimage.2008.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Dell’Acqua F, Scifo P, Rizzo G, Catani M, Simmons A, Scotti G A modified damped Richardson-Lucy algorithm to reduce isotropic background effects in spherical deconvolution. Neuroimage (2010) 49:1446–58. doi:10.1016/j.neuroimage.2009.09.033

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Canales-Rodríguez EJ, Daducci A, Sotiropoulos SN, Caruyer E, Aja-Fernández S, Radua J, et al. Spherical deconvolution of multichannel diffusion MRI data with non-Gaussian noise models and spatial regularization. PLoS One (2015) 10:e0138910. doi:10.1371/journal.pone.0138910

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Canales-Rodríguez EJ, Legarreta JH, Pizzolato M, Rensonnet G, Girard G, Patino JR-, et al. Sparse wars: a survey and comparative study of spherical deconvolution algorithms for diffusion MRI. Neuroimage (2019) 184:140–60. doi:10.1016/j.neuroimage.2018.08.071

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Daducci A, Canales-Rodriguez EJ, Descoteaux M, Garyfallidis E, Gur Y, Lin Y-CC, et al. Quantitative comparison of reconstruction methods for intra-voxel fiber recovery from diffusion MRI. IEEE Trans Med Imaging (2014) 33:384–99. doi:10.1109/TMI.2013.2285500

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Canales-Rodríguez EJ, Melie-García L, Iturria-Medina Y. Mathematical description of q-space in spherical coordinates: exact q-ball imaging. Magn Reson Med (2009) 61:1350–67. doi:10.1002/mrm.21917

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Tuch DS. Q-ball imaging. Magn Reson Med (2004) 52:1358–72. doi:10.1002/mrm.20279

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wedeen VJ, Wang RP, Schmahmann JD, Benner T, Tseng WYII, Dai G, et al. Diffusion spectrum magnetic resonance imaging (DSI) tractography of crossing fibers. Neuroimage (2008) 41:1267–77. doi:10.1016/j.neuroimage.2008.03.036

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Canales-Rodríguez EJ, Iturria-Medina Y, Alemán-Gómez Y, Melie-García L. Deconvolution in diffusion spectrum imaging. Neuroimage (2010) 50:136–49. doi:10.1016/j.neuroimage.2009.11.066

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Lippe S, Poupon C, Cachia A, Archambaud F, Rodrigo S, Dorfmuller G, et al. White matter abnormalities revealed by DTI correlate with interictal grey matter FDG-PET metabolism in focal childhood epilepsies. Epileptic Disord (2012) 14:404–13. doi:10.1684/epd.2012.0547

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Ozarslan E, Shepherd TM, Vemuri BC, Blackband SJ, Mareci TH. Resolution of complex tissue microarchitecture using the diffusion orientation transform (DOT). Neuroimage (2006) 31:1086–103. doi:10.1016/j.neuroimage.2006.01.024

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Behrens TEJ, Berg HJ, Jbabdi S, Rushworth MFS, Woolrich MW. Probabilistic diffusion tractography with multiple fibre orientations: what can we gain? Neuroimage (2007) 34:144–55. doi:10.1016/j.neuroimage.2006.09.018

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Sotiropoulos SN, Behrens TE, Jbabdi S. Ball and rackets: inferring fiber fanning from diffusion-weighted MRI. Neuroimage (2012) 60:1412–25. doi:10.1016/j.neuroimage.2012.01.056

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Ramirez-Manzanares A, Rivera M, Vemuri BC, Carney P, Mareci T. Diffusion basis functions decomposition for estimating white matter intravoxel fiber geometry. IEEE Trans Med Imaging (2007) 26:1091–102. doi:10.1109/TMI.2007.900461

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Melie-García L, Canales-Rodríguez EJ, Alemán-Gómez Y, Lin CP, Iturria-Medina Y, Valdés-Hernández PA A Bayesian framework to identify principal intravoxel diffusion profiles based on diffusion-weighted MR imaging. Neuroimage (2008) 42:750–70. doi:10.1016/j.neuroimage.2008.04.242

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Cheng J, Deriche R, Jiang T, Shen D, Yap PT. Non-negative spherical deconvolution (NNSD) for estimation of fiber orientation distribution function in single-/multi-shell diffusion MRI. Neuroimage (2014) 101:750–64. doi:10.1016/j.neuroimage.2014.07.062

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Alexander DC. Maximum entropy spherical deconvolution for diffusion MRI (2005) In: Information Processing in Medical Imaging. Editors G. E. Christensen, and M. Sonka Springer, Berlin, Heidelberg: Lecture Notes in Computer Science 3565. doi:10.1007/11505730_7

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain. Neuroimage (2012) 61:1000–16. doi:10.1016/j.neuroimage.2012.03.072

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Tariq M, Schneider T, Alexander DC, Gandini Wheeler-Kingshott CA, Zhang H. Bingham-NODDI: mapping anisotropic orientation dispersion of neurites using diffusion MRI. Neuroimage (2016) 133:207–23. doi:10.1016/j.neuroimage.2016.01.046

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Kaden E, Kelm ND, Carson RP, Does MD, Alexander DC. Multi-compartment microscopic diffusion imaging. Neuroimage (2016) 139:346–59. doi:10.1016/j.neuroimage.2016.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Panagiotaki E, Schneider T, Siow B, Hall MG, Lythgoe MF, Alexander DC. Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison. Neuroimage (2012) 59:2241–54. doi:10.1016/j.neuroimage.2011.09.081

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Behrens TE, Woolrich MW, Jenkinson M, Johansen-Berg H, Nunes RG, Clare S, et al. Characterization and propagation of uncertainty in diffusion-weighted MR imaging. Magn Reson Med (2003) 50:1077–88. doi:10.1002/mrm.10609

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Assaf Y, Alexander DC, Jones DK, Bizzi A, Behrens TE, Clark CA, et al. The CONNECT project: combining macro- and micro-structure. Neuroimage (2013) 80:273–82. doi:10.1016/j.neuroimage.2013.05.055

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Veraart J, Nunes D, Rudrapatna U, Fieremans E, Jones DK, Novikov DS, et al. Nonivasive quantification of axon radii using diffusion MRI. Elife (2020) 9:e49855. doi:10.7554/eLife.49855

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Assaf Y, Basser PJ. Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. Neuroimage (2005) 27:48–58. doi:10.1016/j.neuroimage.2005.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Dyrby TB, Søgaard LV, Hall MG, Ptito M, Alexander DC. Contrast and stability of the axon diameter index from microstructure imaging with diffusion MRI. Magn Reson Med (2013) 70:711–21. doi:10.1002/mrm.24501

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Assaf Y, Blumenfeld-Katzir T, Yovel Y, Basser PJ. AxCaliber: a method for measuring axon diameter distribution from diffusion MRI. Magn Reson Med (2008) 59:1347–54. doi:10.1002/mrm.21577

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Alexander DC, Hubbard PL, Hall MG, Moore EA, Ptito M, Parker GJ, et al. Orientationally invariant indices of axon diameter and density from diffusion MRI. Neuroimage (2010) 52:1374–89. doi:10.1016/j.neuroimage.2010.05.043

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Fan Q, Nummenmaa A, Witzel T, Ohringer N, Tian Q, Setsompop K, et al. Axon diameter index estimation independent of fiber orientation distribution using high-gradient diffusion MRI. Neuroimage (2020) 222:117197. doi:10.1016/j.neuroimage.2020.117197

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Pizzolato M, Canales-Rodríguez EJ, Andersson M, Dyrby TB. Axial and radial axonal diffusivities and radii from single encoding strongly diffusion-weighted MRI. Med Image Anal (2023) 86:102767. doi:10.1016/J.MEDIA.2023.102767

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Barakovic M, Girard G, Schiavi S, Romascano D, Descoteaux M, Granziera C, et al. Bundle-specific axon diameter index as a new contrast to differentiate white matter tracts. Front Neurosci (2021) 15:646034. doi:10.3389/FNINS.2021.646034

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Barakovic M, Pizzolato M, Tax CMW, Rudrapatna U, Magon S, Dyrby TB, et al. Estimating axon radius using diffusion-relaxation MRI: calibrating a surface-based relaxation model with histology. Front Neurosci (2023) 17:1209521. doi:10.3389/FNINS.2023.1209521

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Canales-Rodríguez EJ, Pizzolato M, Zhou FL, Barakovic M, Thiran JP, Jones DK, et al. Pore size estimation in axon-mimicking microfibers with diffusion-relaxation MRI. Magn Reson Med (2024) 91:2579–96. doi:10.1002/MRM.29991

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Jelescu IO, Zurek M, Winters KV, Veraart J, Rajaratnam A, Kim NS, et al. In vivo quantification of demyelination and recovery using compartment-specific diffusion MRI metrics validated by electron microscopy. Neuroimage (2016) 132:104–14. doi:10.1016/j.neuroimage.2016.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Ramanna S, Moss HG, McKinnon ET, Yacoub E, Helpern JA, Jensen JH. Triple diffusion encoding MRI predicts intra-axonal and extra-axonal diffusion tensors in white matter. Magn Reson Med (2020) 83:2209–20. doi:10.1002/MRM.28084

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Veraart J, Novikov DS, Fieremans E. TE dependent Diffusion Imaging (TEdDI) distinguishes between compartmental T2 relaxation times. Neuroimage (2018) 182:360–9. doi:10.1016/J.NEUROIMAGE.2017.09.030

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Barakovic M, Tax CMW, Rudrapatna U, Chamberland M, Rafael-Patino J, Granziera C, et al. Resolving bundle-specific intra-axonal T2 values within a voxel using diffusion-relaxation tract-based estimation. Neuroimage (2021) 227:117617. doi:10.1016/j.neuroimage.2020.117617

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Canales-Rodríguez EJ, Pizzolato M, Yu T, Piredda GF, Hilbert T, Radua J, et al. Revisiting the T2 spectrum imaging inverse problem: bayesian regularized non-negative least squares. Neuroimage (2021) 244:118582. doi:10.1016/j.neuroimage.2021.118582

PubMed Abstract | CrossRef Full Text | Google Scholar

52. MacKay A, Laule C, Vavasour I, Bjarnason T, Kolind S, Mädler B. Insights into brain microstructure from the T2 distribution. Magn Reson Imaging (2006) 24:515–25. doi:10.1016/j.mri.2005.12.037

PubMed Abstract | CrossRef Full Text | Google Scholar

53. MacKay AL, Laule C. Magnetic resonance of myelin water: an in vivo marker for myelin. Brain Plast (2016) 2:71–91. doi:10.3233/BPL-160033

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Piredda GF, Hilbert T, Canales-Rodríguez EJ, Pizzolato M, Meuli R, Pfeuffer J, et al. Fast and high-resolution myelin water imaging: Accelerating multi-echo GRASE with CAIPIRINHA. Magn Reson Med (2021) 85(1):209–222. Available from: http://archive.ismrm.org/2019/4400.html.

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Kumar D, Hariharan H, Faizy TD, Borchert P, Siemonsen S, Fiehler J, et al. Using 3D spatial correlations to improve the noise robustness of multi component analysis of 3D multi echo quantitative T2 relaxometry data. Neuroimage (2018) 178:583–601. doi:10.1016/j.neuroimage.2018.05.026

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Dvorak AV, Ljungberg E, Vavasour IM, Lee LE, Abel S, Li DKB, et al. Comparison of multi echo T2 relaxation and steady state approaches for myelin imaging in the central nervous system. Sci Rep (2021) 11:1369. doi:10.1038/s41598-020-80585-7

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Prasloski T, Rauscher A, MacKay AL, Hodgson M, Vavasour IM, Laule C, et al. Rapid whole cerebrum myelin water imaging using a 3D GRASE sequence. Neuroimage (2012) 63:533–9. doi:10.1016/j.neuroimage.2012.06.064

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Raj A, Pandya S, Shen X, LoCastro E, Nguyen TD, Gauthier SA. Multi-compartment T2 relaxometry using a spatially constrained multi-Gaussian model. PLoS One (2014) 9:e98391. doi:10.1371/journal.pone.0098391

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Alonso-Ortiz E, Levesque IR, Pike GB. MRI-based myelin water imaging: a technical review. Magn Reson Med (2015) 73:70–81. doi:10.1002/mrm.25198

PubMed Abstract | CrossRef Full Text | Google Scholar

60. 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:673–7. doi:10.1002/mrm.1910310614

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Canales-Rodríguez EJ, Pizzolato M, Piredda GFGF, Hilbert T, Kunz N, Pot C, et al. Comparison of non-parametric T2 relaxometry methods for myelin water quantification. Med Image Anal (2021) 69:101959. doi:10.1016/j.media.2021.101959

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Canales-Rodríguez EJ, Alonso-Lana S, Verdolini N, Sarró S, Feria I, Montoro I, et al. Age- and gender-related differences in brain tissue microstructure revealed by multi-component T2 relaxometry. Neurobiol Aging (2021) 106:68–79. doi:10.1016/J.NEUROBIOLAGING.2021.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Stanisz GJ, Henkelman RM. Diffusional anisotropy of T2 components in bovine optic nerve. Magn Reson Med (1998) 40:405–10. doi:10.1002/MRM.1910400310

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Peled S, Cory DG, Raymond SA, Kirschner DA, Jolesz FA. Water diffusion, T2, and compartmentation in frog sciatic nerve. Magn Reson Med (1999) 42:911–8. doi:10.1002/(sici)1522-2594(199911)42:5<911::aid-mrm11>3.0.co;2-j

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Andrews TJ, Osborne MT, Does MD. Diffusion of myelin water. Magn Reson Med (2006) 56:381–5. doi:10.1002/MRM.20945

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Avram AV, Guidon A, Song AW. Myelin water weighted diffusion tensor imaging. Neuroimage (2010) 53:132–8. doi:10.1016/j.neuroimage.2010.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Jones DK, Alexander DC, Bowtell R, Cercignani M, Dell’Acqua F, McHugh DJ, et al. Microstructural imaging of the human brain with a “super-scanner”: 10 key advantages of ultra-strong gradients for diffusion MRI. Neuroimage (2018) 182:8–38. doi:10.1016/j.neuroimage.2018.05.047

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Huang SY, Witzel T, Keil B, Scholz A, Davids M, Dietz P, et al. Connectome 2.0: developing the next-generation ultra-high gradient strength human MRI scanner for bridging studies of the micro-meso- and macro-connectome. Neuroimage (2021) 243:118530. doi:10.1016/J.NEUROIMAGE.2021.118530

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Tax CMW, Kleban E, Chamberland M, Baraković M, Rudrapatna U, Jones DK. Measuring compartmental T2-orientational dependence in human brain white matter using a tiltable RF coil and diffusion-T2 correlation MRI. Neuroimage (2021) 236:117967. doi:10.1016/j.neuroimage.2021.117967

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Mueller L, Tax CMW, Jones DK. Unprecedented echo times for diffusion MRI using connectom gradients, spiral readouts and field monitoring. In: MAGNETOM flash (2019). Available from: https://cdn0.scrvt.com/39b415fb07de4d9656c7b516d8e2d907/1800000006277360/3a1f90ead087/siemens-healthineers-magnetom-flash-74-ismrm-unprecedented-echo-times_1800000006277360.pdf (Accessed December 14, 2023).

Google Scholar

71. Tax CMW, Rudrapatna US, Mueller L. Characterizing diffusion of myelin water in the living human brain using ultra-strong gradients and spiral readout. Proc 27th (2019). Available from: https://cds.ismrm.org/protected/19MProceedings/PDFfiles/1115.html.

Google Scholar

72. Wilm BJ, Hennel F, Roesler MB, Weiger M, Pruessmann KP. Minimizing the echo time in diffusion imaging using spiral readouts and a head gradient system. Magn Reson Med (2020) 84:3117–27. doi:10.1002/MRM.28346

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Inouye H, Kirschner DA. Membrane interactions in nerve myelin. I. Determination of surface charge from effects of pH and ionic strength on period. Biophys J (1988) 53:235–45. doi:10.1016/S0006-3495(88)83085-6

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Nilsson M, Lasič S, Drobnjak I, Topgaard D, Westin C-FF, Lasǐ S, et al. Resolution limit of cylinder diameter estimation by diffusion MRI: the impact of gradient waveform and orientation dispersion. NMR Biomed (2017) 30:e3711. doi:10.1002/nbm.3711

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Assaf Y, Freidlin RZ, Rohde GK, Basser PJ. New modeling and experimental framework to characterize hindered and restricted water diffusion in brain white matter. Magn Reson Med (2004) 52:965–78. doi:10.1002/mrm.20274

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Stejskal EO, Tanner JE. Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J Chem Phys (1965) 42:288–92. doi:10.1063/1.1695690

CrossRef Full Text | Google Scholar

77. Mardia KV, Jupp PE. Directional statistics. J. Wiley (2008). doi:10.1002/9780470316979

CrossRef Full Text | Google Scholar

78. Watson GS. Distributions on the circle and sphere. J Appl Probab (1982) 19:265–80. doi:10.2307/3213566

CrossRef Full Text | Google Scholar

79. Jammalamadaka SR, Kozubowski TJ. A general approach for obtaining wrapped circular distributions via mixtures. Sankhya (2017) 79:133–57. doi:10.1007/s13171-017-0096-4

CrossRef Full Text | Google Scholar

80. Ledesma-Motolinía M, Carrillo-Estrada JL, Escobar A, Donado F, Castro-Villarreal P. Magnetized granular particles running and tumbling on the circle S1. Phys Rev E (2023) 107:024902. doi:10.1103/PhysRevE.107.024902

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Castro-Villarreal P, Villada-Balbuena A, Méndez-Alcaraz JM, Castañeda-Priego R, Estrada-Jiménez S. A Brownian dynamics algorithm for colloids in curved manifolds. J Chem Phys (2014) 140:214115. doi:10.1063/1.4881060

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Stephens MA. Random walk on a circle. Biometrika (1963) 50:385. doi:10.2307/2333907

CrossRef Full Text | Google Scholar

83. Ozarslan E, Koay CG, Basser PJ. Remarks on q-space MR propagator in partially restricted, axially-symmetric, and isotropic environments. Magn Reson Imaging (2009) 27:834–44. doi:10.1016/j.mri.2009.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Anderson AW. Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magn Reson Med (2005) 54:1194–206. doi:10.1002/mrm.20667

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Lori NF, Conturo TE, Le Bihan D. Definition of displacement probability and diffusion time in q-space magnetic resonance measurements that use finite-duration diffusion-encoding gradients. J Magn Reson (2003) 165:185–95. doi:10.1016/j.jmr.2003.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Mattiello J, Basser PJ, Lebihan D. Analytical expressions for the b matrix in NMR diffusion imaging and spectroscopy. J Magn Reson Ser (1994) A(108):131–41. doi:10.1006/JMRA.1994.1103

CrossRef Full Text | Google Scholar

87. Kroenke CD, Ackerman JJH, Yablonskiy DA. On the nature of the NAA diffusion attenuated MR signal in the central nervous system. Magn Reson Med (2004) 52:1052–9. doi:10.1002/mrm.20260

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Edén M. Computer simulations in solid-state NMR. III. Powder averaging. Concepts Magn Reson Part (2003) A(18A):24–55. doi:10.1002/CMR.A.10065

CrossRef Full Text | Google Scholar

89. Lasič S, Szczepankiewicz F, Eriksson S, Nilsson M, Topgaard D. Microanisotropy imaging: quantification of microscopic diffusion anisotropy and orientational order parameter by diffusion MRI with magic-angle spinning of the q-vector. Front Phys (2014) 2. doi:10.3389/fphy.2014.00011

CrossRef Full Text | Google Scholar

90. Jensen JH, Russell Glenn G, Helpern JA. Fiber ball imaging. Neuroimage (2016) 124:824–33. doi:10.1016/j.neuroimage.2015.09.049

PubMed Abstract | CrossRef Full Text | Google Scholar

91. Caminiti R, Ghaziri H, Galuske R, Hof PR, Innocenti GM. Evolution amplified processing with temporally dispersed slow neuronal connectivity in primates. Proc Natl Acad Sci U S A (2009) 106:19551–6. doi:10.1073/PNAS.0907655106

PubMed Abstract | CrossRef Full Text | Google Scholar

92. Rafael-Patino J, Romascano D, Ramirez-Manzanares A, Canales-Rodríguez EJ, Girard G, Thiran JP. Robust monte-carlo simulations in diffusion-MRI: effect of the substrate complexity and parameter choice on the reproducibility of results. Front Neuroinform (2020) 14:8. doi:10.3389/fninf.2020.00008

PubMed Abstract | CrossRef Full Text | Google Scholar

93. Chomiak T, Hu B. What is the optimal value of the g-ratio for myelinated fibers in the rat CNS? A theoretical approach. PLoS One (2009) 4:e7754. doi:10.1371/JOURNAL.PONE.0007754

PubMed Abstract | CrossRef Full Text | Google Scholar

94. Campbell JSW, Leppert IR, Narayanan S, Boudreau M, Duval T, Cohen-Adad J, et al. Promise and pitfalls of g-ratio estimation with MRI. Neuroimage (2017) 182:80–96. doi:10.1016/j.neuroimage.2017.08.038

PubMed Abstract | CrossRef Full Text | Google Scholar

95. Linse P, Söderman O. The validity of the short-gradient-pulse approximation in NMR studies of restricted diffusion. Simulations of molecules diffusing between planes, in cylinders and spheres. J Magn Reson Ser A (1995) 116:77–86. doi:10.1006/JMRA.1995.1192

CrossRef Full Text | Google Scholar

96. Söderman O, Jönsson B. Restricted diffusion in cylindrical geometry. J Magn Reson (1995) 94–7. doi:10.1006/jmra.1995.0014

CrossRef Full Text | Google Scholar

97. Bar-Shir A, Avram L, Ozarslan E, Basser PJ, Cohen Y. The effect of the diffusion time and pulse gradient duration ratio on the diffraction pattern and the structural information estimated from q-space diffusion MR: experiments and simulations. J Magn Reson (2008) 194:230–6. doi:10.1016/j.jmr.2008.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

98. Kaden E, Kruggel F, Alexander DC. Quantitative mapping of the per-axon diffusion coefficients in brain white matter. Magn Reson Med (2016) 75:1752–63. doi:10.1002/mrm.25734

PubMed Abstract | CrossRef Full Text | Google Scholar

99. Andersson M, Pizzolato M, Kjer HM, Skodborg KF, Lundell H, Dyrby TB. Does powder averaging remove dispersion bias in diffusion MRI diameter estimates within real 3D axonal architectures? Neuroimage (2022) 248:118718. doi:10.1016/j.neuroimage.2021.118718

PubMed Abstract | CrossRef Full Text | Google Scholar

100. Daducci A, Canales-Rodríguez EJ, Zhang H, Dyrby TB, Alexander DC, Thiran J-PP. Accelerated microstructure imaging via convex optimization (AMICO) from diffusion MRI data. Neuroimage (2015) 105:32–44. doi:10.1016/j.neuroimage.2014.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

101. Ferizi U, Schneider T, Witzel T, Wald LL, Zhang H, Wheeler-Kingshott CAM, et al. White matter compartment models for in vivo diffusion MRI at 300mT/m. Neuroimage (2015) 118:468–83. doi:10.1016/j.neuroimage.2015.06.027

PubMed Abstract | CrossRef Full Text | Google Scholar

102. Jelescu IO, Veraart J, Fieremans E, Novikov DS. Degeneracy in model parameter estimation for multi-compartmental diffusion in neuronal tissue. NMR Biomed (2016) 29:33–47. doi:10.1002/nbm.3450

PubMed Abstract | CrossRef Full Text | Google Scholar

103. Harkins KD, Does MD. Simulations on the influence of myelin water in diffusion-weighted imaging. Phys Med Biol (2016) 61:4729–45. doi:10.1088/0031-9155/61/13/4729

PubMed Abstract | CrossRef Full Text | Google Scholar

104. Brusini L, Menegaz G, Nilsson M. Monte Carlo simulations of water exchange through myelin wraps: implications for diffusion MRI. IEEE Trans Med Imaging (2019) 38:1438–45. doi:10.1109/TMI.2019.2894398

PubMed Abstract | CrossRef Full Text | Google Scholar

105. Mancini M, Karakuzu A, Cohen-Adad J, Cercignani M, Nichols TE, Stikov N. An interactive meta-analysis of MRI biomarkers of Myelin. Elife (2020) 9:e61523–23. doi:10.7554/elife.61523

PubMed Abstract | CrossRef Full Text | Google Scholar

106. Sled JG. Modelling and interpretation of magnetization transfer imaging in the brain. Neuroimage (2018) 182:128–35. doi:10.1016/J.NEUROIMAGE.2017.11.065

PubMed Abstract | CrossRef Full Text | Google Scholar

107. Caprihan A, Wang LZ, Fukushima E. A multiple-narrow-pulse approximation for restricted diffusion in a time-varying field gradient. J Magn Reson Ser (1996) 118:94–102. doi:10.1006/JMRA.1996.0013

CrossRef Full Text | Google Scholar

108. Callaghan PT. A simple matrix formalism for spin echo analysis of restricted diffusion under generalized gradient waveforms. J Magn Reson (1997) 129:74–84. doi:10.1006/jmre.1997.1233

PubMed Abstract | CrossRef Full Text | Google Scholar

109. Veraart J, Raven EP, Edwards LJ, Weiskopf N, Jones DK. The variability of MR axon radii estimates in the human white matter. Hum Brain Mapp (2021) 42:2201–13. doi:10.1002/hbm.25359

PubMed Abstract | CrossRef Full Text | Google Scholar

110. Sen PN, Basser PJ. A model for diffusion in white matter in the brain. Biophys J (2005) 89:2927–38. doi:10.1529/BIOPHYSJ.105.063016

PubMed Abstract | CrossRef Full Text | Google Scholar

111. Harkins KD, Dula AN, Does MD. Effect of intercompartmental water exchange on the apparent myelin water fraction in multiexponential T2 measurements of rat spinal cord. Magn Reson Med (2012) 67:793–800. doi:10.1002/MRM.23053

PubMed Abstract | CrossRef Full Text | Google Scholar

112. Baxter GT, Frank LR. A computational model for diffusion weighted imaging of myelinated white matter. Neuroimage (2013) 75:204–12. doi:10.1016/J.NEUROIMAGE.2013.02.076

PubMed Abstract | CrossRef Full Text | Google Scholar

113. Lee HH, Jespersen SN, Fieremans E, Novikov DS. The impact of realistic axonal shape on axon diameter estimation using diffusion MRI. Neuroimage (2020) 223:117228. doi:10.1016/j.neuroimage.2020.117228

PubMed Abstract | CrossRef Full Text | Google Scholar

114. Lee HH, Tian Q, Sheft M, Coronado-Leija R, Ramos-Llorden G, Abdollahzadeh A, et al. The effects of axonal beading and undulation on axonal diameter estimation from diffusion MRI: insights from simulations in human axons segmented from three-dimensional electron microscopy. NMR Biomed (2024) 37:e5087. doi:10.1002/NBM.5087

PubMed Abstract | CrossRef Full Text | Google Scholar

115. Andersson M, Kjer HM, Rafael-Patino J, Pacureanu A, Pakkenberg B, Thiran JP, et al. Axon morphology is modulated by the local environment and impacts the noninvasive investigation of its structure–function relationship. Proc Natl Acad Sci U S A (2020) 117:33649–59. doi:10.1073/pnas.2012533117

PubMed Abstract | CrossRef Full Text | Google Scholar

116. Barkovich AJ. Concepts of myelin and myelination in neuroradiology. AJNR Am J Neuroradiol (2000) 21(6):1099–109. Available from: https://pmc.ncbi.nlm.nih.gov/articles/PMC7973874/(Accessed January 6, 2025).

PubMed Abstract | Google Scholar

117. Raine CS, Cannella B, Hauser SL, Genain CP. Demyelination in primate autoimmune encephalomyelitis and acute multiple sclerosis lesions: a case for antigen-specific antibody mediation. Ann Neurol (1999) 46:144–60. doi:10.1002/1531-8249(199908)46:2<144::aid-ana3>3.0.co;2-k

PubMed Abstract | CrossRef Full Text | Google Scholar