Improving the Realism of White Matter Numerical Phantoms: A Step toward a Better Understanding of the Influence of Structural Disorders in Diffusion MRI

1. Introduction

Diffusion-weighted magnetic resonance imaging (dMRI), sensitized to the diffusive motion of water along the direction of an applied magnetic field gradient, has become a well-established technique to non-invasively probe the cellular organization of tissues in vivo. Diffusion NMR measurements embed some information about the inhibition of particles motion due to the presence of barriers in the local environment, and can therefore be exploited to map some specific microstructural features characterizing the brain white matter ultrastructure at cellular scales. While Pulsed Gradient Spin Echo (PGSE) sequences [1] are still widely used in clinical routine, alternative schemes such as Oscillating Gradient Spin Echo (OGSE) sequences seem a promising approach since they enable to explore the diffusion pattern in the frequency domain dual to the diffusion time domain and are able to probe shorter diffusion times compared to conventional PGSE. Some OGSE studies have reported that, at frequencies below 400 Hz, the OG-measured extra-axonal diffusivity transverse to axons in white matter is linearly dependent on the frequency of the employed OGSE sequence [2], whereas state-of-the art multi-compartment models of white matter relying on PGSE or OGSE sequences usually assume a Gaussian diffusion in the extra-axonal space [3–6]. A theoretical explanation was given in Burcaw et al. [7], where the observed frequency-dependence is interpreted as resulting from the extra-axonal 2D short-range disorder of axonal packings in the plane transverse to white matter fibers.

The extra-axonal perpendicular diffusivity transverse to axons was thus written as

$\begin{array}{ll}D(t)=D_{\infty }+A.\frac{ln(t/t_c)}{t},t_c\ll t& (1)\end{array}$

and, equivalently, in the frequency domain:

$\begin{array}{ll}\text{Re}{D}(\omega )~D_{\infty }+A.\frac{\pi }{2}\left|\omega \right|,\left|\omega \right|t_c\ll 1& (2)\end{array}$

where t is the diffusion time, ω is the dual frequency and t_c represents the time to diffuse across the correlation length l_c of the packing geometry (l_c closely follows the mean external radius r_ext of the axon packing [7, 8]).

Similar to recent models accounting for extra-axonal time-dependence in the case of single diffusion encoding sequences [9], a multi-compartment model for cosine OGSE sequences was proposed in Ginsburger et al. [10] which added a frequency-dependent term in the extra-axonal diffusion tensor perpendicular diffusivity based on Equation (2), showing a significant improvement of the model fit quality.

The first contribution of this article is to show the relevance of such a frequency-dependent correction using clinically feasible cosine trapezoidal OGSE sequences. The theoretically predicted linear relationship between the extracellular perpendicular diffusivity and the OGSE frequency was observed using 3D simulations of the diffusion process with different values of signal-to-noise ratios (SNR).

Having introduced a physically plausible frequency-dependent correction in our model, the next step is to study the dependence of its scaling coefficient A (Equation 2) to the geometrical features of the extracellular space. Indeed, to our knowledge, this dependence is very little known. An empirical law $A~l_c^2$ relating A to the correlation length l_c was given in Burcaw et al. [7] and Fieremans et al. [8] but is not sufficient to catch the complexity of the scaling coefficient A. The value of A is a measure of the strength of the structural disorder [7, 11], thus related to the geometrical properties characterizing the spatial organization of white matter at various scales. A possible approach to decipher the complex relationship between A and white matter features is to perform Monte-Carlo simulations of the diffusion process in diffusing media with increasing level of structural disorder.

The simulation of the diffusion process in state-of-the-art Monte-Carlo simulators such as CAMINO [12] and DMS [13] is decomposed into three main steps: (1) the generation of a 3D numerical phantom representing the diffusion medium (2) the Monte-Carlo simulation of the Brownian motion of spins (3) the synthesis of a DW-NMR signal. Simulators like CAMINO or DMS are able to extract triangle meshes from histological samples in step (1) in order to simulate diffusion in ultra-realistic media. However, this approach does not allow to have access to the plethora of possible geometries but only to the limited set of configurations provided by the reduced set of histological samples. Beside the possibility to use geometries extracted from histological samples in step (1), state-of-the-art simulators can only generate a limited number of geometries which might not represent white matter sufficiently well. For instance, the CAMINO and DMS simulators are able to simulate the diffusion process in any triangle mesh, but the algorithm used to construct simulation meshes from input geometrical parameters only generates substrates with straight cylinders of various diameters (including crossing between two populations of fibers). Other simulation tools like Fiberfox [14] rely on analytical models associated to each particular cell geometry, including various combinations of sticks, tensors, zeppelins, balls, dots, and astrosticks. They are inherently limited by the realism of the used geometries and the employed analytical models do not account for the presence of structural disorder in the extracellular space. There is thus a real need to propose alternative generative algorithms able to create more complex geometries while controlling the parameters driving the various sources of geometrical disorder to explore more extensively the vast domain of possible geometries.

The main contribution of this article is therefore the development of a novel algorithm to produce a wide variety of biomimicking numerical phantoms representing more realistic white matter tissue configurations from a reduced set of control parameters. Embedded in the Diffusion Microscopist Simulator (DMS) [13], this algorithm enables to control the degree of complexity of the generated geometrical configurations (induction of global angular dispersion and local tortuosity, presence of Ranvier nodes along the axonal membrane, presence of beading) with few design parameters, and without the necessity of any input histological sample. Such numerical white matter phantoms are then used to perform Monte-Carlo simulations of the diffusion process from which simulated diffusion-weighted NMR signals can be synthesized using trapezoidal OGSE sequences at different values of SNR. The obtained signal is fed into our analytical model to explore the evolution of the structural disorder coefficient A for various well-characterized geometrical configurations of the extracellular space.

2. Materials and Methods

2.1. White Matter Numerical Phantoms

Actual simulation tools rely either on membrane surfaces extracted from manual or automatic segmentations of microscopic imaging data stemming from histological tissue samples or on algorithms generating designs of axonal membranes built from a distribution of rectilinear cylinders (possibly with angular dispersion) which diameters d follow a Gamma distribution h. The function h can be tuned using a α shape parameter and a β inverse scale parameter such that:

$\begin{array}{ll}h(d;\alpha ,\beta )=\frac{{\beta }^{\alpha }{d}^{\alpha -1}{e}^{-\beta d}}{\Gamma (\alpha )}& (3)\end{array}$

where Γ is a complete Gamma function. The first approach is hard to achieve and generally yields a limited collection of membrane configurations, whereas the second approach provides over-simplistic representations of white matter. Here, we propose a novel framework to design more realistic membrane geometries better mimicking white matter structure. One of the challenging issue of such an approach is to design a tool enabling to cover a wide range of actual geometries while using a limited number of parameters. Several observations can be taken into account to improve the development of more realistic white matter geometries:

• Several heterogeneous populations of fibers can populate the field of view (FOV) of interest; whatever the target FOV (from mesoscale to millimeter scale), complex configurations of fibers are likely to happen and several studies in the field of diffusion MRI have reported a percentage of around 60 percent of voxels containing crossing, kissing or splitting fibers at the conventional millimeter resolution of diffusion MRI data [15].

• Each fiber population is composed of myelinated or unmyelinated axons which diameters follow the previous Gamma distribution (Equation 3) but the shape and inverse scale parameters can vary from a population to another; in addition, each population is characterized by its mean orientation in the 3D space and by its volume fraction.

• Myelinated axons are regularly interrupted by Ranvier nodes along the axon main direction [16]; the internode distance d has been extensively studied in Rushton [17] leading to the maximum conduction relationship

$\begin{array}{ll}\frac{d}{D}=kg\sqrt{log\left(\frac{1}{g}\right)}& (4)\end{array}$

where k is a constant, D is the external diameter of the axon (including the myelin sheath) and g is the g-ratio defined as the ratio between the axonal membrane and the external myelin sheath outer membrane diameters.

• The fibers of a given population depict a macroscopic angular dispersion that corresponds to the global misalignment of the axons which has previously been modeled both in ActiveAx [4] and NODDI [5] models using Watson's or Bingham's distributions relying on the knowledge of the principal orientation and of one or two concentration parameters respectively, thus imposing a strong assumption on the nature of angular dispersion.

