Ross Callaghan^{1}, Noam Shemesh^{2}, Daniel Alexander^{1}, Hui Zhang^{1}, and Marco Palombo^{1}

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.

*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 simulators^{3,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 Camino^{3 }using 500,000 diffusing spins uniformly distributed with bulk-diffusivity D_{0}=2µm^{2}/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 water^{5,6}.

To this end, 300 fibres were grown:

- with diameter d
_{0 }drawn from a gamma distribution with mean d_{0}=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; - with variable fibre diameter d, gamma distributed with mean d
_{0}=1 µm and σ=0.5d_{0},1d_{0},2d_{0 }and no undulation.

The diffusion of 500,000 spins distributed proportionally to the volume of each fibre was simulated using Camino^{3}. From the spins’ trajectories, the axial diffusivity D_{z }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 values^{7}. 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 D_{0}=2µm^{2}/ms) while beading a larger effect at long t (t>25ms, for D_{0}=2µm^{2}/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.

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 μm^{3} 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/μm^{2} 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, from^{7 }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 D_{z} at values lower than D_{0}. At short t undulation leads to time-dependence. b) At long t, increased beading of fibres leads to time-dependence in the D_{z}. 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/µm^{2}, D=6.85,15.85,25.85,50.85,75.85,100.85 ms, 15 gradient directions, resolution 100x100x250 mm^{3}, 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)..