Robust Monte-Carlo Simulations in Diffusion-MRI: Effect of the Substrate Complexity and Parameter Choice on the Reproducibility of Results

Monte-Carlo Diffusion Simulations (MCDS) have been used extensively as a ground truth tool for the validation of microstructure models for Diffusion-Weighted MRI. However, methodological pitfalls in the design of the biomimicking geometrical configurations and the simulation parameters can lead to approximation biases. Such pitfalls affect the reliability of the estimated signal, as well as its validity and reproducibility as ground truth data. In this work, we first present a set of experiments in order to study three critical pitfalls encountered in the design of MCDS in the literature, namely, the number of simulated particles and time steps, simplifications in the intra-axonal substrate representation, and the impact of the substrate's size on the signal stemming from the extra-axonal space. The results obtained show important changes in the simulated signals and the recovered microstructure features when changes in those parameters are introduced. Thereupon, driven by our findings from the first studies, we outline a general framework able to generate complex substrates. We show the framework's capability to overcome the aforementioned simplifications by generating a complex crossing substrate, which preserves the volume in the crossing area and achieves a high packing density. The results presented in this work, along with the simulator developed, pave the way toward more realistic and reproducible Monte-Carlo simulations for Diffusion-Weighted MRI.

1. Introduction

Diffusion-Weighted Magnetic Resonance Imaging (DW-MRI) is a non-invasive technique with enormous potential for the study of the brain's microstructure by measuring the diffusion properties of biological tissue. For instance, state-of-the-art methods can use these measurements to estimate tissue properties of the brain white matter, e.g., axonal diameter estimations (Assaf et al., 2008), orientation and volume fraction of the axonal bundles (Zhang et al., 2011; Daducci et al., 2015b) and neurite dispersion (Zhang et al., 2012). The previous information is valuable to understand the brain's maturation (Nilsson et al., 2012; Sexton et al., 2014) as well as the degeneration process associated with neuronal diseases like axonal degeneration (Lovas et al., 2000) and multiple sclerosis (Trapp et al., 1998).

In DW-MRI, an attenuated signal is usually recovered via a Pulsed Gradient Spin-Echo protocol (PGSE) (Stejskal and Tanner, 1965) sensitive to the displacement of the water molecules. Analytical solutions of the signal attenuation can be derived for simple geometrical shapes such as impermeable planes, cylinders, and spheres (Neuman, 1974). However, for applications where the signal attenuation of complex cellular structures or non-homogeneous media is needed, e.g., to generate ground truth data, an analytical solution is no longer feasible to pursue due to its inherent complexity. Because of this, simplifications of the diffusion media have been used as the backbone of most of the microstructure models in the literature (Bammer, 2003; Panagiotaki et al., 2012; Zhang et al., 2012; Ferizi et al., 2015).

Monte-Carlo Diffusion Simulations (MCDS) provide a fundamental approach to study the diffusion phenomena in scenarios where the analytical solutions cannot be computed due to their complexity. In contrast with other numerical methods, MCDS does not require an explicit model of the diffusion signal for a given geometry. Instead, MCDS require an accurate geometrical and physical representation of the diffusion media (called the substrate), a large number of samples, and an acquisition protocol like the classical Pulsed Gradient Spin-Echo (PGSE). In general, an accurate approximation of the diffusion signals, mimicking the ones obtained from the brain's white matter, can be computed if the substrate captures the relevant white matter microstructure features and the simulation parameters are tuned properly (number of samples and their step sizes). Despite this, many simplifications are usually used in order to decrease the computational burden. The most common ones include the use of substrates of small size, the use of a limited number of samples, the use of simple geometries, and the restriction of the 3D diffusion to the 2D case (Lipinski, 1990; Szafer et al., 1995; Fieremans et al., 2010; Dyrby et al., 2013).

Lipinski (1990) presented the first work, to the best of our knowledge, which employs MCDS to study the extracellular diffusion in brain tissue. In this work, 2D histological data was used to draw binary contours to be used as irregular intracellular barriers. From this study on, most studies simplified the representation of the extracellular space as a collection of restricted corridors in the orthogonal plane of the axonal direction and unrestricted parallel to them (Novikov et al., 2011; Dyrby et al., 2013; Sanguinetti and Deriche, 2014; Lin et al., 2016). In addition, the intracellular compartment is usually idealized as a collection of parallel hollow cylinders with constant radii or radii sampled from a distribution estimated from histological data (Fieremans et al., 2008; Hall and Alexander, 2009; Alexander et al., 2010; Raffelt et al., 2012). Recent studies have suggested that such simplification cannot capture the complexity of the axonal structures of white matter, and thus its diffusion characteristics (Nilsson et al., 2013; Ginsburger et al., 2018). For instance, changes in the diffusion signal and parameters derived from the diffusion tensor, such as the fractional anisotropy (FA) and mean diffusivity (MD), were obtained by introducing regular undulations in the intra-axonal compartment (Nilsson et al., 2012).

Because of the aforementioned problems, a number of works proposed experiments where non-trivial structures were used as intracellular substrates, e.g., dispersed axons (Ginsburger et al., 2018), the presence of abutting cylinders (Yeh et al., 2013) and arbitrarily generated meshes (Panagiotaki et al., 2010). However, to this day, such approaches have not been thoroughly adopted by the DW-MRI community because of the high computational burden they demand and the lack of available tools. Because of this, more realistic diffusion simulations remain virtually unexploited.

In this work, we study three important pitfalls encountered in the design of MCDS in the literature used to reduce the computational burden, namely the number of simulated particles and the number of time steps, the intracellular geometrical representation, and the generated extra-axonal space in terms of the substrate's size. Each experiment presented below illustrates a possible bias induced in the computed signal when such simplifications are not properly addressed, which affects its reproducibility. Finally, driven by the results from the previous experiments, we outline a general framework that can be used to generate complex substrates in order to overcome the limitations of previous studies.

2. Theory

The obtained signal from a PGSE DW-MRI measurement, at a time t, is given by (Price, 1997)

$\begin{array}{ll}S(t=TE)=S_0{\displaystyle \int }P(\varphi ,t)e^{-i\varphi }d\varphi ,& (1)\end{array}$

where S₀ denotes the signal obtained in the absence of a diffusion gradient magnetic field, TE is the echo time, P(ϕ, t) is the phase distribution function of the spin ensemble at time t = TE, and ϕ is the accumulated phase shift of the spin.

The amount of attenuation of a single diffusing spin on the measured PGSE signal is proportional to the dephasing due to the effect of the time-dependent magnetic field G(t), and the spin's displacement. For a single spin and a given magnetic gradient vector, the phase shift due to the applied gradient over time can be numerically formulated as in Price (1997)

$\begin{array}{ll}\varphi (t)=a(t)\gamma G(t)·x(t),& (2)\end{array}$

where γ is the gyromagnetic ratio, G(t) is the applied magnetic diffusion gradient at time t, x(t) is the spin's displacement from the starting position, and a(t) is a function that shifts the sign of the gradient vector due to the refocusing Radio Frequency (RF) pulse. In a classic PGSE experiment: a(t) is equal to +1 for all time t before the RF pulse and −1 after. The produced attenuated signal is then the result of the accumulated phase shift of the full assembly of spins at the TE, given by

$\begin{array}{ll}S(t=TE)/S_0={\displaystyle \langle e^{-i{\int }_0^{t=TE}\varphi (t)d{t}^{\prime }}\rangle .}& (3)\end{array}$

2.1. Simulation Fundamentals

Equation (3) can be approximated using a finite number of spin samples N_s over a discrete time lapse, following an approach as in Szafer et al. (1995):

$\begin{array}{ll}S/S_0=\frac{1}{N_s}{\displaystyle \sum ^{N_s}}{e}^{-i{\sum }_t^{Nt}\varphi (t)dt},& (4)\end{array}$