• The fibers of a given population also depict local tortuosity than can be simply measured by the ratio between the geodesic distance along the curvilinear frame defined by the centroid axis of the fiber and the Euclidean distance between the two extremities of the fiber.

• It is not clear whether the axon diameter and myelin sheath thickness remain constant along the axon; several studies have assumed this absence of variation [18, 19] whereas there is no clear assessment of such a property; in particular, it is known that membrane injury can induce axonal beading for instance due to cytoskeletal damage. According to Budde and Frank [20], beading-induced changes in cell-membrane morphology are sufficient to significantly hinder water mobility and thereby decrease the apparent diffusion coefficient; it is therefore recommended to account for this and allow axon diameter variation.

Accounting for all these observations, we propose an algorithm relying on a six-fold strategy to design white matter mimicking numerical phantoms, that do not present the actual limitations of existing phantom design tools. Such phantoms will allow to go further into the study of the impact of both intra- and extra-axonal compartments on the diffusion signal. In particular, current achievable numerical phantoms do not provide the possibility to induce local fiber tortuosity, which might have a significant impact on the signal stemming from the extra-axonal compartment by modifying the parallel and transverse diffusivities along and perpendicular to the fibers.

The phantom generation algorithm takes a maximal number of 14N_Populations + 2 parameters to generate complex axonal geometries, where N_Populations is the number of fiber populations (N_Populations > 1 in the case of crossing fibers). The list of control parameters used to design a fiber population is summarized in Table 1.

TABLE 1

Table 1. List of control parameters names associated to each step of the phantom generation algorithm.

Step 1—During the first step of the phantom construction, similarly to state-of-the-art simulation tools, a set of over-simplistic fiber populations is constructed, each fiber population corresponding to a set of rectilinear and parallel outer envelopes. To each population corresponds an orientation, a Gamma distribution of fiber envelope diameters (defined by a mean diameter D and a standard deviation σ_D related to the shape α and scale β of the Gamma distribution, such that αβ = D and $\alpha {\beta }^2={\sigma }_D$), and an intracellular fraction, amounting to 4 parameters (Figure 1A). In the case of multiple populations, the degree of interweaving of axons from different populations is ruled by 2 scale parameters which control the distance between axons within the same population, thus enabling to create aggregated structures or a sheet organization [21] where axons find their way amongst other populations. At each step of the algorithm, the absence of intersection between the outer envelopes of fibers populating the phantom is ensured using a 3D collision algorithm.

FIGURE 1

Figure 1. (A) Phantom generation for two populations (blue and green) with orientations u_i (B). Global angular dispersion is created by selecting randomly an axon among a given fiber and a point on this axon which will be the center of rotation of the fiber. The selection of the rotation center follows a Gaussian distribution with mean μ and variance σ. A perturbation vector g whose components follow a Gaussian distribution with a variance proportional to the target angular dispersion is added to the orientation vector u, resulting in a rotation around the selected fixed point and in a new orientation vector u′. (C) Local angular dispersion is induced by deforming each axon separately. A point on the axon (in red) and a direct trieder (u_x, u_y, u_z) are chosen randomly. u_y defines the direction of the tortuosity deformation whose amplitude follows a Gaussian distribution with a variance depending on the number of points which are affected by the deformation around the central red point. (D) Creation of the myelin sheath. Inside each cylinder of radius R_tot, an inner cylinder of radius R = g.R_tot (g is the g-ratio) is created which represents the axonal membrane, and the external cylinder represents the outer layer of myelin sheath. (E) Creation of Ranvier nodes. The resolution of the fiber mesh around each Ranvier node is refined to better account for the exponential decay of the myelin thickness around the node. (F) Beading generation. Both the axonal contour (inner mesh) and the myelin sheath (external mesh) are swollen with a sine function. The myelin sheath thickness is preserved since beading comes from the swelling of the axonal membrane due to injury.

Step 2—Once the fiber populations composed of parallel, rectilinear axons are constructed, the second step of the algorithm consists in the induction of global orientation dispersion requiring one further parameter per fiber population. As depicted in Figure 1B, global angular dispersion is created by selecting randomly one fiber among a fiber population. The fiber population is randomly selected among those that did not reach their target angular dispersion yet. Then, a center of rotation is selected along this fiber, following a Gaussian distribution to ensure that most of the selected rotation centers belong to the central part of the fibers (see Figure 1B). A perturbation vector g is then added to the orientation vector u of the considered fiber, resulting in a rotation around the center and in a new orientation vector u'. Each component of g is obtained randomly, following a Gaussian distribution whose variance is proportional to the target angular dispersion. The proposed rotation is accepted if and only if the modified fiber does not collide with another fiber. The global angular dispersion AD is computed as follows:

$\begin{array}{ll}\text{AD}=\frac{1}{{\displaystyle \sum _{P\in {P}}}Card(P)}{\displaystyle \sum _{P\in {P}}}{\displaystyle \sum _{f\in P}}\theta (u_f^{\prime },u_P)& (5)\end{array}$

where $u_f^{\prime }$ is the new orientation vector of fiber f in fiber population P, u_P is the principal orientation vector of fiber population P, $\theta (u_f^{\prime },u_P)$ is the angle between the two lines supported by $u_f^{\prime }$ and u_P vectors, ${P}$ denotes the ensemble of all fiber populations in the considered field of view and Card(P) corresponds to the number of fibers in each fiber population P.

Step 3—The third step of the phantom generation algorithm consists in the induction of local tortuosity in the geometry. To this aim, one population is randomly selected from the set of populations that did not reach their target angular dispersion yet. One fiber of this population and a point along this fiber are then selected randomly. A random orthonormal trieder (x, y, z) is built such that z corresponds to the local direction u_f of the fiber f. The y-axis of the direct trieder defines the direction along which is applied the tortuosity deformation on the fiber (see Figure 1C). This deformation follows a Gaussian distribution with a zero mean and a variance proportional to the tortuosity perturbation value provided as an input. It moves all the points of the fiber in the neighborhood of the selected point that is defined using a tortuosity neighborhood size given as an input parameter of the algorithm for each fiber population. The neighborhood size controls the frequency of the undulations, e.g., the larger this size, the smoother the undulations. Provided that all the deformed points remain in the field of view and that the modified fiber does not collide with any other fiber, the Gaussian deformation is accepted. The induction of local tortuosity is a computationally complex problem. Indeed, the total number of points N_Points in the field of view is given by

$\begin{array}{ll}{N}_{Points}=N_{Populations}.{\displaystyle \sum _{P\in {P}}}{\displaystyle \sum _{f\in P}}Card(f)& (6)\end{array}$

where N_Populations is the number of fiber populations and Card(f) stands for the number of control points of the centroid of a given fiber f. A typical order of magnitude of N_Points is 10⁴ which means that at each step of the induction of local tortuosity, the algorithm must check at 10⁴ points whether the fibers are colliding, which can be computationally heavy. Thus, our algorithm is based on the construction of a look-up table (LUT) which is updated at each step of the tortuosity induction and whose size is optimized so that we check the intersection at each point of the fiber only with all the points in a neighborhood of this very point. The LUT size must be adequately chosen to optimize the computation time. If the LUT size is too small, collisions might be missed whereas a too large LUT size will drastically increase the computation time.

Step 4—Our phantom generation algorithm also enables to distinguish the myelin sheath from the axon. The axon membrane is created within the fiber envelope with a radius computed from a predefined g-ratio, corresponding to the ratio between the axon membrane diameter and the outer fiber diameter. For each population, the g-ratio follows a Gamma distribution, thus adding two further control parameters for each fiber population. The myelin sheath corresponds to the space between the axon membrane and the outer envelope of the fiber (see Figure 1D).

Step 5—The algorithm also accounts for the presence of Ranvier nodes along the myelin sheath (Figure 1E). The internodal distance d is set using Equation (4) (maximal conduction relationship [17]). The width w_R of each Ranvier node corresponds to a fraction α_R of the internodal length d such that w_R = α_R.d, with α_R typically equal to 10⁻³, as described in Salzer [16]. This fraction follows a Gamma distribution adding two further control parameters.

