Towards a more realistic and flexible white matter numerical phantom generator for diffusion MRI simulation
Ross Callaghan1, Noam Shemesh2, Daniel Alexander1, Hui Zhang1, and Marco Palombo1

1Centre for Medical Image Computing, Department of Computer Science, University College London, London, United Kingdom, 2Champalimaud Research, Champalimaud Centre for the Unknown, Lisbon, Portugal


We present a novel method for the generation of realistic white matter numerical phantoms with controllable morphology, and more realistic orientation-dispersion and packing-density. A shift of paradigm is proposed: rather than ‘packing fibres’, our algorithm ‘grows fibres’ contextually and efficiently, avoiding intersection between fibres. The potential of the method is demonstrated by reaching the highest orientation-dispersion and packing-density ever. The flexibility of the method is demonstrated with simulations of diffusion-weighted MR signal in three example substrates with differing orientation-dispersions, packing-densities and permeabilities. A proof-of-concept application is presented to investigate the impacts of undulation and beading on the intra-axonal axial diffusivity.


The aim of this work is to present a novel method for generating realistic white matter (WM) numerical phantoms, with controllable morphology and more realistic orientation dispersion and packing-density. Generating realistic WM numerical phantoms (dispersion/undulation/beading/etc.) at high packing-densities is a major open challenge for the diffusion MRI community. While densely packing straight, parallel, fibres is relatively easy, only a few groups have attempted to densely pack irregular, non-parallel, fibres1,2, with the highest packing-density achieved under modest dispersion reaching only 20%. Here, we propose a completely different method: rather than densely ‘packing’ irregular fibres, we ‘grow’ fibres contextually, mimicking natural fibre genesis. We provide proof-of-concept demonstration and example of application of interest for DW-MRI.


Fibre-Growth Algorithm

The fibre-growth algorithm grows fibres one-by-one on a fully connected network of nodes avoiding intersection between fibres and is illustrated in Fig.1. Each fibre has a specified starting point and target point towards which it will grow. Fibres grow along the edges of a fully connected 3D network defined by the Delaunay triangulation of a set of randomly distributed nodes. This acts to limit the number of paths tested for a given fibre, making this approach faster than brute force growth of fibres. The connected network holds memory of the growth path of each fibre meaning collision checks are computationally efficient. Then, a meshing process in Blender (https://blender.org), creates 3D meshes compatible with available DW-MRI simulators3,4.

Demonstration of the Fibre-Growth Algorithm

To demonstrate the potential of the method, three substrates at different (dispersion, packing-density) conditions were generated: (0°,60%), (15°,30%); (35°,25%) (Fig.2). For each substrate, the Pulsed-Gradient-Spin-Echo (PGSE) signal was simulated in Camino3 using 500,000 diffusing spins uniformly distributed with bulk-diffusivity D0=2µm2/ms. To show the range of simulation possibilities available, three different membrane permeabilities (κ=0,0.0025,0.0050 µm/ms) were also imposed.

Application: Undulation vs. Beading

As a proof-of-concept application, we use our algorithm to investigate the specific impact of axonal undulation and axonal size variability, namely beading, on the diffusion-time t dependence of water diffusivity and whether it would be possible to disentangle the two effects. Indeed, it has been shown that both axonal beading and undulation have an effect on the t-dependence of the axial apparent diffusion coefficient (ADC||) of intra-axonal water5,6.

To this end, 300 fibres were grown:

  1. with diameter d0 drawn from a gamma distribution with mean d0=1 µm and standard deviation σ=0.2 µm fixed along the fibre and variable undulation amplitude, a=2,4,8 µm and wavelength λ=10,50,100 µm;
  2. with variable fibre diameter d, gamma distributed with mean d0=1 µm and σ=0.5d0,1d0,2d0 and no undulation.

The diffusion of 500,000 spins distributed proportionally to the volume of each fibre was simulated using Camino3. From the spins’ trajectories, the axial diffusivity Dz was computed for t=[0:0.1:500]ms


Fibre-Growth Algorithm

The versatility of the algorithm is shown in Fig.2a,with substrates with a range of (dispersion, packing-density) conditions (0°,60%), (15°,30%); (35°,25%). The corresponding simulated PGSE signals at different permeabilities (κ=0,0.0025,0.0050 µm/ms) are shown in Fig.2b&c.

Application: Undulation vs. Beading

Fig.3 shows some examples of the fibres used in this application, together with comparison with published histological values7. The t-dependent axial diffusivity in the undulating vs. beading substrates are reported in Fig.4. Results suggest that the effects of undulation and beading impact the t-dependence of the intra-axonal ADC|| in different diffusion time regimes: undulation has a larger effect at short t (t<5ms, for D0=2µm2/ms) while beading a larger effect at long t (t>25ms, for D0=2µm2/ms). The experimental data in Fig.4c follow a combination of undulation and beading, suggesting that at t>25ms, the measured t-dependence of ADC|| may be mostly driven by beading, rather than undulation or dispersion.

Discussion and Conclusion

Here we present a novel method for generating more realistic synthetic WM substrates for diffusion simulation based on a shift of paradigm: from ‘packing fibres’ to ‘growing fibres’. The framework is flexible, allowing for growth of a range of controllable fibre configurations. The proposed algorithm represents a major step towards very high fibre packing, enabling us to reach the highest dispersion at the highest packing-density ever, to our knowledge. Our (15°,30%) and (35°,25%) represent an average ~50% and ~200% improvement, respectively, over the best previously reported results of (10°,20%)2.

A proof-of-concept experiment is proposed to demonstrate the usefulness of the new algorithm in investigating the specific contribution of selected microstructural features to DW-MRI metrics. We show, for the first time to our knowledge, that the effects of undulation and beading impact the ADC in different diffusion time regimes offering the potential to disentangle the two effects.


This work is supported by the EPSRC-funded UCL Centre for Doctoral Training in Medical Imaging (EP/L016478/1) and the Department of Health’s NIHR-funded Biomedical Research Centre at UCLH.

This work is supported by EPSRC (EP/G007748, EP/I027084/01, EP/L022680/1, EP/M020533/1, N018702), EPSRC EP/M507970/1 and European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Starting Grant, agreement No. 679058).