where dt is the step duration, defined as the total diffusion time divided by the number of steps taken (N_t). The value of dt can be fixed as in Hall and Alexander (2009) and Yeh et al. (2013) or normally distributed as in Balls and Frank (2009). In our exploration, we made use of fixed step sizes derived from Einstein's equation: $r=\sqrt{(6Ddt)}$, where D is the diffusion coefficient and r is the expected mean displacement. A fixed step size have been shown by Hall and Alexander (2009) to “reduce the fluctuation in the mean-square displacement of the spins and improve the convergence in the model.” Moreover, Barlett et al. (2013) have shown the fixed step size to be better suited for non-homogeneous systems.

The idea behind a Monte-Carlo Diffusion Simulator is to compute Equation (4) by simulating the particles' Brownian motion and their interaction with respect to a defined substrate. At the beginning of the simulation, the particles are uniformly placed inside the defined substrate's voxel, or substrate's limits. This way, the number of particles in all compartments is proportional to the defined volume fractions. If necessary, the local position of each particle can be tracked to separate the signal contribution of each compartment by, for example, tracking if the particle is inside a given compartment. Over the duration of the simulation, the simulated particles collide and bounce with the substrate's barriers, depending on the barrier properties. Finally, the accumulated phase shift is tracked depending on the spin-echo protocol using Equation (4).

Overall, the formulation above presents an accurate numerical approximation of the diffusion signal based solely in the phase shift distribution. However, is worth noticing that many other effects such as noise levels or the magnetization relaxation should be considered in order to approximate a more realistic DWI-MRI signal.

3. Materials and Methods