Step 6—Finally, our algorithm gives the possibility to represent beading caused by cytoskeletal damage of the axon membrane: the contours of both axonal and outer myelin membranes are swollen using adequate bell-shaped functions like sine functions (see Figure 1F). The amplitude and spacing of those lobes both follow a Gamma distribution (adding four control parameters for each fiber population).

2.2. Multi-Compartmental Model for Trapezoidal OG Measurements

Oscillating gradient spin echo (OGSE) sequences are sensitive to diffusion on the time scale of the oscillation period rather than the interval between the pulses and can thus enhance the sensitivity to small axonal restrictions. Recently, an ActiveAx OGSE model was proposed in Ginsburger et al. [10] which accounts for the frequency-dependence of the diffusivity transverse to axons in the extra-axonal space using cosine OGSE schemes. A frequency-dependence correction was proposed for the extra-axonal tensor and showed significant reduction of the fit error. However, the correction derived in Ginsburger et al. [10] is only valid for sinusoidal waveforms, while square (or in practice trapezoidal) wave oscillating gradients maximize sensitivity to smaller pore sizes in comparison with sinusoidal sequences. Indeed, they yield the highest diffusion weighting within one period compared to other periodic waveforms [22], which also makes trapezoidal sequences more clinically feasible. The aim of this section is thus to adapt the previous model proposed in Ginsburger et al. [10] to trapezoidal waveforms.

2.2.1. General Model

White matter tissues are modeled using three tissue compartments embedding three types of micro-structural environments: intra-cellular, extra-cellular, and cerebro-spinal fluid (CSF) compartments. A common assumption of effectively impermeable axonal walls is used [3–5]. Thus each compartment provides a separate normalized MR signal and no exchange between the populations of water molecules occurs. The resulting model for the diffusion MR signal S can thus be written as

$\begin{array}{ll}S=(1-{\nu }_{iso})({\nu }_{ic}{S}_{ic}+(1-{\nu }_{ic})S_{ec})+{\nu }_{iso}{S}_{iso}& (7)\end{array}$

where S_ic and ν_ic are the normalized signal and volume fraction of the intra-cellular compartment, S_ec is the normalized signal of the extra-cellular compartment, and S_iso and ν_iso are the normalized signal and volume fraction of the CSF compartment. The model for S_iso assumes an isotropic Gaussian distribution of displacements [23]. Water diffusion in the intra-cellular compartment is restricted by axonal walls and further restricted by myelin sheath in case of myelinated fibers. Fibers are assumed to be parallel and oriented along a single direction n. Hence, one computes the intra-cellular signal S_ic using the Gaussian Phase Distribution approximation of the signal from particles trapped inside a cylinder, which has been derived for trapezoidal OGSE sequences [6, 24]. Diffusion in the extra-axonal compartment is assumed to be hindered. This compartment is generally characterized by a 3D Gaussian displacement distribution:

$\begin{array}{ll}{S}_{ec}=e^{b{\text{G}}^T{D}_{ec}(\text{n},{\nu }_{ic})\text{G}}& (8)\end{array}$

where b is the diffusion sensitization for a given tuning of a diffusion-weighted NMR sequence, G represents the gradient magnitude and D_ec is the extracellular diffusion tensor. Assuming that the diffusion tensor is cylindrically symmetric, D_ec is defined as [5, 25]

$\begin{array}{ll}{D}_{ec}(\text{n},{\nu }_{ic})=(d-d_{\perp }({\nu }_{ic}))\text{n}{\text{n}}^T+d_{\perp }({\nu }_{ic})I& (9)\end{array}$

where d_⊥ is the apparent diffusion coefficient perpendicular to axons and I is the identity tensor [5].

2.2.2. Frequency Dependence of Extra-Axonal Space with Trapezoidal OGSE Sequences

In Ginsburger et al. [10], a correction to the d_⊥ component of the extracellular diffusion tensor D_ec which describes diffusion perpendicular to the fibers was proposed for cosine OGSE, making the diffusion transverse to the fiber bundle in the extra-axonal space d_⊥ dependent on the frequency of the OGSE sequence ω₀ :

$\begin{array}{ll}{d}_{\perp }({\omega }_0)=d_{\perp ,\infty }+A\frac{\pi }{2}\left|{\omega }_0\right|& (10)\end{array}$

where d_{⊥, ∞} is the bulk diffusion constant [5].

Equation (10) was obtained combining two observations. First, the signal attenuation can be written as a function of the dispersive diffusivity ${D}(\omega )$ [26]

$\begin{array}{ll}-lnS~\int \frac{d\omega }{2\pi }{q}_{-\omega }{D}(\omega )q_{\omega }& (11)\end{array}$

where q_ω is the Fourier transform of the integral ${\int }_0^{t}d{t}^{\prime }g(t^{\prime })$ of the applied Larmor frequency gradient g(t). Second, the Fourier transform q_ω of the integral ${\int }_0^{t}d{t}^{\prime }g(t^{\prime })$ of the applied Larmor frequency gradient g(t) can be written as follows

$\begin{array}{ll}{q}_{\omega }=\frac{i\pi \gamma g_0}{{\omega }_0}\left[\delta (\omega -{\omega }_0)-\delta (\omega +{\omega }_0)\right]& (12)\end{array}$

where δ is the Dirac distribution, for an OG profile g(t) = g₀cos(ω₀t) and for a sufficiently large number of oscillations [$\frac{{\omega }_{0}T}{2\pi }\gg 1$, T being the total duration of the gradient train g(t)].

Indeed, the presence of these Dirac distributions in q_ω leads to the much simplified expression of the signal attenuation

$\begin{array}{ll}-lnS(T)~b.Re{D}({\omega }_0)& (13)\end{array}$

which, combined with Equation (2), results in Equation (10).

In the case of cosine trapezoidal OGSE (OGSE-CT), the peak frequency (frequency at which q_ω reaches its maximum) will be the same as for a cosine OGSE (OGSE-C) sequence with the same frequency ω, since the Fourier expansion of the OGSE-CT f_trap(t) is a infinite sum of odd cosine functions

$\begin{array}{ll}{f}_{trap}(\omega t)={\displaystyle \sum _{n=0}^{\infty }}{a}_{2n+1}cos((2n+1)\omega t)& (14)\end{array}$

where |a_2n+1| is a decreasing sequence. However, it is not a priori clear whether the frequency selectivity (which is characterized by both the full-width-half-maximum of the main lobe and the maximum ratio between the side lobes and main lobe amplitudes) will be preserved using OGSE-CT, which might potentially prevent the use of Equation (12) for OGSE-CT. This question was disambiguated in Van et al. [27] who compared the encoding spectrum F(ω) of both sequences. The encoding spectrum for a OGSE-C sequence is:

$\begin{array}{ll}|F(\omega ){|}^2={\left(\frac{\pi \gamma G}{{\omega }_0}\right)}^2\left[\delta (\omega +{\omega }_0)+\delta (\omega -{\omega }_0)\right]& (15)\end{array}$

It was shown that the differences between the encoding spectra of OGSE-CT and OGSE-C waveforms with equal frequency are minimal, thus justifying as well the use of Equation (12) for OGSE-CT sequences having a high selectivity around the frequency ω₀ of interest (see also Figure 2).

FIGURE 2

Figure 2. Power modulation spectra for trapezoidal cosine OGSE gradient waveforms at 64 Hz for different number of half periods (or lobes), showing the influence of the number of lobes on the frequency selectivity. The theoretical peak frequency is denoted as f_th. We observe that the difference between this theoretical peak and the actual frequency peak of the sequence decreases with increasing number of lobes. The actual frequency peak (not the theoretical one) of the OGSE sequences employed in this article is fed in our model for better precision.

Consequently, the same frequency-dependence correction of the extra-axonal space tensor (Equation 10) can be used for both OGSE-C and OGSE-CT waveforms.

2.3. Monte-Carlo Simulations and NMR Signal Synthesis

The Diffusion Microscopist Simulator (DMS) [13] relies on a three-fold architecture:

• a 3D phantom generation engine

• a Monte-Carlo simulator engine

• a DW-NMR signal synthesizer

