This work proposes a novel Microstructure Environment Designer with Unified Sphere Atoms (MEDUSA) to construct realistic white matter cellular phantoms ensuring the absence of overlapping cells. It relies on a CUDA-based collision solver that enables to reach high values of packing density and angular dispersion, even in the case of multiple fiber populations. MEDUSA thus achieves the fast construction of biomimicking white matter phantoms with fully controlled geometrical properties, including axons, astrocytes and oligodendrocytes.
Cells construction
MEDUSA constructs axonal fibers as sets of overlapping sphere atoms (fig.1), according to the target volume fraction for each fiber population. Each population has Gamma-distributed diameters. Within a given fiber, the radii of each sphere atom can be modified and a myelin sheath, Ranvier nodes and beadings can be created (fig.1). The amount of global angular dispersion is controlled by drawing fiber orientations from a Watson distribution and local tortuosity is obtained by applying random Gaussian deformations.
Astrocytes are generated by constructing a Minimum Spanning Tree (MST) from a 3D point cloud inside a bounding sphere (fig.2), using a distance cost function controlling the amount of branching of the cell processes13. The processes are constructed as sets of overlapping spheres whose radii $$$r$$$ decrease with the distance $$$d$$$ to the soma:
$$r = r_0 e^{(-\alpha.d / R)} $$
with $$$R$$$ the total radius of the cell and $$$\alpha$$$ a tunable parameter. Similar to axons, tortuosity is induced on the processes using random Gaussian deformations.
Oligodendrocytes creation is similar to astrocytes. Additionnally, a search radius defining the area in which each oligodendrocyte myelinates axons is specified. A set of spheres composing connection points for the oligodendrocyte are selected from the outer membranes of axons, preferentially around the center of the internode regions.
Collision solver
MEDUSA starts with the generation of overlapping individual cells (axons, glial cells, neurons) from their geometrical characteristics. Each generated item belongs to a given cell type and is made of overlapping spheres, which may also overlap with spheres from other items. This set of spheres is fed into a CUDA-based collision solver which iteratively applies repulsion forces between spheres from distinct items14. The algorithm is optimized using a 1D sweep-and-prune algorithm15 that reduces the number of collision checks by sorting the spheres in the scene using their position along an arbitrary axis, enabling a 10-fold speed-up. In order to preserve the topological properties of each item, a smoothing procedure is employed for axonal fibers after each application of repulsion forces. For glial cell processes, nodal spheres corresponding to connection points between different branches or between a branch and the soma are defined. At each step, the position of all spheres between two nodal spheres are smoothed under the constraint that each branch of a given process is attached to its nodal spheres on both sides.
Fig.2 (bottom) shows that MEDUSA controls the amount of branching and creates realistic glial cells with varying radii along the processes and tortuosity. Fig.3 illustrates the ability of MEDUSA to reach a high packing density of 0.7 even with multiple fiber populations (up to 3), and shows various degrees of detail that can be achieved (presence of a myelin sheath, Ranvier nodes, local and global angular dispersion, beadings). Fig.4 shows a biomimicking virtual tissue example comprising axonal fibers, astrocytes and oligodendrocytes using realistic geometrical parameters. Fig.5 illustrates potential applications of MEDUSA to model gray matter by showing an example packing of simple stellate neural cells.
1. Hall MG, Alexander D. Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI. IEEE transactions on medical imaging, 28(9):1354–1364, 2009.
2. Yeh CH et al. Diffusion microscopist simulator: a general Monte-Carlo simulation system for diffusion magnetic resonance imaging. PLoS one, 8(10):e76626, 2013.
3. Balls GT and Frank LR. A simulation environment for diffusion weighted mr experiments in complex media. Magnetic Resonance in Medicine, 62(3):771–778, 2009.
4. Budde MD and Frank JA. Neurite beading is sufficient to decrease the apparent diffusion coefficient after ischemic stroke. Proceedings of the National Academy of Sciences, 107(32):14472–14477, 2010.
5. Fieremans E et al. Monte carlo study of a two-compartment exchange model of diffusion. NMR in Biomedicine, 23(7):711–724, 2010.
6. Harkins KD and Does MD. Simulations on the influence of myelin water in diffusion-weighted imaging. Physics in Medicine & Biology, 61(13):4729, 2016.
7. Neher PF et al. Fiberfox: facilitating the creation of realistic white matter software phantoms. Magnetic resonance in medicine, 72(5):1460–1470, 2014.
8. Landman BA et al. Complex geometric models of diffusion and relaxation in healthy and damaged white matter. NMR in Biomedicine, 23(2):152–162, 2010.
9. Lin M et al. Simulation of changes in diffusion related to different pathologies at cellular level after traumatic brain injury. Magnetic resonance in medicine, 76(1):290–300, 2016.
10. Rensonnet G et al. Towards microstructure fingerprinting: Estimation of tissue properties from a dictionary of monte carlo diffusion mri simulations. NeuroImage, 2018.
11. 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. Frontiers in Physics,6:12, 2018.
12. Palombo M et al. Can we detect the effect of spines and leaflets on the diffusion of brain intracellular metabolites? NeuroImage, 2017.
13. Cuntz H et al. One rule to grow them all: a general theory of neuronal branching and its practical application. PLoS computational biology, 6(8):e1000877,2010.
14. Altendorf H and Jeulin D. Random-walk-based stochastic modeling of three-dimensional fiber systems. Physical Review E, 83(4):041804, 2011.
15. Avril Q et al. A broad phase collision detection algorithm adapted to multi-core architectures. In Vric 2010 Proceedings, volume 12, page 95, 2010.
Fig.3: Top: Example phantoms containing 1, 2 and 3 fiber populations at a volume fraction of 0.7. A mean diameter of $$$2.0\mu m$$$ was employed for each population. The voxel size is $$$100 \mu m^3$$$. The MEDUSA computation time on an Nvidia Quadro M1000M GPU was equal to 113s, 125s and 144s respectively for each phantom.
Bottom: Axonal fibers (1 population) generated with MEDUSA with different microstructural details (LAD=Local angular dispersion/GAD=Global angular dispersion). Left: tortuous axons with myelin sheath and Ranvier nodes. Middle: 10 degrees of GAD are added according to a Watson distribution. Right: axonal beadings (local swelling of the axonal membrane) are rendered.