MEDUSA: a GPU-based tool to create realistic phantoms of the brain white matter microstructure using tiny spheres.
Kévin GINSBURGER1, Felix MATUSCHKE2, Fabrice POUPON3, Jean-François MANGIN3, Markus AXER2, and Cyril POUPON1

1UNIRS, CEA/Joliot/Neurospin, GIf-sur-Yvette, France, 2Research Centre Jülich, Institute of Neuroscience and Medicine, Jülich, Germany, Juelich, Germany, 3UNATI, CEA/Joliot/Neurospin, GIf-sur-Yvette, France


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.


One key objective of diffusion MRI is to understand the relationship between the diffusion-weighted signal and the underlying microscopic properties of the observed tissue. To investigate white matter, a possible approach consists in synthesizing diffusion-weighted signals for a large variety of numerical phantoms. However, state-of-the art phantom generators1-11 are mainly based on the representation of axonal fibers as cylinders, which considerably limits the achievable range of packing densities and angular dispersions. Moreover, such generators cannot create synthetic glial cells such as oligodendrocytes and astrocytes, which are prevalent in white matter and might thus have an important impact on the observed diffusion-weighted MRI signal12. This work proposes a novel Microstructure Environment Designer with Unified Sphere Atoms (MEDUSA) to tackle these issues and construct realistic white matter cellular environments 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 construction of biomimicking white matter phantoms with fully controlled geometrical properties, including axons, astrocytes and oligodendrocytes.

Material & Methods

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.

Results & Discussion

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.


MEDUSA enables the fast generation of realistic white matter phantoms with controlled geometrical properties, comprising different types of cells, while avoiding collisions thanks to a spherical decomposition of the shapes. It is an efficient tool to generate a wide variety of configurations, which can be seen as a dictionary of all possible samples of brain tissues in a given voxel. Future work will consist in widening the MEDUSA framework to a larger variety of neural cells, and creating phantoms of multiple voxels catching the structure and organization of the brain at a mesoscopic scale.


This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement No. 785907 (HBP SGA2).


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.1: Illustration of the control parameters ($$$<.>$$$: mean, $$$\sigma$$$: standard deviation ) for the fiber phantoms construction. MEDUSA creates fibers as sets of spheres. Each sphere in the scene is represented by its center, its radius, and the IDs of the cell population and item it belongs to. After the collision solver, the inner axonal membrane is created by duplicating all spheres with a radius computed using the $$$g$$$-ratio specified by the user. Ranvier nodes are modeled by setting the radius of outer membrane spheres to the value of the corresponding inner membrane sphere. Beadings are created by locally applying bell-shaped functions to the sphere radii.

Fig.2: Illustration of the astrocyte creation procedure. The top scheme shows that each astrocyte is created from a random point cloud by computing an Euclidean minimum spanning tree. A balancing factor in the distance function is used to control the level of branching. Astrocytes generated with the MEDUSA framework at 2 values of the balancing factor (BF) are shown (bottom left). The oligodendrocyte generated with MEDUSA (bottom right) has its processes connected to the myelin sheath of neighboring axons.

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.

Fig.4: Realistic scene generated with the MEDUSA framework with astrocytes (in green) and fibers (top, computation time of 12s), astrocytes, oligodendrocytes (in purple) and fibers (bottom, computation time of 15s). For illustration purposes, a small volume fraction of 0.3 was employed and 200 astrocytes were generated in both pictures. 200 oligodendrocytes were additionnally generated in the second picture. A mean diameter of $$$2.0\mu m$$$ was used for axonal fibers, and the voxel size was set to $$$100\mu m^3$$$.

Fig.5: Example phantom containing stellate neural cells at a volume fraction of 0.62 (computation time 248s), showing the potential application of the MEDUSA framework to go beyond white matter and generate gray matter phantoms.

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