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

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)

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)

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

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):

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:

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.

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:

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

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.

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:

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.

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.

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.

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.

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.

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.

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

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