The 3D phantom generation engine, improved using our novel algorithm to produce biomimicking numerical phantoms, creates the 3D triangulated surfaces from the fiber descriptions. It includes two non-linearly sampled curves (also called axon and outer myelin sheath membrane centroids), and two profiles of radii to represent the axon membrane and the outer myelin sheath membrane, respectively.

The Monte-Carlo simulator engine is composed of several element that have to be individually tuned:

• A scene modeler that contains the description of the evolving simulated space, including the dimensions of the simulation domain corresponding to a global bounding box set to (−60, +60, −60, +60, −60, +60 μm) in all our simulations, the temporal resolution set to 10 μs, the number of simulation steps set to 50,000, the set of membranes generated by the 3D phantom generation engine and the set of particles used to perform the Monte-Carlo simulation (10⁶ particles in our simulations).

• A motion model that drives the motion of particles, set to a Brownian random walk model tuned using a parameter corresponding to the diffusivity of the medium (2.0 × 10⁻⁹ m²/s in our simulations).

• A membrane model built for each individual axon or outer myelin membrane that integrates:

- the triangulated surface itself

- a particle-to-membrane interaction model set to the total reflexion interaction in our case

- a vertex evolution function set to static in our case (DMS is also able to deal with a temporal evolution of the membranes that is not used in this work)

- a polygon cache that facilitates and speeds up the computation of the list of the closest membrane triangles likely to interact with a particle of arbitrary position in the simulation domain

- a particle or random-walker model that represents water molecules moving in the simulation domain. Particles are randomly distributed either over the whole simulation domain, or only in the intra- or extra-cellular space. In this work, two types of simulations were performed, either with particles randomly distributed in the whole domain or restricted to the extra-cellular space (see Figure 3), since the present study is focused on the characterization of the extracellular space signal for different levels of structural disorder.

FIGURE 3

Figure 3. Illustration of the two simulation modalities available in DMS on an example mesh (A) with intracellular volume fraction of 0.2 and 5° angular dispersion induced by tortuosity. In (B) corresponding to a cross-section of (A), random walkers for Monte-Carlo simulation of the diffusion process are placed only in the extra-axonal space, with impermeable membranes. The obtained signal thus stems exclusively from the extracellular space. In (C), particles are initially placed in both intra- and extra-cellular compartments, with impermeable membranes. The diffusion process is thus simulated in both compartments, without exchange between them.

The trajectories of random walkers computed using the Monte-Carlo simulator engine are fed into the diffusion-weighted NMR signal synthesizer, which synthesizes a volume of DW-MRI signals for a given DW-MRI sequence and for a specific tuning of the sequence parameters. The DW-NMR signal synthesizer is also composed of several elements that have to be individually tuned:

• A DW-NMR sequence factory that allows to simulate the chronograms of gradients and radio-frequency pulses; several schemes are available (bipolar double STEAM sequence, bipolar STEAM sequence, multiple PGSE sequence, single PGSE sequence, twice refocused spin echo sequence, OGSE sequence); in this work, a trapezoidal OGSE sequence was tuned (see Figure 4), requiring to define the period of oscillating gradients (varying in this study), the gradient time resolution (10.0 μs), the maximum gradient slew rate and the maximum gradient magnitude (respectively 200 T/m/s and 80 mT/m corresponding to the latest Connectome gradient coils available for 3T MRI systems on the market).

• A Cartesian grid defined within the MC simulation domain using a local bounding box and the 3D volume size; in our case, to avoid boundary effects, the local bounding box was chosen slightly smaller than the global bounding box of the MC domain and set to (−55, +55, −55, +55, −55, +55 μm).

• A noise model to simulate the actual level of noise corrupting the signal of real acquisitions; the analysis of the impact of the noise on the signal of the extra-cellular space and on the inference of the structural white matter disorder has been done for four different values of SNR in this study: an infinite value corresponding to the absence of noise, a SNR of 30 corresponding to state-of-the-art experimental conditions with the latest 3T clinical MRI systems available on the market, and SNR of 20 and 10 corresponding to worse experimental conditions.

• A spin model associated to each random walker in charge of accumulating the net phase evolution induced by the trapezoidal OGSE sequence and seen by the random walker along its trajectory. Although they clearly influence the SNR when long echo times are chosen, the effects of T1/T2 relaxation have not been taken into account by the spin model in the performed simulations.

FIGURE 4

Figure 4. Schematic representation of the employed trapezoidal cosine OGSE gradient waveform, here with six lobes (half-periods) before and after the 180° refocusing pulse. The separation between gradient waveforms, t_sep, is required to accomodate the 180° RF pulse, and has been set to allow a continuous single frequency oscillating gradient to be drawn between the two waveforms to obtain a narrower peak at the desired frequency [28]. The duration of the shorter lobes are increased by half the gradient ramp time, t_ramp, to ensure zero cumulative gradient area.

The synthesized DW-MR image consists of a 3D volume corresponding to the T₂ reference volume at b = 0 s/mm² (set to 10,000 in our case) and a 4D volume corresponding to the employed trapezoidal OGSE sequence along a set of uniformly distributed diffusion directions over the unit sphere (300 directions in our case).

2.4. Estimation of the Scaling Coefficient A with Simulations

Monte-Carlo Simulations were launched by placing 10⁶ particles in the extra-axonal space of each studied white matter numerical phantom. The employed phantoms composed of a single fiber population, listed in Table 2, have an intracellular volume fraction of 0.2 and the fiber radii follow a Gamma distribution whose shape and scale parameters α and β are set to obtain a mean diameter αβ = 2.0 μm and a standard deviation αβ² = 0.2 μm. The signal obtained from diffusion MRI data synthesis comes only from the extra-axonal space, since the fiber membranes were set to be impermeable. Thus, the intra-cellular volume fraction was set to 0 in the employed multi-compartment model. For each studied geometrical configuration, 4D DW-volumes were synthesized for five distinct frequency values of the employed trapezoidal OGSE sequence : 60, 70, 80, 90, and 100 Hz. In order to have a constant b-value of 200 s/mm² and a constant TE of 116ms, the gradient magnitude was varied up to 80 mT/m (corresponding to the maximum achievable gradient magnitude of modern clinical 3T MRI systems) and the number of lobes (half-periods of the OGSE sequence) on each side of the refocusing pulse was also varied, as follows: 6, 8, 8, 10, and 12 lobes/50, 57, 70, 75, and 80 mT/m for encoding frequencies of 60, 70, 80, 90, and 100 Hz. Figure 4 gives an illustration of the employed trapezoidal OGSE sequence. The multiple-frequency sampling enables to perform a linear fitting procedure on Equation (10) in order to obtain reliable estimates of the scaling coefficient A of the extra-axonal perpendicular diffusivity linear-in-frequency term. The error associated to the linear fit was assessed using the covariance matrix returned by the linear least-square procedure employed to estimate the scaling coefficient A in Equation (10). These fitting errors are represented by error bars in figures giving estimated values of A for various geometric configurations (see section 3).

TABLE 2

Table 2. Values of the control parameters of the phantom generation algorithm for each studied configuration (C1 to C5).

The evolution of the scaling coefficient A with respect to the various sources of disorder is explored by simulating the diffusion process within the geometric configurations presented in Table 2 (configurations C1–C5, as well as some variants of those configurations). Table 2 provides the list of all the control parameters as well as their values employed for the construction of the numerical phantoms. The purpose is to catch the complex dependence of A with respect to geometrical characteristics of the diffusion medium. The degree of complexity of the employed numerical phantom is fully monitored using the DMS simulator, thus enabling a very precise study of the influence of each parameter on the fitted coefficient A.

3. Results

3.1. Numerical Phantoms Mimicking White Matter