All the simulated signals presented below were computed using the sum of the accumulated phase shift approximation showed in Equation (4) implemented in an in-house Monte-Carlo simulator. The simulator employs a similar approach to compute the diffusion signal as the ones presented in Hall and Alexander (2009) and Yeh et al. (2013). The simulator uses a hybrid GPU/Multi-CPU framework, implemented in C++11. It includes routines to optimize the collision detection and the memory consumption based on the complexity analysis of Appendix A; making the software able to handle simulations of 3D meshed substrates with millions of triangles and particles. The simulator was initially validated by verifying that the generated signals from particles within impermeable planes, cylinders, and spheres were equal to those obtained from their corresponding analytical solutions. Moreover, results in more complex domains including the extra-axonal space of brain tissue, were comparable to those obtained from an alternative and independent Finite Element Method approach described in Rafael-Patino et al. (2017). The substrates' data, meshes, and the simulator are available from the corresponding author upon request on the paper's Git-Hub repository: (https://github.com/jonhrafe/Robust-Monte-Carlo-Simulations).

3.1. Confidence Level Estimation

In Monte-Carlo based methods, the number of samples is critical for the confidence level of the estimated results. However, the number of particles has the most significant impact on the computational burden. To highlight the importance of the number of simulated spins, an experiment was performed in order to quantify the variance of the estimated signal as a function of the number of particles sampled in a substrate. To do this, the errors of a set of simulated signals with different numbers of samples were measured. The measured errors were compared against the expected analytical solution in the intra-axonal space and for a gold-standard estimation of the extracellular space. A substrate with 10, 000 parallel cylinders with diameters sampled from a Gamma distribution, Γ(κ, θ), with shape, κ = 4.0, and scale, θ = 4.5 × 10⁻⁷, was used, resulting in a mean diameter μ = 1.8 μm with a standard deviation of σ = 0.9 μm, using a packing algorithm as the one described in Hall and Alexander (2009), which results in a distribution of radii comparable to the ones found in the literature (Zhang et al., 2011; Dyrby et al., 2013; Benjamini et al., 2016).

This substrate was used since the analytical signal of the intra-axonal space can be computed using the volume-weighted sum of the individual signals:

$\begin{array}{ll}{S}_i=\frac{v_{1}S{c}_{i,1}+\cdots +v_{n}S{c}_{i,n}}{\sum v_j},& (5)\end{array}$

where S_i is the ith acquisition, v_j is the volume of the jth cylinder and Sc_{i, 1} is the analytical signal of the cylinder obtained using the Gaussian Phase Distribution (GPD) approximation of the signal in cylinders for a given radius (Van Gelderen et al., 1994). Figure 1 shows the resulting distribution of radii as well as the computed ground-truth intra-axonal signal. For the extracellular signal, as there is no analytical model, the gold-standard was estimated using a very high number of particles: 20 × 10⁶ particles, and time-steps: 2 × 10⁴ steps. These parameters were chosen based on previous results (Rafael-Patino et al., 2017) and by studying the convergence properties for higher numbers of particles and time-steps (Hall and Alexander, 2009). In fact, we verified that the signal converges for even less demanding simulation parameters (i.e., 1 × 10⁶ particles, and 5 × 10³ steps). In order to keep results as accurate as possible, however, we decided to use simulation parameters higher than the minimum required.

FIGURE 1

Figure 1. Gamma distributed radii and corresponding intra-axonal diffusion signal. (Left) The distribution of the sampled diameters, the dotted line marks the sampled distribution mean. (Right) The computed ground-truth along with the simulated signal used for the intra-axonal space representation. A total of four curves are plotted corresponding to each b-value = 1925, 1932, 3093, and 13191 s/mm². The curves corresponding to a b-value = 1925 and 1932 s/mm² are completely overlapped and corresponds to the lowest decay. The signals of each shell are ordered by the normalized Z coefficient of the gradient direction.

The estimated signals were computed varying the number of particles from 1 × 10³ to 1 × 10⁶ particles, and the time-steps from 1 × 10² to 2 × 10⁴ steps. The diffusion coefficient was fixed to D = 0.6 × 10⁻³ mm²/s (corresponding to an ex-vivo diffusivity), and TE = 0.054 s, for both, the simulations and the ground-truth data. The original optimized ActiveAx PGSE protocol (Alexander et al., 2010) was used, which consist of a four shell HARDI acquisition with 90 orientations per shell, each shell with the following parameters, respectively, (i) b = 1, 930 s/mm², G = 140 mT/m, δ = 0.010 s, and Δ = 0.016 s; (ii) b = 1, 930 s/mm², G = 140 mT/m, δ = 0.010 s, and Δ = 0.016 s; (iii) b = 3, 090 s/mm², G = 131 mT/m, δ = 0.007 s, and Δ = 0.045 s; (iv) b = 13, 190 s/mm², G = 140 mT/m, δ = 0.017 s, and Δ = 0.035 s. Figure 1 shows the plot of a diffusion signal obtained with this protocol separated by shell and ordered with respect to the angle with the main fiber axis (Z-axis).

A bootstrapping analysis was performed to evaluate the variance of the error between the estimations with different samples sizes: 1 × 10³, 2 × 10³, 5 × 10³, 1 × 10⁴, 2 × 10⁴, 5 × 10⁴, 1 × 10⁵, 2 × 10⁵, 1 × 10⁶, and 2 × 10⁶ samples; and time-steps: 1 × 10², 5 × 10², 1 × 10³, 5 × 10³, 1 × 10⁴, and 2 × 10⁴. For each combination of the sample sized and time-steps, the signals from 50 repetitions were generated. The error between the ground-truth and each estimated signal was computed using the Relative Mean Absolute Error (RMAE), expressed as a percentage:

$\begin{array}{ll}RMAE(S_{gt},S_c)=\frac{100}{N_g}{\displaystyle \sum _i^N}\frac{|S_{gt}(i)-S_c(i)|}{\left|S_{gt}(i)\right|},& (6)\end{array}$

where S_gt is the ground-truth signal, S_c is the estimated signal and N_g is the number of acquisitions. The result is a total of 50 estimated error points for each sample size.

3.2. Intra-Axonal Space Representation

In our second study, we look into the effect of using curved or angled geometries against straight cylinders as representations of the intra-axonal space. Such effect is of special interest on the computation of axonal diameter indexes when it is assumed that straight cylinders capture the diffusion properties of the intra-axonal compartment.

To understand this effect, an experiment extending the previous work from Nilsson et al. (2012) was performed, where the diffusion properties of undulating axonal substrates is studied. In our experiment, we quantified the difference on the diameter fitting estimation between parallel cylinders of constant radius and undulating cylindrical substrates.

To create curved cylindrical substrates for MCDS, a helical undulation parametrization along z was used

$\begin{array}{ll}U(z)=\left(A_{x}cos(\frac{2\pi z}{L}),A_{y}sin(\frac{2\pi z}{L}),z\right),& (7)\end{array}$

where L is the wavelength and A_x, A_y denote the amplitude in the X and Y axis, respectively (Nilsson et al., 2012). The amplitudes Ax and Ay were set to be equal to obtain helical undulations. Using the formulation above, a set of substrates was created by deforming cylinders with diameters 1, 2, and 3 μm. The wavelength and amplitude of the undulations ranged from 4 to 32 μm and from 0.2 to 2.6 μm, respectively; which covers a range of values of interest in the literature (Haninec, 1986; Bain et al., 2004; Nilsson et al., 2012). The resulting undulating cylindrical shapes were triangulated to use them as mesh substrates suited for MCDS. Figure 2 presents three different substrate examples.

FIGURE 2

Figure 2. Examples of the curved meshes used as intra-axonal substrates in this study, for three different diameters and different undulation parameters.

To compute the diameter estimation error in the intra-axonal signal, a fitting procedure was performed using an exhaustive search approach. The exhaustive search computes the RMAE between the resulting simulated signal of each undulating substrate and the analytical signal of a range of cylinders with different diameters, sampled between 0.4 and 8 μm with a step size of 0.01 μm. The analytical signals were computed using the GPD approximation for the signal in cylinders (Van Gelderen et al., 1994). The fitting procedure returns the range of plausible diameters such that the computed error between them is below a given threshold. For each undulating substrate, the threshold was fixed to a 1% difference from the minimum fitting error, based on the results of the confidence study from section 4.1.

Two different acquisition protocols were used to perform the fitting procedure. First, the original ex-vivo ActiveAx PGSE protocol (Alexander et al., 2010) explained in section 3.1 was used. Second, we used an optimized PGSE protocol for ex-vivo axonal diameter estimation presented in Dyrby et al. (2013). The protocol consists of a three shell acquisition with 90 orientations per shell, and a TE = 0.0359 s. The relevant parameters of each shell are as follows, (i) b = 2081 s/mm², G = 300 mT/m, δ = 0.0056 s, and Δ = 0.0121 s; (ii) b = 3038 s/mm², G = 219 mT/m, δ = 0.007 s, and Δ = 0.0204 s; and (iii) b = 9542 s/mm², G = 300 mT/m, δ = 0.0105 s, and Δ = 0.0169 s. Since the RMAE difference between cylinders of similar diameter depends on the protocol used an analysis of the sensitivity for each protocol was carried out.

Finally, the MC simulation parameters were chosen using a similar analysis as the one presented in section 3.1 (not shown). The confidence estimation was computed ranging the number of particles and time steps on the substrate with higher curvature (higher amplitude and smaller wavelength) and choosing a the parameters that shows almost no variance on the estimations. A total of 5 × 10⁴ particles and 5 × 10⁴ steps were chosen to compute the signal for each individual substrate separately.

3.3. Extra-Axonal Space Representation

In the case of macroscopically homogeneous substrates, e.g., with randomly packed cylinders and in absence of bundle dispersion, it has been shown that extra-axonal spins exhibit an effective diffusivity that can be described by an axi-symmetric tensor, if the volume size of the sample is high enough (Hrabe et al., 2004). Models to estimate white matter microstructure from DW-MRI therefore assume that the extra-axonal radial contribution does not change for any direction aligned to the bundle's perpendicular plane (Assaf et al., 2008; Alexander et al., 2010; Zhang et al., 2011, 2012; Panagiotaki et al., 2012; Daducci et al., 2015a; Benjamini et al., 2016).

Such an assumption seems to fit the validations. However, the importance of the design of the extra-axonal space has been underestimated in MCDS by assuming that substrates with any hindered configuration would match the model. To show the importance of the sample size, in terms of the number of cylinders used to construct a substrate, an analysis of the extra-axonal radial contribution in simulated signals was performed.

The radial extra-axonal DW-MRI signal was simulated for a selection of voxels with different numbers of cylinders, and a fixed distribution of diameters and intra-axonal volume fractions. To do so, N diameters (N = 100, 1, 000, 10, 000, 50, 000 and 100, 000) were sampled from a gamma distribution with parameters Γ(4.0, 4.5 × 10⁻⁷), as in our first study. The corresponding cylinders were randomly positioned in substrates with voxel size adapted such that the intra-axonal volume fraction was 60% and ensuring periodicity at the voxel boundaries as is described in Hall and Alexander (2009). The extra-axonal signal was simulated with the following settings: 1 × 10⁶ particles in the extra-axonal space with diffusivity of 0.6 × 10⁻⁴ mm²/s, TE = 0.075 s, and 1 × 10³ steps. This parameters where chooseng from the previous results showed in section 3.1. The diffusion protocol was set to highlight the radial contribution of the diffusion signal in different diffusion time regimes as follows: G = 300 mT/m, δ = 0.010 s and Δ from 0.015 to 0.060 s, acquired in 180 directions evenly distributed over the xy-plane. The anisotropy of the simulated noiseless signal was quantified by computing the standard deviation of the signal divided by its mean, giving an estimate of how much the signal deviates from a perfectly radially isotropic signal.

3.4. Framework for Complex Substrates Generation

Driven by the results from the previous experiments, and based on a previously published algorithm to generate tractography phantoms (Close et al., 2009) the following section outlines a general framework in order to generate complex substrates. We show that such framework overcomes some of the simplifications presented in the previous sections. To illustrate such capabilities, a crossing of axons bundles was generated as a study case. A qualitative evaluation was performed over the representation of crossing fibers in terms of the resulting intra-axonal volume fraction and diffusion properties in different resolutions.

The framework is a tailored extension of the work presented in Close et al. (2009). The original framework is based in the optimization of a objective function that penalizes the overlap, curvature and length of a set of initial splines called as strands. Each strand has a constant radius used to ensure no overlapping. The optimization cost-function has the following form:

$\begin{array}{ll}E(\cup S)={\displaystyle \sum _i^{\#S}}{w}_o{J}_o(S_i)+w_c{J}_c(S_i)+w_l{J}_l(S_i),& (8)\end{array}$

where the set ∪S of size #S, represents the set of all initialized strands. S_i represents the strand i for i = 1, ⋯ , #S. The functions J_o(·), J_c(·), J_l(·) are the overlap, curvature, and length penalization functions; and the coefficients w_o, w_c, and w_l are their, respectively, weights. Each strand S_i is parametrized using a discrete set of control points that define the backbone of the strand i and constant given radius; the transversal area associated with this radius is later subdivided to form sub-strands. Finally, the DW-MRI signal is then simulated by assigning a symmetric tensor along each sub-strand trajectory i.e., a simplistic local model of the micro-environment (Daducci et al., 2015b). The reader is referred to Close et al. (2009) for more details.

In our study, the aforementioned framework was modified and used to map a gamma distributed set of diameters inside the resulting strands' trajectories. The 3D-overlapping algorithm between strands implemented in the cost-function J_o, was also modified to make it more suitable for creating 3D meshes. This was done by computing the analytical intersection between two strands' control points, using the cylinder to cylinder collision detection described in Verth and Bishop (2008).

The result is a gamma distributed crossing configuration of deformed cylinders. The main advantage of this configuration is that the bundles inside a common area do not overlap or intersect, but interdigitate, which means that the volume is preserved in the crossing region. In addition, the curvature and length penalizations promotes a higher packing density. Finally, the proposed framework computes the DW-MRI signal by a Monte-Carlo simulation using a mesh substrate created from the configuration obtained above, instead of assigning a symmetric tensor along the sub-strands. Figure 3 shows the crossing configuration before and after the optimization procedure.

FIGURE 3

Figure 3. Optimization procedure of initial trajectories. (Left) initial trajectories parametrized as a set of control points with constant radius. (Right) the resulting trajectories after the optimization procedure which ensures that there is no overlapping between the resulting strands.

In the presented study case, the diameters from a gamma distribution with parameters Γ(1.2, 1.5 × 10⁻⁶) were sampled, resulting in a mean diameter of μ = 1.8 μm and standard deviation σ = 1.6 μm, which are in the range of anatomical interest (Alexander et al., 2010). The resulting values were truncated to avoid strands with diameters smaller than 0.2 μm. The dimensions of the resulting enclosing volume were 1,200 × 240 × 480 μm; the resulting 3D geometrical crossing is shown in Figure 4. The 3D mesh model consists of 1,698,328 triangular faces after a post-processing of decimation and smoothing to reduce the triangle density. The total length end-to-end of the most extended strand is 1.58 mm. The resulting diameter distribution of the overall structure is displayed in the bottom panel of Figure 4.

FIGURE 4

Figure 4. (Top panel) Shows a visualization of the resulting fiber crossing substrate after the strand refinement and the smoothing and decimation of the triangular faces. (Left-bottom panel) Shows the resulting sub-strand configuration of one of the crossings bundles. (Right-bottom panel) Shows the overall diameter distribution of the displayed bundle on the (left). A rendered video of the full crossing is included as Supplementary Material.

To compute the simulated MRI signal, the total volume was divided in three voxel resolutions: 80 × 16 × 32, 40 × 8 × 16, and 20 × 4 × 8 voxels. A total of 105 × 10⁶ particles, and 5, 000 steps were used to compute the signal for the three resolutions. The original ActiveAx protocol (Alexander et al., 2010) from the first study was used with a diffusivity coefficient equal to 0.6 × 10⁻³ mm²/s and a total diffusion time of 0.053 s.

To show qualitative results on the generated signals, the Diffusion Tensor (DT) estimation and the corresponding FA were computed using Dipy (Garyfallidis et al., 2014), as well as the ICVF maps for each of the three resolutions. Only the shell with b = 3, 080 s/mm² was used to compute the DT in each voxel. Given the lack of an analytic representation of the substrate, the ICVF was approximated by tracking the local position of the uniformly random located particles and labeling them as inside or outside the meshed substrate.

Finally, an evaluation of the axon diameter estimation within the crossing area was performed for the three different voxel resolutions. The axon diameter estimation was performed using the same exhaustive search method described in section 3.2. Only one single bundle orientation was used to compute the analytical GPD approximation; which was selected from the DT estimation at each voxel. The fitting procedure was performed using solely the intra-axonal signal and in voxels with FA greater than 0.25, in order to separate the effect of the extra-axonal space regarding the diameter mis-estimation.

4. Results

4.1. Confidence Level Estimation

The overall results of the bootstrapping analysis are summarized in Figures 5, 6 for the intra- and extra-axonal space, and for both, the number of samples, and the number of time-steps. Figure 5 shows the mean error of the 50 samples for each one of the possible combination of the selected parameters, color-coded in a heat-map. In Figure 5, we show the error of each repetition by (i) fixing the number of steps to the maximum value (2 × 10⁴) and varying the number of particles (left column), and (ii) fixing the number of particles to the maximum value (2 × 10⁶) and varying the number of steps (right column). Each data point represents one repetition of a given sample size. A total of 50 points are plotted in each row, and the mean error for each sample is highlighted with a red asterisk. The total simulation time for each repetition ranged from few seconds for the simulation with a total of 1 × 10³ particles to 918 s for the simulation with 1 × 10⁶ particles. Each simulation was performed in a single node of Fidis EPFL's cluster with 14 cores, 2.6 GHz, and 528 MB of RAM.

FIGURE 5

Figure 5. RMAE for each repetition and sampled size for (Left) the number of samples and (Right) number of time-steps. The two panels on the top row correspond to the intra-axonal results, and the bottom row to the extra-axonal. The X-axis shows each sample size, and the Y-Axis shows the RMAE of all the repetition in same color. The mean RMAE of all the repetitions is depicted with a red marker.

FIGURE 6

Figure 6. Heat map of the mean RMAE for all the combinations between the number of steps and the number of samples. Each cell corresponds to the mean RMAE of the 50 repetitions.

For the study regarding the number of particles in the intra-axonal space, the mean RMAE between the analytical ground truth and the set of repetitions with the biggest sample size of 2 × 10⁶ particles was of 0.47%. For the extra-axonal space, the mean RMAE between the computed gold-standard with 20 × 10⁶ and the set with 1 × 10⁶ particles was equal to 0.71%. For the analysis varying the number of time-steps, the minimum mean RMAE achieved was of 0.44% for the intra-axonal space and 0.38% for the extra-axonal. The difference between the mean RMAE between 5 × 10³ and 2 × 10⁴ was less than 0.2% for both the intra- and extra-axonal space.

4.2. Intra-Axonal Space Representation

The range of diameters, computed from our fitting procedure on both protocols, are displayed in Figure 7. Each cell is colored according to its minimum RMAE. An amplitude (amp) of 0 μm corresponds to a straight cylinder which presented the minimum fitting error achievable for each diameter. Values with the highest amplitude and lowest wavelength (wl) corresponds to the axons with the highest undulation (amp = 2.6 μm, wl = 4 μm); on the other hand, values with the lowest amplitude and highest wavelength (amp = 0.2 μm, wl = 32 μm) corresponds to almost straight axons.

FIGURE 7

Figure 7. Tables of the fitting results. Left column shows the fitted intervals of the original ex-vivo ActiveAx protocol (Alexander et al., 2010) and the right column of the optimized ex-vivo protocol from Dyrby et al. (2013), for the three simulated diameters. The min and max diameters ( μm) of the fitted range are listed between the square brackets for each simulated amplitude and wavelength. The color of each cell is encoded with respect the minimum RMAE in the fitted range accordingly to the color-bar on the right.

The protocols' sensitivity analysis is shown in Figure 8 which presents a visualization of the RMAE between the analytical signal of straight cylinders with different diameters. Regions with homogeneous values are difficult to differentiate between each other, e.g., the region with diameters between 0 μm and 2.0 μm. The colored line shown in both plots marks the 1% level curve. In this plot, the protocols' contrast in function of the cylinder's diameter can be visualized, which correlates with the intervals showed in Figure 7.

FIGURE 8

Figure 8. The in-between RMAE of the analytical signal of a cylinder, obtained using the GPD approximation, for the range of diameters used in this study. The original ex-vivo ActiveAx acquisition protocol (Alexander et al., 2010) is displayed on the left panel, and the optimized ex-vivo acquisition protocol from Dyrby et al. (2013) on the right panel. Values of the diagonal correspond to the RMAE of two straight cylinders of the same diameter and therefore equals to 0. The colored line shown in both plots marks the 1% difference level curve.

4.3. Extra-Axonal Space Representation

Three different substrates (with 100, 1, 000 and 10, 000 cylinders, corresponding to voxel sizes of 23 × 23 μm, 71 × 71 μm, and 230 × 230 μm, respectively) and their corresponding radial DW-MRI signals, are shown in Figure 9. The shown voxel sizes were chosen to highlight the radial anisotropy in three representative sizes. The substrate with 10, 000 cylinders, i.e., with the biggest voxel size, had the most isotropic radial DW-MRI signal. On the other hand, the most anisotropic signal was observed for the substrate with the smallest number of cylinders. Figure 10 shows the mean and standard deviation of the radial extra-axonal signal as a function of the voxel size. The same experiment (not shown) was conducted using cylinders with higher diameter. Results indicated that the number of cylinders was the limiting factor. Indeed, the mean of the radial extra-axonal signal also converged for 10,000 cylinders, but this time a voxel size of approximately 400 × 400 μm was required to generate isotropic profiles instead of the 230 × 230 μm limit observed for a distribution with smaller cylinders.

FIGURE 9

Figure 9. Results for 3 substrates, with 100, 1,000, and 10,000 cylinders, respectively. First row: sampled diameter distributions for each voxel-size, getting closer to the desired distribution law as the voxel-size increases. Second row: cylinder positions in each substrate. White scale bar corresponds to 10 μm. Third row: radial DW-MRI signal simulated from the respective substrates. Each colored line corresponds to one different Δ duration. Dotted lines correspond to the mean radial signal for each diffusion time.

FIGURE 10

Figure 10. Mean and standard deviation of the radial DW-MRI signal as a function of substrate size. The signal is shown for each of the different Δ.

4.4. Framework for Complex Substrates Generation

The resulting crossing with two fiber populations is outlined in Figure 4. The total optimization time to create the substrate was around 42 h, where most of the optimization time (about 35 h) was needed in the second optimization iteration, after the subdivision on gamma distributed radii, that ensure that no small overlaps were introduced due to the subdivision and abrupt angular changes. The optimization was performed using a single core 2.8 GHz CPU. On the other hand, the total simulation time for the full geometry with 105 × 10⁶ particles was less than 24 h using a total of 8 nodes with 28 cores on fidis EPFL's cluster with 6GB of RAM per node (48GB in total).

The resulting diffusion tensor and FA maps are shown in Figure 11 for the three different resolutions. Local diffusivity changes, as well as signal alterations related to the curvature of the individual axons, can be observed.

FIGURE 11

Figure 11. From the leftmost to the right: diffusion tensor map, the resulting fractional anisotropy and the two highlighted ROIs in each map, respectively. Each image corresponds to the same volume slice in the XZ-plane. The ROI's highlights one area where different compartments result from the optimization procedure.

Figure 12 shows the intra-axonal volume fraction in all resolutions. In the highest resolution, small water compartments can be seen in the crossing sections; this is an effect of the optimization procedure which ensures no overlapping fibers. In the lowest resolution, such compartments are no longer visible, but they are reflected in the decrease of the intra-axonal volume fractions.

FIGURE 12

Figure 12. The ICVF maps of one volume slice in the XZ-plane in three different resolutions. The highest achieved ICVF value for each resolution were: 0.8013, 0.5792, 0.4825, from top to bottom, respectively. The two green areas highlighted in the two lowest resolutions were used to evaluate the axon diameter estimation.

Figure 13 shows a visualization of one plane of the axon diameter estimation maps of the volumetric region highlighted in Figure 12, and the obtained diameter distribution for the three resolutions. The higher resolution (80 × 16 × 32) estimation includes a total of 2,848 voxels, while the lowest resolution contains a total of 112 voxels. Figure 4 bottom-right panel shows the resulting sampled diameters inside the crossing configuration, which is noticeably skewed to smaller diameters; this is an effect of the packing algorithm inside individual circular strands which under-represent the tail of the distribution because of the difficulty of packing strands with big diameters. This effect will irremediably affect the effective apparent radius $r_{eff}\approx {(\frac{<r^6>}{<r^2>})}^{1/4}$ (Burcaw et al., 2015) given by the intra-axonal contribution of the signal. The resulting effective diameter of the conjoint assemble of strands was 2*r_eff = 3.48 μm, which in average agrees with the estimated mean diameters shown in Figure 13 on each resolution.

FIGURE 13

Figure 13. Axon diameter estimation maps (Left column) of the regions highlighted in Figure 12, and diameter histograms (Right column) estimated on the full volume enclosed by the highlighted regions. Top row shows the axon diameter map and the diameter estimation histogram for the 80 × 16 × 32 nominal resolution; middle row shows the same maps for the 40 × 8 × 16 nominal resolution, and the bottom row shows the same maps for the 20 × 4 × 8 nominal resolution. The dotted line indicates the histograms' mean diameter within the regions, to be compared with the effective apparent diameter (2*r_eff) of 3.48 um.

5. Discussion

In the past two decades, the research community has used MCDS to generate and validate MR diffusion data and microstructure models (Lipinski, 1990; Hall and Alexander, 2009; Fieremans et al., 2010; Panagiotaki et al., 2010; Nilsson et al., 2012; Barlett et al., 2013; Baxter and Frank, 2013; Plante and Cucinotta, 2013). However, questions have been raised on the accuracy of the simplified geometries used to create the diffusion substrates (Balls and Frank, 2009; Panagiotaki et al., 2010; Nilsson et al., 2012, 2013), emphasizing the need of highly-validated and reproducible simulations. Such oversimplifications have been proven not to capture the complexity of the axonal structures of white matter, and thus its diffusion characteristics (Nilsson et al., 2013). Moreover, it can be argued that the use of such elementary geometries-used as backbone in the microstructure models-as a ground-truth, not only introduce a systematic bias that inherently supports the evaluated method, but also misapplies the very purpose of using Monte-Carlo simulations. In this work, we outlined pitfalls encountered in the design of such simulations. Our experiments showed how the design of each substrate compartment is likely to introduce an estimation bias if it is not addressed appropriately.

Our first study specifically shows the effect of an inappropriate selection of parameters on the reproducibility of a estimated signal, which could also skew an analysis toward inaccurate results. Differently to previous studies (Hall and Alexander, 2009), we compute the extra-axonal ground-truth from high-quality simulations, avoiding the use of tortuosity models that could introduce a bias because of their oversimplifications. The error in the estimation presented in Figure 5 illustrates the great amount of possible estimation variability for a relatively simple substrate. We found that the signal on each compartment showed a high variability for simulations with less than 5 × 10⁵ particles and 1 × 10⁴ steps. We can extrapolate from this that any estimation from more complicated substrates, such as the ones with undulation or crossings, or even higher diffusivity, will likely entail even higher variability. In order to avoid such uncertainty on the estimations for more complicated substrates a similar analysis as the one presented should be procured.

In our second study, we explored the effect of breaking the assumption of straight cylinders as the intra-axonal representation in function of the apparent diameter estimation. The helical representation used in this study, while reported to appear in the nervous system, maybe not be an accurate representation of the axonal angular variations along the longitudinal direction in the brain white matter, specially in the micro-scale. However, it gives us a convenient starting point to study the effect of angular variations in the intra-axonal compartment over the diffusion signal, a theoretical analysis on this type of structures can be found as well in Brabec et al. (2019). From this study, we found a considerable mis-estimation in the presence of undulation for both protocols and in the three studied diameters. The relative fitting error for the smaller diameter (1 μm) was the higher among the three cases (more than 300% for some cases). Previously, Nilsson et al. (2017) proposed a formulation to compute the minimal diameter of a parallel cylinder able to produce a signal attenuation larger than that from a cylinder with a diameter of zero, using standard single-shell PGSE sequences. According to this formalism, the minimum differentiable diameter is $d_{min}={(768\sigma D/7{\gamma }^2\delta |G{|}^2)}^{1/4}$, where σ is the significance level, defined as the minimum tolerated percentage of signal change. For a fixed value σ = 1% change, the resolution limit predicted for both protocols used in this study were d_min = 2.29 μm for |G| = 140mT/m, and d_min = 1.76 μm for |G| = 300mT/m. However, such estimates are based on a number of assumptions which does not hold in our experimental conditions. For example, the formulation is valid for parallel and straight cylinders and for acquisition protocols with a single shell with parameters δ = Δ. As in this experiment we are studying non-parallel and curved cylinders with multi-shell protocols with Δ>>δ, we performed a numerical sensitivity analysis to obtain more accurate resolution limits. From the resulting plots showed in Figure 8, it can be seen that the signal originated from cylinders with diameters below 2.5 μm for the first protocol, and 2.0 μm for the second, are virtually indistinguishable. On the other hand, diameters above 3.0 μm have more significant RMAE, which make them easier to differentiate. We also observed that the range of diameters from our fitting method did not follow a simple trend between protocols; that is to say, increments on the undulation parameters, which effect can be summarized in terms of the tortuosity factor (Nilsson et al., 2012) $\lambda =\sqrt{{(2\pi A/L)}^2+1}$, does not follow a simple relationship between protocols (horizontal axis of the results in Figure 7). This is likely to be an effect of the parameters of the acquisition protocol (δ, Δ, and the TE), which vary between shells and thus changing the effective diffusion time. From a comparison of both protocols, we corroborated that the optimized protocol showed better results in terms of the fitted diameter and range of similar diameters. However, there was still a considerable mis-estimation, especially for the undulation of 1 μm diameters. We consider this experiment to be of great interest for any future protocol optimization or diameter estimation framework, since it illustrates how sensible the estimation of the axon's diameter based on the cylindrical model are, even for regular and smooth angular deviations.

Our third experiment showed that a sufficiently rich sampling is required for the simulated signal to converge. Indeed, small substrates have a limited number of cylinders, limiting the variability of hindered micro-environments sampled by the spins during the M.C. simulation—yielding anisotropic patterns in the radial DW-MRI signal. The results also showed a bias in the mean amplitude, with small voxels having lower signal than bigger voxels. Our results suggest that, for a given diameter distribution, substrates with an area smaller than 200 × 200 μm will present biased extra-axonal signals. Such results are in accordance with previously results (Hall et al., 2017) in terms of the voxels' size. However, this lower bound probably depends on the distribution of diameters and cylinder packing on one side, as well as the typical diffusion length of the spins, given by their diffusivity and the diffusion time of the experiment.

Finally, as part of our effort to create more realistic substrates, we outlined a framework to tackle the challenging problem of creating non-overlapping crossing configurations that preserves the volume fractions between the non-crossing and crossing area, while enforcing a high packing density. Configurations which mimic better real tissue (Shacklock, 2007), are important since they provide a more challenging environment to test and validate microstructure models and even tractography methods, in contrast with naive crossings which have been proven to be indistinguishable from a simple superposition of individual fascicles (Rensonnet et al., 2018). From the diffusion tensor and FA maps shown in Figure 11 we can observe the presence of multiple compartments as an effect of the volume preservation condition. Also, Figure 12 shows how the intra-axonal volume fraction changes as the resolution decreases. Such information can be used to study the microstructure information in the presence of several diffusion compartments and volume fractions in different resolutions without using an explicit interpolation. This decrease of the ICVF is an effect of the presence of dispersion and deformation of the fiber bundles. However, even in the lowest resolution, the intra-axonal volume fraction achieved was over 48%, which is considerably higher than the icvf (of 20%) of a previously presented framework for generating realistic numerical phantoms for crossing fascicles (Ginsburger et al., 2018). By optimizing the penalization term of the strands' curvature in our framework, we expect to be able to achieve even higher packing densities—closer to the expected ones from the brain's white matter tissue. On the other hand, the diameter estimations computed over the merging area of the two fiber populations, showed a overestimation in accordance with the results of section 3.2. Such mis-estimation can be explained by angular perturbations in the fiber trajectory in both the micro- and meso- scale of the simulated fibers. In previous studies, the axon diameters were overestimated by factors 3-5 in clinical scanners (Alexander et al., 2010; Zhang et al., 2011). This bias was attributed to the insensitivity of the measurement schemes to small axons (Dyrby et al., 2013), the noise, or the commonly neglected time-dependence of diffusion in the extra-axonal space (De Santis et al., 2016). The diameters reported in this study were estimated by using only the intra-axonal signals, thus the overestimation can be explained only by the dMRI signal insensitivity to the smaller axons and by the signal changes due to axon undulations and microscopic dispersion. This renders our estimations as a best case scenario.

5.1. Considerations and Future Work

The generalisability of the results presented above is subject to certain limitations. For instance, in-vivo diffusion and protocol settings, the use of non-regular deformations in the intra-axonal substrates, and the joint study of the intra- and extra- axonal space, may affect the results toward higher variability or mis-estimations of the axon diameters. In addition, a number of structural features present in white matter tissue -such as the axonal myelin sheath, Ranvier nodes, or diameter changes along the axons trajectory-are missing. Because of this, the results presented above should be taken as a type of lower bound in terms of the minimum parameters needed (for the number of samples and time-steps) and possible mis-estimations (in terms of our axon diameters estimates).

Notwithstanding these limitations, we consider that the aforementioned framework, complemented with the optimized simulator developed, are able to overcome the simulations pitfalls presented in this work. In addition, the parameter selection analysis presented in this work provides a way to ensure the reproducibly of the Monte-Carlo simulations. A thorough study of the properties of more complex substrates generated with the proposed framework is beyond the scope of this study. Future research should therefore concentrate on the generation and study of such configurations, which may help the DW-MRI research community to generate more reliable ground-truth data.

6. Conclusions

The main contribution of this work can be summarized in three main aspects. First, this paper outlines and investigates a set of pitfalls encountered on the parameter selection and substrates' design for Monte-Carlo simulations. Our results over the effect of the number of particles and time-steps, as well as our quantification over the effect of the substrate's size on the extra-axonal space can be immediately taken to evaluate the design of future experiments. In overall, we found that for experiments with parameters in the range used in this study—which are in the range of interest in the literature—simulations with less than 5 × 10⁵ particles and 1 × 10⁴ steps carried a significant variance between the computed signals for both, the intra- and extra-axonal compartments. In addition, we found that simulations substrates with less than 10,000 sampled cylinders induced an important bias on the directional symmetry of the diffusion signal in directions transversal to the main fiber direction. Such parameters are almost one order of magnitude bigger than previously used on the literature, which inherently affects the reproducibility of such results (Hall and Alexander, 2009; Alexander et al., 2010; Rensonnet et al., 2018, 2019). Second, our evaluation over the effect of introducing angular perturbations in the intra-axonal space representation—by means of the estimated axon diameter based on the cylindrical model—showed a considerable deviation from the expected results. This results are somehow in agreement with previous findings and contributes additional evidence that suggests that performing whole brain axon diameter estimation is still far from being straightforward using simplified models, such as the straight cylindrical diffusion model. Finally, this paper presents a framework able to generate complex fiber configurations with desired microstructure information based on a previous algorithm used to create tractography phantoms. We showed the framework's capabilities to generate complex fibers configurations which, along with the simulator developed in this work, are able to generate more challenging and composite Monte-Carlo simulations.

We consider that the results presented in this work, along with the reported procedure to evaluate the estimations' variability, the substrate generation framework, and the simulator developed, pave the way toward more realistic and reproducible Monte-Carlo simulations for Diffusion-Weighted MRI.

Data Availability Statement

The datasets generated for this study are available in: https://github.com/jonhrafe/Robust-Monte-Carlo-Simulations.

Author Contributions

JR-P: original idea conception and experimental design, overall writing, plots and figures generation, and lead programmer of the simulator presented in this work. DR: section 3.3 experiment setup, generation of section 4.3 results and writing based on previously published work, and active contributions on the experimental setup and selected parameters. AR-M: contributed on sections 3.2 and 4.2 experiments design and writing, and overall review of the work and editing. EC-R: contributed on sections 4.1 and 3.1 experimental setup and parameters choice, overall review of the work and editing. GG: generation of the maps presented in section 4.4 coding and experimental setup, Figures 10, 11 generation based on previous work, and overall review of the work and editing. J-PT: supervised all the work done, contributed to the general writing, and editing of the work.

Funding

This project has received funding from the Swiss National Science Foundation under grant number 205320_175974. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan Xp GPU used for this research. This work was supported by EPFL through the use of the facilities of its Scientific IT and Application Support Center. This project has received funding from the European Union's Horizon 2020 Framework Programme for Research and Innovation under grant agreement 665667 (call 2015) and from the Natural Sciences and Engineering Research Council of Canada (NSERC).

Conflict of Interest

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.

Acknowledgments

We would like to acknowledge Dr. Gary Zhang for his advice that highly contributed to the development of this work. We also thank Thomas Yu for his generous help proofreading this manuscript.

Supplementary Material

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

Supplementary Video 1. 3d rendering of the resulting fiber-crossing substrate after refinement, and the smoothing and decimation of the triangular faces.

References

Alexander, D. C., Hubbard, P. L., Hall, M. G., Moore, E. A., Ptito, M., Parker, G. J., et al. (2010). Orientationally invariant indices of axon diameter and density from diffusion MRI. NeuroImage 52:1374. doi: 10.1016/j.neuroimage.2010.05.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Assaf, Y., Blumenfeld-Katzir, T., Yovel, Y., and Basser, P. J. (2008). Axcaliber: a method for measuring axon diameter distribution from diffusion MRI. Magn. Reson. Med. 59:15334406. doi: 10.1002/mrm.21577

PubMed Abstract | CrossRef Full Text | Google Scholar

Bain, A. C., Shreiber, D. I., and Meaney, D. F. (2004). Modeling of microstructural kinematics during simple elongation of central nervous system tissue. J. Biomech. Eng. 125, 798–804. doi: 10.1115/1.1632627

PubMed Abstract | CrossRef Full Text | Google Scholar

Balls, G. T., and Frank, L. R. (2009). A simulation environment for diffusion weighted MR experiments in complex media. Magn. Reson. Med. 62:771. doi: 10.1002/mrm.22033

PubMed Abstract | CrossRef Full Text | Google Scholar

Bammer, R. (2003). Basic principles of diffusion-weighted imaging. Eur. J. Radiol. 45:169. doi: 10.1016/S0720-048X(02)00303-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Barlett, V. R., Hoyuelos, M., and Martin, H. O. (2013). Monte carlo simulation with fixed steplength for diffusion processes in nonhomogeneous media. J. Comput. Phys. 239, 51–56. doi: 10.1016/j.jcp.2012.12.029

CrossRef Full Text | Google Scholar

Baxter, G. T., and Frank, L. R. (2013). A computational model for diffusion weighted imaging of myelinated white matter. NeuroImage 75, 204–212. doi: 10.1016/j.neuroimage.2013.02.076

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, D., Komlosh, M. E., Holtzclaw, L. A., Nevo, U., and Basser, P. J. (2016). White matter microstructure from nonparametric axon diameter distribution mapping. Neuroimage 135, 333–344. doi: 10.1016/j.neuroimage.2016.04.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Brabec, J., Lasic, S., and Nilsson, M. (2019). Time-dependent diffusion in undulating structures: impact on axon diameter estimation. NMR Biomed. 33:e4187. doi: 10.1002/nbm.4187

CrossRef Full Text | Google Scholar

Burcaw, L. M., Fieremans, E., and Novikov, D. S. (2015). Mesoscopic structure of neuronal tracts from time-dependent diffusion. NeuroImage 114:15334406. doi: 10.1016/j.neuroimage.2015.03.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Close, T. G., Tournier, J. D., Calamante, F., Johnston, L. A., Mareels, I., and Connelly, A. (2009). A software tool to generate simulated white matter structures for the assessment of fibre-tracking algorithms. NeuroImage 47:1288. doi: 10.1016/j.neuroimage.2009.03.077

PubMed Abstract | CrossRef Full Text | Google Scholar

Daducci, A., Canales-Rodríguez, E. J., Zhang, H., Dyrby, T. B., Alexander, D. C., and Thiran, J. P. (2015a). Accelerated microstructure imaging via convex optimization (amico) from diffusion MRI data. NeuroImage 105:32. doi: 10.1016/j.neuroimage.2014.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Daducci, A., Palù, A. D., Lemkaddem, A., and Thiran, J. P. (2015b). Commit: convex optimization modeling for microstructure informed tractography. IEEE Trans. Med. Imaging 34:246. doi: 10.1109/TMI.2014.2352414

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Dyrby, T. B., Sagaard, L. V., Hall, M. G., Ptito, M., and Alexander, D. C. (2013). Contrast and stability of the axon diameter index from microstructure imaging with diffusion MRI. Magn. Reson. Med. 70:711. doi: 10.1002/mrm.24501

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferizi, U., Schneider, T., Witzel, T., Wald, L. L., Zhang, H., Wheeler-Kingshott, C. A., et al. (2015). White matter compartment models for in vivo diffusion MRI at 300 mt/m. NeuroImage 118, 468–483. doi: 10.1016/j.neuroimage.2015.06.027

CrossRef Full Text | Google Scholar

Fieremans, E., Deene, Y. D., Delputte, S., Özdemir, M. S., D'Asseler, Y., Vlassenbroeck, J., et al. (2008). Simulation and experimental verification of the diffusion in an anisotropic fiber phantom. J. Magn. Reson. 190:189. doi: 10.1016/j.jmr.2007.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Fieremans, E., Novikov, D. S., Jensen, J. H., and Helpern, J. A. (2010). Monte carlo study of a two-compartment exchange model of diffusion. NMR Biomed. 23:711. doi: 10.1002/nbm.1577

PubMed Abstract | CrossRef Full Text | Google Scholar

Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., van der Walt, S., Descoteaux, M., et al. (2014). Dipy, a library for the analysis of diffusion MRI data. Front. Neuroinform. 8:8. doi: 10.3389/fninf.2014.00008

PubMed Abstract | CrossRef Full Text | Google Scholar

Ginsburger, K., Poupon, F., Beaujoin, J., Estournet, D., Matuschke, F., Mangin, J.-F., et al. (2018). Improving the realism of white matter numerical phantoms: a step toward a better understanding of the influence of structural disorders in diffusion MRI. Front. Phys. 6:12. doi: 10.3389/fphy.2018.00012

CrossRef Full Text | Google Scholar

Hall, M. G., and Alexander, D. C. (2009). Convergence and parameter choice for monte-carlo simulations of diffusion MRI. IEEE Trans. Med. Imaging 28, 1354–1364. doi: 10.1109/TMI.2009.2015756

PubMed Abstract | CrossRef Full Text | Google Scholar

Hall, M. G., Nedjati-Gilani, G., and Alexander, D. C. (2017). Realistic voxel sizes and reduced signal variation in monte-carlo simulation for diffusion mr data synthesis. arXiv: 1701.03634.

Google Scholar

Haninec, P. (1986). Undulating course of nerve fibres and bands of Fontana in peripheral nerves of the rat. Anat. Embryol. 174, 407–411. doi: 10.1007/BF00698791

PubMed Abstract | CrossRef Full Text | Google Scholar

Hrabe, J., Hrabetová, S., and Segeth, K. (2004). A model of effective diffusion and tortuosity in the extracellular space of the brain. Biophys. J. 87, 1606–1617. doi: 10.1529/biophysj.103.039495

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, M., He, H., Schifitto, G., and Zhong, J. (2016). Simulation of changes in diffusion related to different pathologies at cellular level after traumatic brain injury. Magn. Reson. Med. 76:290. doi: 10.1002/mrm.25816

PubMed Abstract | CrossRef Full Text | Google Scholar

Lipinski, H. G. (1990). Monte carlo simulation of extracellular diffusion in brain tissues. Phys. Med. Biol. 35:441. doi: 10.1088/0031-9155/35/3/012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lovas, G., Szilágyi, N., Majtényi, K., Palkovits, M., and Komoly, S. (2000). Axonal changes in chronic demyelinated cervical spinal cord plaques. Brain 123:308. doi: 10.1093/brain/123.2.308

PubMed Abstract | CrossRef Full Text | Google Scholar

Neuman, C. H. (1974). Spin echo of spins diffusing in a bounded medium. J. Chem. Phys. 60:4508. doi: 10.1063/1.1680931

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Nilsson, M., Latt, J., Staahlberg, F., van Westen, D., and Hagslatt, H. (2012). The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study. NMR Biomed. 25:795. doi: 10.1002/nbm.1795

PubMed Abstract | CrossRef Full Text | Google Scholar

Nilsson, M., Westen, D. V., Staahlberg, F., Sundgren, P. C., and Latt, J. (2013). The role of tissue microstructure and water exchange in biophysical modelling of diffusion in white matter. Magn. Reson. Mater. Phys. Biol. Med. 26:345. doi: 10.1007/s10334-013-0371-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Novikov, D. S., Fieremans, E., Jensen, J. H., and Helpern, J. A. (2011). Random walks with barriers. Nat. Phys. 7:1004.2701. doi: 10.1038/nphys1936

CrossRef Full Text | Google Scholar

Panagiotaki, E., Hall, M. G., Zhang, H., Siow, B., Lythgoe, M. F., and Alexander, D. C. (2010). “High-fidelity meshes from tissue samples for diffusion mri simulations,” in Medical Image Computing and Computer-Assisted Intervention - MICCAI 2010 (Beijing), 404.

PubMed Abstract | Google Scholar

Panagiotaki, E., Schneider, T., Siow, B., Hall, M. G., Lythgoe, M. F., and Alexander, D. C. (2012). Compartment models of the diffusion mr signal in brain white matter: a taxonomy and comparison. NeuroImage 59:2241. doi: 10.1016/j.neuroimage.2011.09.081

PubMed Abstract | CrossRef Full Text | Google Scholar

Plante, I., and Cucinotta, F. A. (2013). “Monte-carlo simulation of particle diffusion in various geometries and application to chemistry and biology,” in Theory and Applications of Monte Carlo Simulations, ed V. Chan, 193. doi: 10.5772/53203

CrossRef Full Text | Google Scholar

Price, W. S. (1997). Pulsed-field gradient nuclear magnetic resonance as a tool for studying translational diffusion: part 1. Basic theory. Concept. Magn. Reson. 9:299.

Google Scholar

Rafael-Patino, J., Ramirez-Manzanares, A., Peña Acevedo, J., and Zhang, H. (2017). “Validating particle dynamics in monte carlo diffusion simulation using the finite element method,” in Proceedings of the 25th Scientific Meeting of the International Society for Magnetic Resonance in Medicine (Honolulu, HI), 1849.

Google Scholar

Raffelt, D., Tournier, J.-D., Rose, S., Ridgway, G. R., Henderson, R., Crozier, S., et al. (2012). Apparent fibre density: a novel measure for the analysis of diffusion-weighted magnetic resonance images. NeuroImage 59, 3976–3994. doi: 10.1016/j.neuroimage.2011.10.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Rensonnet, G., Scherrer, B., Girard, G., Jankovski, A., Warfield, S. K., Macq, B., et al. (2019). Towards microstructure fingerprinting: estimation of tissue properties from a dictionary of Monte Carlo diffusion MRI simulations. NeuroImage 184, 964–980. doi: 10.1016/j.neuroimage.2018.09.076

PubMed Abstract | CrossRef Full Text | Google Scholar

Rensonnet, G., Scherrer, B., Warfield, S. K., Macq, B., and Taquet, M. (2018). Assessing the validity of the approximation of diffusion-weighted-MRI signals from crossing fascicles by sums of signals from single fascicles. Magn. Reson. Med. 79, 2332–2345. doi: 10.1002/mrm.26832

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanguinetti, G., and Deriche, R. (2014). “Mapping average axon diameters under long diffusion time,” in 2014 IEEE 11th International Symposium on Biomedical Imaging (ISBI) (Beijing), 242–245.

Google Scholar

Sexton, C. E., Walhovd, K. B., Storsve, A. B., Tamnes, C. K., Westlye, L. T., Johansen-Berg, H., et al. (2014). Accelerated changes in white matter microstructure during aging: a longitudinal diffusion tensor imaging study. J. Neurosci. 46, 15425–15436. doi: 10.1523/JNEUROSCI.0203-14.2014

CrossRef Full Text | Google Scholar

Shacklock, M. O. (2007). Biomechanics of the Nervous System: Breig Revisited. Adelaide, SA: Neurodynamic Solutions.

Google Scholar

Stejskal, E. O., and Tanner, J. E. (1965). Spin diffusion measurements: spin echoes in the presence of a time-Dependent field gradient. J. Chem. Phys. 42, 288–292. doi: 10.1063/1.1695690

CrossRef Full Text | Google Scholar

Szafer, A., Zhong, J., and Gore, J. C. (1995). Theoretical model for water diffusion in tissues. Magn. Reson. Med. 33:697. doi: 10.1002/mrm.1910330516

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapp, B. D., Peterson, J., Ransohoff, R. M., Rudick, R., Mörk, S., and Bö, L. (1998). Axonal transection in the lesions of multiple sclerosis. N. Engl. J. Med. 338, 278–285. doi: 10.1056/NEJM199801293380502

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Gelderen, P., Des Pres, D., Van Zijl, P. C., and Moonen, C. T. (1994). Evaluation of restricted diffusion in cylinders. Phosphocreatine in rabbit leg muscle. J. Magn. Reson. Ser. B 103, 255–260. doi: 10.1006/jmrb.1994.1038

PubMed Abstract | CrossRef Full Text | Google Scholar

Verth, J. M. V., and Bishop, L. M. (2008). “Chapter 12 - intersection testing,” in Essential Mathematics for Games and Interactive Applications, 2nd Edn., eds J. M. V. Verth and L. M. Bishop (Boston, MA: Morgan Kaufmann), 541–599.

Google Scholar

Yeh, C.-H., Schmitt, B., Bihan, D. L., Li-Schlittgen, J.-R., Lin, C.-P., and Poupon, C. (2013). Diffusion microscopist simulator: a general monte carlo simulation system for diffusion magnetic resonance imaging. PLoS ONE 8:e76626. doi: 10.1371/journal.pone.0076626

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Hubbard, P. L., Parker, G. J. M., and Alexander, D. C. (2011). Axon diameter mapping in the presence of orientation dispersion with diffusion MRI. Neuroimage 56, 1301–1315. doi: 10.1016/j.neuroimage.2011.01.084

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Schneider, T., Wheeler-Kingshott, C. A., and Alexander, D. C. (2012). Noddi: practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage 61:1000. doi: 10.1016/j.neuroimage.2012.03.072

PubMed Abstract | CrossRef Full Text | Google Scholar