1. Rafael-Patino, J. et al. Realistic 3D Fiber Crossing Phantom Models for Monte Carlo Diffusion Simulations. in 26th Annual Meeting of the International Society for Magnetic Resonance in Medicine (ISMRM)(2018).

2. Ginsburger, K. et al. 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.5,1–18 (2018).

3. Hall, M. G. & Alexander, D. C. Convergence and Parameter Choice for Monte-Carlo Simulations of Diffusion MRI. IEEE Trans. Med. Imaging 28,1354–1364 (2009).

4. Yeh, C.-H. et al. Diffusion Microscopist Simulator: A General Monte Carlo Simulation System for Diffusion Magnetic Resonance Imaging. PLoS One8,e76626 (2013).

5. Brabec, J., Lasic, S. & Nilsson, M. Simulation of diffusion in axons with harmonic and stochastic trajectories. in 26th Annual Meeting of the International Society for Magnetic Resonance in Medicine (ISMRM)(2018).

6. Lee, H.-H., Novikov, D. S. & Fieremans, E. Monte Carlo simulations of diffusion in mouse corpus callosum reconstructed from 3-d electron microscopy validate the time-dependence along axons. in 26th Annual Meeting of the International Society for Magnetic Resonance in Medicine (ISMRM)(2018).

7. Lee, H.-H. et al. Electron microscopy 3-dimensional segmentation and quantification of axonal dispersion and diameter distribution in mouse brain corpus callosum. bioRxiv, 357491,(2018).


The growth algorithm. a) from the current node (highlighted in green), candidate nodes for the next step are found as those sharing an edge. b) For each edge a cost is calculated based on how well aligned with the target that node is as well as how much the diameter of the fibre would have to shrink at that node. c) The lowest cost node is chosen and a fibre section grown between that and the current node. d) Node weights are calculated based on the distance from the new fibre section. These weights will then control the shrinkage of future fibres.

a) Example substrates (cut into 30x30x30 μm3 cube) from the fibre growth algorithm, left to right: Zero macroscopic dispersion (60% density), 15° macroscopic dispersion (30% density), 35° dispersed (25% density). b) Simulations for each substrate for varying permeabilities with SNR = $$$\infty$$$ and c) SNR = 20. Simulated PGSE measurement parameters: /∆=1/40 ms and 50 b-values from 0 to 9 ms/μm2 along x-, y- and z-directions. Signal is averaged across directions.

Examples of individual fibres used in simulation of undulation vs variable radius. a) zero undulation, variable diameter b) undulation with fixed diameter and c) both undulation and variable radius. d) and e) Measured axonal diameter histogram and distribution along a number of fibres respectively, from7 f) and g) Diameter histogram and distribution from a number of fibres from the examples in a).

Simulation and experimental results. a) At long t, undulation leads to almost plateau of Dz at values lower than D0. At short t undulation leads to time-dependence. b) At long t, increased beading of fibres leads to time-dependence in the Dz. c) Measured Pulsed-Gradients-Stimulated-Echo data from ex-vivo mouse corpus-callosum follow a mix of undulation and beading. Experiments at 16.4 T (scanner Bruker-BioSpin, Germany) with: TE/TR=46/2500 ms; d=2.5 ms, b=1.5,3 ms/µm2, D=6.85,15.85,25.85,50.85,75.85,100.85 ms, 15 gradient directions, resolution 100x100x250 mm3, 4 averages. The first eigenvalues of the diffusion tensor, ADC||, at each diffusion time Δ-δ/3 were computed using FSL (https://fsl.fmrib.ox.ac.uk/fsl)..

Proc. Intl. Soc. Mag. Reson. Med. 27 (2019)