We present in this section numerical phantoms generated with our novel algorithm. Figure 5 represents various configurations with one, two, and three fiber populations and with or without global angular dispersion. Figure 6-1 represents a set of straight parallel fibers randomly placed in the phantom volume, with mean fiber diameter of 2.0 μm (diameter variance of 0.2 μm) and volume fraction of 0.1. In Figure 6-2, global angular dispersion is induced, enabling to reach 5.6° of angular dispersion (for a target value of 10°). The tortuosity induction (Figure 6-3) brings this angular dispersion up to the 10° target. Figures 6-4,5 illustrate the creation of myelin sheath and Ranvier nodes which account for the actual structure of myelinated fibers. Beading -consisting in a swelling of both axonal and myelin sheath membranes—is also handled (Figure 6-6). In all presented surface renderings, there is no collision between the membranes. Figure 6-6 presents a realistic geometry mimicking a complex white matter environment, taking into account all the putative deformations of membranes observed in real tissues (angular dispersion, tortuosity, myelination and creation of Ranvier nodes, beadings). We note that the generated phantoms have been presented here in the case of multiple fiber populations for more generality. However, the simulations performed to study the structural disorder of white matter in next sections only used phantoms with a single fiber population. While essential to make our approach directly applicable to most white matter regions, the study of the structural disorder induced in crossing areas of multiple fiber populations is out of the scope of this study. This topic is addressed in the section 4.

FIGURE 5

Figure 5. (Top) Surface renderings of outer fiber envelopes corresponding to one (left), two (middle), and three (right) fiber populations with intracellular fraction of 0.06 for one population, 0.03 each for the two populations, 0.02 each for the three populations and mean diameters 2.5 μm (one population), 1.5 and 2.5 μm (two populations), 1.5, 2.0, and 2.5 μm (three populations). Bottom same as top. with global angular dispersion of 10° for the single population, 5.5 and 8.5° for the two populations, 5.5, 6.7, and 8.5° for the three populations (for a target angular dispersion of 10° each).

FIGURE 6

Figure 6. (1) Surface renderings corresponding to three fiber populations with mean diameters 1.5, 2.0, and 3.0 μm, respectively with intracellular volume fraction of 0.05 each. N is the number of populations in general (here N = 3). (2) Global angular dispersion is induced (4.7, 2.5, and 2.8° per population). (3) Tortuosity is induced, enabling to reach 10° of angular dispersion for each population. (4) Creation of the myelin sheath. (5) Creation of Ranvier nodes with mean ratio 1,000 between internodal length and node width. (6) Creation of beadings with 20.0 μm mean inter-beading length and beading magnitude ratios of 1.7, 1.5, and 1.2, respectively.

3.2. Measuring Structural Disorder in the Extra-Axonal Space

3.2.1. Validation of the Employed Trapezoidal OGSE Model

A linear relationship between the extra-axonal perpendicular diffusivity and the frequency of the employed trapezoidal OGSE sequence is shown in Figure 8, where the frequency was varied from 60 to 100 Hz at a constant b-value of 200 s/mm². This result validates the use of Equation (10).

3.2.2. Studying Different Types of Structural Disorders

Monte-Carlo simulations were run in geometries characterized by increasing structural complexity. In what follows, the different geometrical configurations were classified according to the type of structural disorder they represent, from configurations C1 to C5. Table 2 lists all the employed geometrical configurations and gives the associated set of control parameters.

In configuration C1 (Figure 7), parallel axons are randomly placed in the simulation volume. The centers of the cross section of the cylinders representing the axons follow a uniform distribution. The diameters follow a Gamma distribution, with a mean axonal diameter of 2.0 μm. The value of the scaling coefficient A reaches its maximum value of 9.09 μm² (all the values of A are given for a SNR of 30) in this configuration which corresponds exactly to the 2D short-range disorder geometry described in Burcaw et al. [7] and Novikov et al. [11]. From C1 to C2, the induction of global angular dispersion (fibers remain straight but are rotated to induce angular dispersion) of 3.5° yields to a substantial diminution of the scaling coefficient from 9.09 down to 8.54 μm². Local tortuosity is induced in configuration C3 (fibers are deformed to induce angular dispersion), enabling to reach an angular dispersion value of 15° and yielding a moderate decrease of the scaling coefficient from 8.54 down to 8.26 μm². The addition of Ranvier nodes along the myelin sheath (with a g-ratio of 0.6) in configuration C4 does not significantly change the value of the estimated A coefficient with respect to the previous configuration (A = 8.24 μm² vs. A = 8.26 μm² in configuration C3, which is negligible given the fitting uncertainty at a SNR of 30). From C4 to C5, beading is introduced (swelling of the axonal membrane and myelin sheath) with amplitudes equal to 1.5 times the fiber radii in average, and mean spacing of 20.0 μm, yielding a significant decrease of A (up to 19% with A going from 8.32 to 6.76 μm²).

FIGURE 7

Figure 7. From left to right, simulation domains corresponding to different geometrical configurations (from C1 to C5) are shown. A 3D rendering of the root-mean-squared difference (represented by the minus sign) between the diffusion signal stemming from each configuration and the “reference” configuration C1 (most left configuration composed of straight parallel cylinders) is shown on a spherical surface. (A) The RMS signal difference is computed between configuration C1 and C2. The red area where the diffusion signal difference is the strongest corresponds to diffusion sensitization directions perpendicular to fibers. Blue areas correspond to directions parallel to fibers where the signal differences are weaker but not null, and originate from the variations of diffusion properties around those directions when structural disorder is added. (B–D) Represent the same RMS signal differences between configurations C1 and C3, C1 and C4_R4 (corresponding to configuration C4 with a demyelination ratio of 25%, see Table 2), C1 and C5, respectively.

It appears from Figures 7A–D. that the diffusion signal variation between configuration C1 and configurations C2–C5 is stronger for diffusion sensitization directions perpendicular to fibers, meaning that most of the information stemming from the increasing complexity of the studied geometrical configurations is included in the evolution of perpendicular diffusivity. Figure 7 also clearly indicates that beading (corresponding to configuration C5) has the strongest effect on the perpendicular diffusivity. This result can be directly related to the important variations of A observed in configuration C5, as reported previously.

3.3. Understanding the Quantitative Influence of Each Disorder Parameter

The different configurations studied in Figure 9 correspond to qualitatively different geometric configurations. The quantitative evolution of the structural disorder coefficient A with respect to each geometrical parameter has also been studied (see Figure 10), by varying one geometrical parameter in each configuration separately. The evolution of each parameter is explained in Table 2.

The induction of global angular dispersion from 0° to 4.5° in configuration C2 (see Table 2) results in a decrease of the scaling coefficient A (see Figure 10A) from 9.09 to 8.43 μm² for the biggest value of global angular dispersion of 4.5°. The diminution of A gets stronger for increasing values of global angular dispersion. However, the study of the influence of global angular dispersion on A was limited to small values of angular dispersion, owing to the fact that higher values of angular dispersion were not reachable for the specific fiber density and radii distribution of C2. These values of angular dispersion are far from the values of microscopic misalignments of axons estimated up to 18° [29]. One possibility to reach this target is to decrease the intracellular volume fraction. Another option consists in increasing angular dispersion using local tortuosity, as was done in configuration C3. The induction of tortuosity in this configuration causes a moderate decrease of the scaling coefficient from 8.54 μm² corresponding to configurations C2 (global angular dispersion of 3.5°) to 8.24 μm² for the biggest tortuosity value (20° of angular dispersion), as shown in Figure 10B.

In configuration C4 (corresponding to a global angular dispersion of 3.5° and to a total angular dispersion of 15° after the induction of tortuosity, with the presence of Ranvier nodes), a small increase of the scaling coefficient is observed due to demyelination, from 8.24 up to 8.32 μm² for a demyelination ratio of 25% (see Figure 10C), which is not significant given the uncertainties of the estimated values of A at a SNR of 30.

The introduction of beading in configuration C5 causes important variations of the scaling coefficient A (Figure 10D). As mentioned previously, a mean beading spacing of 20.0 μm yields a significant decrease of A (up to 19% with A going from 8.32 to 6.76 μm² The scaling coefficient A = 7.95 μm² is still 4.5% smaller for a beading spacing of 100.0 μm, corresponding to a low density of beads, than in the absence of beading where A = 8.32 μm².

4. Discussion

4.1. White Matter Biomimicking Numerical Phantoms

Designing realistic numerical phantoms of white matter tissue seems to be a promising approach to study the influence of various structural properties of white matter on the measured diffusion signal. Being able to construct biomimicking simulations geometries without using histological samples is an essential step toward the comprehension of the specific effect of each geometrical characteristic of the diffusing medium on the obtained signal. An important aspect of our phantom generation algorithm is its ability to deal with multiple fiber populations, which was reported to represent up to 60% of the number of voxels of a mask of the brain white matter at a spatial resolution of 2 mm [15]. However, generating phantoms with multiple fiber populations come with additional difficulties, notably related to the generation of global angular dispersion. Phantoms presented in Figure 5 exhibit geometrical configurations with multiple fiber populations. In Figure 5, the target angular dispersion value of 10° can be reached only in the single population case for an intra-cellular fraction of 0.2. The maximum reachable global angular dispersion strongly depends on the number of fiber populations, on the distribution of radii and on the target intra-cellular volume fractions of these populations (for instance, the lower the intra-cellular fraction, the higher the reachable angular dispersion). The use of tortuosity is essential to reach high values of angular dispersion in multiple population configurations. Indeed, in Figure 6-2, the induction of global dispersion enables to reach 5.6° of angular dispersion (for a target of 10°). The tortuosity induction (Figure 6-3) brings this angular dispersion up to the 10° target.

4.2. Characterizing Structural Disorder with A

As shown in Figure 7, the differences between diffusion signals obtained by simulating the diffusion process in the extra-axonal space of various geometrical configurations are predominant in the diffusion direction perpendicular to fibers, although there exists differences in all directions. This observation shows that most of the structural disorder effects are caught by the diffusion signal around the equator perpendicular to the mean fiber orientation. Modeling structural disorder using an additional term in the extra-axonal perpendicular diffusivity thus appears to be physically reasonable.

Across all configurations, at a SNR of 30, the values of the scaling coefficient A vary between 9.09 and 6.77 μm² which is consistent with previous studies [7, 8] reporting the empirical law $A~l_c^2$ were the correlation length l_c closely follows the mean external radius of the fibers. Indeed, in our case the mean axonal diameter is equal to 2.0 μm yielding l_c ~ 1.0 μm, which leads to A ~ 1.0 μm². This value corresponds to the order of magnitude of the fitted values of A.

As reported in Fieremans et al. [8], the effect of angular dispersion introduced in configurations C2 and C3 on the estimated perpendicular diffusivity can be understood by considering the diffusion process along each fiber in the presence of orientational dispersion. When the orientation of a given fiber differs from the mean orientation of the fibers population, the diffusion process along this fiber yields a local parallel diffusivity whose projection on the plane perpendicular to the mean fiber population direction is not null. Due to this projection effect, the existence of longitudinal frequency-dependence along each elementary fiber will yield a frequency-dependence of the global perpendicular diffusivity. Thus, the frequency-dependence observed in our work might not only originate from the 2D short-range disorder in the plane transverse to axons, but might also be partly explained by the contamination of the perpendicular diffusivity (and thus of the scaling coefficient A) with longitudinal diffusion frequency-dependence. As a consequence, in addition to modeling the 2D short-range disorder in the plane perpendicular to fibers, the estimated scaling coefficient A might also embed some information related to physical properties along the fibers. While this hypothesis makes the physical interpretation of the coefficient A tricky, it would enable to obtain information from both parallel and perpendicular extra-axonal diffusion processes by measuring only one coefficient. The interest of such an approach is that the estimation of A relies on a simple and robust fitting procedure. It only requires to perform data acquisition using trapezoidal OGSE sequences with a sufficient number of frequencies to be able to fit the data properly, which is clinically feasible.

The introduction of Ranvier nodes in configuration C4 does not yield to a significant change in the scaling coefficient A. This is not surprising since Ranvier nodes correspond to a very low volume fraction of the extracellular space, due to their low width, reported to correspond on average to a few thousandths of the internodal length [16], which is amplified by the low intracellular fraction of 0.2 employed in our simulations. In Figure 10C, demyelination is mimicked by progressively increasing the width of the Ranvier nodes, from 1 to 25% of the internodal length. The observed increase of the scaling coefficient A is again small (up to 3% for the maximal demyelination ratio of 25%, roughly corresponding to the percentage of demyelinated areas observed in the cerebral cortex of multiple sclerosis patients [30]).

In this study, only the structural effects of Ranvier nodes and demyelination on diffusion properties in the extra-axonal space were studied. However, from the diffusion point of view, the most interesting feature of Ranvier nodes and unmyelinated regions areas is that they represent those regions along a myelinated axon where the exchange between intra- and extra-axonal water is the fastest. The analysis of the effect of such an increased exchange would be of great interest to thoroughly study the effect of Ranvier nodes and demyelination on the scaling constant A. Indeed, in the Appendix F of Burcaw et al. [7], the authors theoretically tackled this problem and their prediction is that the values of the scaling constant A may or may not be affected by the exchange, depending on the exchange regime (slow, intermediate, or fast). Their theoretical analysis takes into account a general uniform exchange of molecules between intra- and extra-axonal space. However, the Ranvier nodes and demyelination around them would introduce a local exchange linked to the disorder with which the Ranvier nodes occurs within the voxel. In this condition, it is not clear if the argument in Burcaw et al. [7] still holds. This question will be investigated in a future work.

4.3. The Strong Influence of Beading on the Scaling Coefficient A

The value of the scaling coefficient A appears to be mostly influenced by the local enlargement of both axonal and myelin membranes, also called beading. Indeed, Figure 9 shows a significant decrease of A (up to 19%) in the presence of beading.

According to Figure 10D the scaling coefficient is still 4.5% smaller for a beading spacing of 100.0 μm than in the absence of beading, suggesting that A could be a putative marker of the presence of beading, since A is significantly reduced even for a lower density of beadings within the phantom.

The influence of beading on the diffusion signal has been extensively studied in Budde and Frank [20], where it was emphasized that neurite beading might explain the decrease of the apparent diffusion coefficient after ischemic stroke. Our results point in the same direction, since the significant decrease of the coefficient A in case of beading induces a net decrease of the transverse diffusivity which results in a diminution of the apparent diffusion coefficient. However, in Budde and Frank [20], simulations were run on numerical phantoms with an hexagonal packing of fibers and periodic restrictions along the fiber. The periodicity of the employed simulation domain might affect the realism of the obtained diffusion signal since it does not reflect the 1D short-range disorder along white matter fibers nor the 2D short-range disorder in the plane transverse to fibers [7, 11]. In our simulations, short-range disorder effects are expected to be accounted for, since the spacing between each beading is highly variable (as shown in Table 2, the variance of the spacing distribution is equal to one fourth of the spacing value) and the fibers are randomly placed in the phantom, thus mimicking transverse short-range disorder better representing actual brain white matter tissues.

The extracellular diffusion signal obtained from our simulations thus gives a more realistic view of the effect of beading on the transverse diffusion coefficient, which appears to be quantitatively significant. As discussed earlier, the effect of beading on the scaling coefficient A might partly originate from the hindrance of extra-axonal diffusion along the beaded fibers, due to the projection of parallel diffusivity on the transverse plane, as a consequence of angular dispersion.

The effect of beadings on the intra-cellular diffusion process -which is not studied in this work- might also be quantitatively substantial, as suggested in Marco et al. [31]. In this work, the sensitivity of the diffusion signal of intracellular metabolites with respect to beaded structures was studied using Monte-Carlo simulation of brain metabolites dynamics, which can be compared, from the numerical simulation point of view, with the waters one after proper scaling. A clear dependence of both radial and axial intracellular dispersive diffusivities with respect to the frequency of the employed OGSE sequence in the presence of beadings was observed. Moreover, results of MC molecular diffusion simulations in complex synthetic substrates mimicking the presence of beads showed a clear $1/\sqrt{t}$ dependence of the axial intracellular Apparent Diffusion Coefficient due to 1D short-range disorder introduced in the axial direction by the randomly placed beads, in good accordance with theoretical predictions in Burcaw et al. [7] and Novikov et al. in [11]. Indeed, beadings are supposed to primarily affect the diffusion process along the fibers and a $1/\sqrt{t}~\sqrt{\omega }$ frequency-dependence of the parallel diffusivity in the intra-cellular space of beaded axons is expected. Similarly to what was done in Marco et al. [31] for intracellular metabolites, an interesting approach would be to capture directly the effect of beading on the intra-cellular parallel diffusivity $\sqrt{\omega }$ term, by performing 3D Monte-Carlo simulations of spin dynamics in the intra-cellular space and fitting the scaling factor of the $\sqrt{\omega }$ term for various beaded geometrical configurations. The amount of variation of this “intra-cellular disorder” scaling factor in the presence of beading could be compared to the variation of the scaling factor A studied in this work. Moreover, since it has been observed in this study that A depends (mostly) both on angular dispersion and beading, further information from an intracellular model accounting for beading would enable to disentangle the influence of angular dispersion and beading on the coefficient A.

In any case, the presented approach provides a reliable way to detect beading-induced modifications of the diffusion process, which are measurable for a SNR greater or equal to 20, as shown in Figure 9.

The estimation of A was performed using a linear fit with five values of the perpendicular diffusivity corresponding to five distinct OGSE frequencies. The degradation of the estimation of A while reducing the number of measurements (by using only three OGSE frequencies for instance) should also be studied in order to reduce acquisition time. However, a protocol relying on a single b-value at five different OGSE frequencies and along 60 directions already meets the requirements of a clinical research protocol and will be used in the future to assess all the findings presented in this work. The range of explored OGSE frequencies from 60 to 100 Hz allows to obtain a sufficiently short echo time (<120 ms) that enables to preserve at least 23% of the magnetization before diffusion-weighting, when considering an average T₂ value of 80 ms at 3T. A diffusion sensitization of 200 s/mm² still preserves 20% of the signal after T₂ relaxation and diffusion decay for a diffusion coefficient of 0.7 × 10⁻⁹ m²/s. Therefore, it seems possible, if the voxel spatial resolution is kept on the order of 2mm, to apply this imaging protocol on a clinical 3T MRI system in vivo in human subjects.

4.4. Effect of SNR

In order to invoke practical conclusions from the numerical simulation results reported here, the present study addressed the impact of noise on the quality of the A scaling coefficient fit. Indeed, Gaussian noises with equal standard deviations were added to the real and imaginary parts of the complex NMR signal before computing its magnitude which corresponds to the simulated diffusion-weighted signal, resulting in a Rician noise corruption, as expected. Three different values of SNR were employed: SNR = 30 (corresponding to good experimental conditions with the latest 3T clinical MRI systems available on the market), SNR = 20 (intermediate experimental conditions), and SNR = 10 (worse experimental conditions). As shown in Figure 8, our noise analysis suggests that a SNR of 10 does not enable to fit the scaling coefficient A properly since the estimation of the perpendicular diffusivity for each frequency value has a too large uncertainty: the differences between the obtained values of A for the different geometrical configurations (shown in Figures 9, 10) are strongly mitigated by noise. However, this analysis suggests that it is possible to reliably detect changes in perpendicular diffusivity and estimate the corresponding scaling coefficient A for a SNR greater or equal to 20 (see again Figures 8–10). The algorithm employed to estimate the parameters of the model presented in this work is part of a framework which maximizes a Rician log-likelihood function using a robust Expectation Maximization algorithm. In the presence of noise in DW data, the use of such a framework enables to alleviate the fit error of the model parameters.

FIGURE 8

Figure 8. Frequency-dependent perpendicular diffusivity in the extracellular space measured by performing Monte-Carlo simulations with diffusing particles in the extracellular space of configuration C1 (see Table 2), plotted against the frequency of the employed OGSE-CT sequence. A linear fit is also plotted which shows the linear dependence of diffusivity to frequency.

FIGURE 9

Figure 9. Value of the scaling coefficient A in m² for geometric configurations with increasing structural disorder. The simulation are performed in configurations C1–C5, whose design parameters are summarized in Table 2.

FIGURE 10

Figure 10. (A) Value of the scaling coefficient A in m² plotted against the global angular dispersion. The simulations were performed in variants of configuration C2 by varying the global angular dispersion value (see Table 2). (B) Value of the scaling coefficient A plotted against the angular dispersion. The simulations were performed in variants of configuration C3 by varying the local tortuosity value (see Table 2). (C) Value of the scaling coefficient A plotted against the percentage of demyelination. The simulations were performed in variants of configuration C4 by varying the percentage of demyelination, directly related to the width of Ranvier nodes (see Table 2). (D) Value of the scaling coefficient A plotted against the mean beading spacing. The simulations were performed in variants of configuration C5 by varying the value of the mean beading spacing (see Table 2). A bigger spacing yields to a lower beading density.

4.5. Limitations

An intracellular fraction of 0.2 was employed to generate our numerical phantoms, which is not realistic since values of intracellular fraction are reported in the range (0.6–0.8). The choice of such a low fraction was deliberately made because it enables an important variation of both global and local angular dispersion as well as beading amplitude which is not possible at higher intracellular fractions. Those geometrical considerations are a major difficulty when trying to design realistic numerical phantoms, which becomes even stronger in the case of multiple populations.

The simulations performed in this article only used phantoms with a single fiber population, thus omitting the effect of crossing fibers on the scaling coefficient A. Studying the effect of multiple fiber populations on our model is essential and will be possible since our algorithm to design numerical phantoms enables to generate multiple populations geometries. However, we chose to restrict this study to single populations because (1) this study has never been performed, even in the case of a single population; understanding the various structural disorder effects for one population is, in our opinion, already a major challenge (2) it enables a more important variation of both global and local angular dispersion whose effects are of interest and (3) choosing a proper adaptation of our model to multiple fiber configurations requires a thorough investigation which should be the object of a complete study. Indeed, there are at least two distinct approaches to adapt our model to fiber crossing configurations. The first approach relies on the hypothesis formulated in Burcaw et al. [7] that the behavior of Equation (2) will persist in fiber crossing regions because as long as the neurites of each tract are randomly positioned, the dynamical exponent will remain equal to 1, thus still leading to a |ω| frequency dependence in the extra-axonal space. However, in the case of fiber crossing, Equation (2) will no longer describe the frequency-dependence of the sole perpendicular diffusivity. Indeed, the fact of introducing fibers with multiple directions leads to a 3D disorder, while in the case of parallel fibers, the disorder was 2D, in the plane perpendicular to fibers. Thus, Equation (2) will apply not only to the diffusion coefficient transverse to axons—the perpendicular diffusivity—but to the overall diffusion coefficient. The second approach assumes that there is not a qualitative change of the underlying physics when introducing fiber crossings. In this case, an adaptation of our model to fiber crossings would draw from similar adaptations of state-of-the art microstructure models to deal with multiple fiber configurations, such as AMICOx [32] which estimates axon diameter indices in two fiber orientations (synthetic data only, using ActiveAx model in two orientations [5]). A second approach [33] introduced the spherical mean technique, capable of factoring out the effects of fiber crossing to estimate per-axon parallel and perpendicular effective diffusion coefficients, and subsequently extract fiber orientation using spherical deconvolution. Similarly, estimation of NODDI in two directions [34] for tractography uses fiber orientation estimates from neighboring voxels.

5. Conclusion

In this article, a novel tool to design more realistic phantoms of white matter was presented, enabling to study the influence of different geometrical features on the linear-in-frequency dependence of the extra-axonal perpendicular diffusivity, weighted by a scaling coefficient A.

By performing Monte-Carlo simulations in the extracellular space of numerical phantoms with increasing geometrical complexity, it was observed that this scaling coefficient A is sensitive to the modification of geometrical properties of the diffusing medium, such as the introduction of global angular dispersion and tortuosity. The presence of Ranvier nodes and demyelinated areas along the axons in the numerical phantom did not seem to significantly change the fitted value of A and further simulations in both intra- and extra-axonal spaces taking into account the high level of exchange around unmyelinated areas have to be performed to possibly observe as stronger effect on the coefficient A. The introduction of beading in the numerical phantoms was by far the most impacting geometrical modification, with a strong deviation of the fitted scaling coefficient A from geometries without beaded structures.

Future work will consist in studying the effect of multiple fibers populations on the estimation of the scaling coefficient A, since crossing configurations represent at least 60% of white matter regions. The effect of further geometrical characteristics on the structural disorder coefficient A, such as the presence of astrocytes and oligodendrocytes which could slow down the diffusion in the extra-axonal space, should also be considered by introducing those geometries in our numerical phantom generation algorithm. Further developments are also needed to be able to reach higher values of angular dispersion and beading at high and realistic intracellular volume fractions (>0.7).

This simulation study shows the importance of the generation of more realistic numerical phantoms in order to catch the complexity of the underlying diffusion biophysics. Analytical models such as the one employed in this article enable to assess the degree of realism needed to perform Monte-Carlo simulations reflecting the actual diffusion process in white matter without adding dispensable and computationally costly details in the phantoms geometry. This is a necessary step towards the construction of dictionaries of simulated biomimicking geometries to inversely decode white matter microstructure.

Author Contributions

KG, CP, J-FM, and MA: Designed research; KG, JB, and CP: Performed analytical calculations; KG, FP, DE, and FM: Performed numerical calculations; KG, JB, and CP: Analyzed data; KG and CP: Wrote the paper.

Funding

This project has received funding from the European Union's Horizon 2020 Framework Programme for Research and Innovation under Grant Agreement No 720270 (Human Brain Project SGA1).

Conflict of Interest Statement

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

The handling editor declared a shared affiliation, though no other collaboration, with the authors KG, FP, JB, DE, JF-M, CP.

References

1. 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.

Google Scholar

2. Portnoy S, Flint J, Blackband S, Stanisz G. Oscillating and pulsed gradient diffusion magnetic resonance microscopy over an extended b-value range: implications for the characterization of tissue microstructure. Magn Reson Med. (2013) 69:1131–45. doi: 10.1002/mrm.24325

PubMed Abstract | CrossRef Full Text | Google Scholar

3. 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

4. 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

5. 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

6. Ianuş A, Siow B, Drobnjak I, Zhang H, Alexander DC. Gaussian phase distribution approximations for oscillating gradient spin echo diffusion MRI. J Magn Reson. (2013) 227:25–34. doi: 10.1016/j.jmr.2012.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Burcaw LM, Fieremans E, Novikov DS. Mesoscopic structure of neuronal tracts from time-dependent diffusion. Neuroimage (2015) 114:18–37. doi: 10.1016/j.neuroimage.2015.03.061

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Fieremans E, Burcaw LM, Lee HH, Lemberskiy G, Veraart J, Novikov DS. In vivo observation and biophysical interpretation of time-dependent diffusion in human white matter. Neuroimage (2016) 129:414–27. doi: 10.1016/j.neuroimage.2016.01.018

PubMed Abstract | CrossRef Full Text | Google Scholar

9. De Santis S, Jones DK, Roebroeck A. Including diffusion time dependence in the extra-axonal space improves in vivo estimates of axonal diameter and density in human white matter. Neuroimage (2016) 130:91–103. doi: 10.1016/j.neuroimage.2016.01.047

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ginsburger K, Poupon F, Mangin JF, Poupon C. “Frequency dependence of the extra-axonal diffusion coefficient for OGSE sequences,” in Proceedings of ESMRMB 2017. Barcelona (2017).

11. Novikov DS, Jensen JH, Helpern JA, Fieremans E. Revealing mesoscopic structural universality with diffusion. Proc Natl Acad Sci USA. (2014) 111:5088–93. doi: 10.1073/pnas.1316944111

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Hall MG, Alexander DC. Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI. IEEE Trans Med Imaging (2009) 28:1354–64. doi: 10.1109/TMI.2009.2015756

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Yeh CH, Schmitt B, Le Bihan D, Li-Schlittgen JR, Lin CP, Poupon C. Diffusion microscopist simulator: a general Monte Carlo simulation system for diffusion magnetic resonance imaging. PLoS ONE (2013) 8:e76626. doi: 10.1371/journal.pone.0076626

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Neher PF, Laun FB, Stieltjes B, Maier-Hein KH. Fiberfox: facilitating the creation of realistic white matter software phantoms. Magn Reson Med. (2014) 72:1460–70. doi: 10.1002/mrm.25045

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Behrens TE, Berg HJ, Jbabdi S, Rushworth MF, 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

16. Salzer JL. Clustering sodium channels at the node of Ranvier: close encounters of the axon-glia kind. Neuron (1997) 18:843–46.

PubMed Abstract | Google Scholar

17. Rushton WAH. A theory of the effects of fibre size in medullated nerve. J Physiol. (1951) 115:101–22.

PubMed Abstract | Google Scholar

18. Beaulieu C. The basis of anisotropic water diffusion in the nervous system-a technical review. NMR Biomed. (2002) 15:435–55. doi: 10.1002/nbm.782

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Daducci A, Dal Palú A, Descoteaux M, Thiran JP. Microstructure informed tractography: pitfalls and open challenges. Front Neurosci. (2016) 10:247. doi: 10.3389/fnins.2016.00247

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Budde MD, Frank JA. Neurite beading is sufficient to decrease the apparent diffusion coefficient after ischemic stroke. Proc Natl Acad Sci USA. (2010) 107:14472–7. doi: 10.1073/pnas.1004841107

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wedeen VJ, Rosene DL, Wang R, Dai G, Mortazavi F, Hagmann P, et al. The geometric structure of the brain fiber pathways. Science (2012) 335:1628–34. doi: 10.1126/science.1215280

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Drobnjak I, Zhang H, Ianuş A, Kaden E, Alexander DC. PGSE, OGSE, and sensitivity to axon diameter in diffusion MRI: insight from a simulation study. Magn Reson Med. (2016) 75:688–700. doi: 10.1002/mrm.25631

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Barazany D, Basser PJ, Assaf Y. In vivo measurement of axon diameter distribution in the corpus callosum of rat brain. Brain (2009) 132:1210–20. doi: 10.1093/brain/awp042

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Xu J, Does MD, Gore JC. Sensitivity of MR diffusion measurements to variations in intracellular structure: effects of nuclear size. Magn Reson Med. (2009) 61:828–33. doi: 10.1002/mrm.21793

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Jespersen SN, Kroenke CD, Østergaard L, Ackerman JJ, Yablonskiy DA. Modeling dendrite density from magnetic resonance diffusion measurements. Neuroimage (2007) 34:1473–86. doi: 10.1016/j.neuroimage.2006.10.037

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Novikov DS, Kiselev VG. Surface-to-volume ratio with oscillating gradients. J Magn Reson. (2011) 210:141–5. doi: 10.1016/j.jmr.2011.02.011

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Van AT, Holdsworth SJ, Bammer R. In vivo investigation of restricted diffusion in the human brain with optimized oscillating diffusion gradient encoding. Magn Reson Med. (2014) 71:83–94. doi: 10.1002/mrm.24632

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Baron CA, Beaulieu C. Oscillating gradient spin-echo (OGSE) diffusion tensor imaging of the human brain. Magn Reson Med. (2014) 72:726–36. doi: 10.1002/mrm.24987

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Ronen I, Budde M, Ercan E, Annese J, Techawiboonwong A, Webb A. Microstructural organization of axons in the human corpus callosum quantified by diffusion-weighted magnetic resonance spectroscopy of N-acetylaspartate and post-mortem histology. Brain Struct Funct. (2014) 219:1773–85. doi: 10.1007/s00429-013-0600-0

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Bø L, Vedeler CA, Nyland HI, Trapp BD, Mørk SJ. Subpial demyelination in the cerebral cortex of multiple sclerosis patients. J Neuropathol Exp Neurol. (2003) 62:723–32. doi: 10.1093/jnen/62.7.723

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Palombo M, Ligneul C, Valette J, Hernández-Garzón E. Can we detect the effect of spines and leaflets on the diffusion of brain intracellular metabolites? Neuroimage (2017). doi: 10.1016/j.neuroimage.2017.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Auría A, Romascano D, Canales-Rodriguen E, Wiaux Y, Dirby T, Alexander D, et al. “Accelerated microstructure imaging via convex optimisation for regions with multiple fibres (AMICOx),” in 2015 IEEE International Conference on Image Processing (ICIP). Quebec City, QC (2015) p. 1673–76.

33. 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

34. Reddy CP, Rathi Y. Joint multi-fiber NODDI parameter estimation and tractography using the unscented information filter. Front Neurosci. (2016) 10:166. doi: 10.3389/fnins.2016.00166

PubMed Abstract | CrossRef Full Text | Google Scholar