Phantoms for Diffusion Simulations: Multi-Objective Differential Evolution for Realistic Numerical (MODERN) Phantoms.
Jonathan Rafael-Patino1, Thomas Yu1, Mariam Andersson2,3, Hans Martin Kjer2,3, Vedrana Andersen Dahl3, Alexandra Pacureanu4, Anders Bjorholm Dahl3, Tim B. Dyrby2,3, and Jean-Philippe Thiran1,5

1Signal Processing Lab (LTS5), École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, 2Danish Research Centre for Magnetic Resonance, University Hospital Hvidovre, Hvidovre, Denmark, 3Department of Applied Mathematics and Computer Science, Technical University of Denmark, Kongens Lyngby, Denmark, 4European Synchrotron Radiation Facility, Genoble, France, 5Radiology Department, Centre Hospitalier Universitaire Vaudois and University of Lausanne, Lausanne, Switzerland


The following work presents a novel, generally applicable framework (MODERN) which takes, as input, statistics from histology and automatically generates biologically plausible numerical phantoms whose morphological features are optimized to match the input statistics. As a proof of concept, MODERN is used to generate three-dimensional geometrical meshes representing a bundle configuration of axons which can be immediately integrated into Monte Carlo simulations for diffusion MRI acquisitions. The obtained features of the optimised phantoms are shown to match those from the input values. The statistics used were obtained from axonal segmentations from synchrotron imaging.


Monte-Carlo Diffusion Simulations have proved to be a powerful approach to study Diffusion-Weighted Magnetic Resonance Imaging, from generating simulated ground-truth data to studying the diffusion process in complex media. However, current numerical phantoms are based on simple, 3D geometries which can lead to biased simulation results[1,2]. We present a novel framework for automatically generating 3D numerical phantoms for use in Monte-Carlo simulations. The framework takes as input arbitrary parameters encoding desired morphological features such as the mean curvature and uses multi-objective optimization with differential-evolution(DE) to generate a phantom whose features match the input parameters. To generate realistic numerical phantoms, the framework(MODERN) takes as input a list of statistical parameters from 3D-histology and iteratively deforms an initial configuration by minimizing, using DE, a multiobjective cost function. The cost function penalizes the mismatch between the phantom’s features and the inputs. As a proof of concept, we present a phantom with a realistic axon bundle configuration. We propose to use three statistics taken from 3D-histology of axons using synchrotron imaging: tortuosity and local/global mean angular deviation. Realistic axon phantoms for diffusion-MRI are especially important since novel methods require ground-truths for validation which mimic the actual tissue properties[3,4,5].


Samples from the perfusion fixed splenium of a 32-month-old female Vervet monkey brain were stained with 0.5% osmium tetroxide(OsO4) to give contrast to the myelin. These were embedded in EPON and imaged with x-ray phase contrast tomography(PCT) at beamline ID16A of the European Synchrotron Research Facility(ESRF) in Grenoble, France. The PCI experiments produced 75nm isotropic resolution 3D-volumes, with a field-of-view(FOV) of 150µm in all dimensions. Fifty-four axons were semi-manually segmented using ITK-Snap[6], and the axonal centerlines were generated using a layered surface segmentation method developed at the DTU and the DRCMR. In MODERN, we utilize three statistics commonly used in characterization of vascularities(which can similarly be applied to axons): the tortuosity$$$(tt)$$$ and the local/global mean angular deviations$$$(\theta^l$$$and $$$\theta^g$$$resp.). Each metric was computed using the centerlines of each segmented axon. The metrics are defined as follows:$$tt=\sum_{i=0}^{\#S}\frac{|S_{\#S}-S_1|}{|S_i-S_{i-1}|}$$ $$\theta^I=\sum_{i=2}^{\#S-1}\frac{ang((S_{i+}-S_i ), (S_i -S_{i-1}))}{\#S}$$ $$\theta^g=\sum_{i=2}^{\#S}\frac{ang((S_{\#S}-S_1), (S_i-S_{i-1}))}{\#S}$$

Where$$$ang(v,w)$$$ computes the angle between two vectors,$$$\#S$$$ is the number of segments, and$$$S_i$$$ denotes the i-th point of a given axon skeleton. MODERN steps: The workflow of MODERN requires four steps as is shown in Figure1.The input parameters are the desired values for the ICVF, the number of cylinders to generate, the tortuosity, and the mean local/global deviations. 1)Gamma distributed radii fitting. Taking as input the number of cylinders, the parameters of a gamma distribution ($$$\alpha$$$ $$$\beta$$$), and the desired intra-cellular volume fraction(ICVF), a set of cylinders is sampled and positioned inside a circular boundary to satisfy the desired ICVF until a user set stopping criterion is reached(Figure1-a). 2)Initial bundle configuration. An initial configuration of regular, undulating trajectories is chosen such that the obtained tortuosity is close to the objective value (Figure1-b). 3)DE Optimization. The full set of trajectories is then optimized using a multi-objective cost function defined as follows:$$E({\cup}S)=\lambda_1(tt'-tt({\cup}S))^2+\lambda_2(\theta^{l'}-\theta^{l}({\cup}S))^2+\lambda_3(\theta^{g'}-\theta^{g}({\cup}S))^2+\lambda_4{\tau}({\cup}S,lr)$$

Where $$${\cup}S$$$ is the full set of generated trajectories,$$$tt'$$$,$$$\theta^{l'}$$$,$$$\theta^{g'}$$$ are the desired (input) values to achieve, and$$$tt(\cup S)$$$,$$$\theta^{l}({\cup}S),$$$,$$$\theta^{g'}({\cup}S)$$$ are the mean values of the descriptors computed from the full set of trajectories. Finally, $$$\tau({\cup}S,lr)$$$ is a penalization term that checks that no axon control point is outside a limit radius ($$$lr$$$) defined by the initial configuration and that there are no intersections between trajectories. The objective function is then optimized using a bespoke implementation of differential evolution which saves any set of parameters that fulfil the no overlapping condition$$$(\tau=0)$$$ and is close to the objective value by a given threshold(Figure1-c). 4)Post-processing and added perturbations. The resulting trajectories are meshed and decimated. Optionally, local perturbations can be added to each strand to add local complexity along the trajectory, in contrast with the usual cylindrical representation of the axons(Figure1-d).

Results and Discussion

Figure2 shows the distribution of features from the segmentation. The mean for each feature: tortuosity, local angular deviation and global deviation were$$$1.07$$$,$$$13.02°$$$, and$$$13.39°$$$ respectively. The previous values of each feature were used as an objective value for the optimization. Figure3 shows the resulting full bundle with ~5000 axons, along with the feature distribution for the set of individual axons. The mean obtained parameters for the tortuosity, local angular deviation and global deviation are$$$1.058$$$,$$$13.68°$$$, and$$$14.68°$$$ respectively.


In this work, we present an in-silico framework(MODERN) to automatically generate realistic, numerical phantoms for axon bundles with a single, general orientation and angular dispersion in the micro/meso scale. Our results show the capacity of MODERN to generate complex, 3D geometrical configurations which match a set of feature descriptors obtained from histological data.


Capital Region Research Foundation (grant number: A5657)

Thomas Yu is supported by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie project TRABIT (agreement No 765148).


1. Hall M, Alexander DC. Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI data. (2009).

2. Novikov, Dmitry S. and Fieremans, Els and Jensen, Jens H. and Helpern, Joseph A. Random walks with barriers. Nature Publishing Group. (2011)

3. Nedjati-Gilani et. al. Machine learning based compartment models with permeability for white matter (2017)

4. Close, Thomas G. and Tournier, Jacques Donald and Calamante, Fernando and Johnston, Leigh A. and Mareels, Iven and Connelly, Alan. A software tool to generate simulated white matter structures for the assessment of fibre-tracking algorithm. (2009)

5. Ginsburger, Kevin and Poupon, Fabrice and Beaujoin, Justine and Estournet, Delphine and Matuschke, Felix and Mangin, Jean-Francois and Axer, Markus and Poupon, Cyril. Improving the Realism of White Matter Numerical Phantoms: A Step toward a Better Understanding of the Influence of Structural Disorders in Diffusion MRI. (2018)

6. Yushkevich, Paul A. and Piven, Joseph and Hazlett, Heather Cody and Smith, Rachel Gimpel and Ho, Sean and Gee, James C. and Gerig, Guido. User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability. (2006)


Figure1: The four steps of the optimization framework. a) The plot of the radii configuration inside the bounding circle. b) Illustration of a regularly undulating configuration used as the initial configuration. Each undulated cylinder is parametrized as 3D set of control points $$$p_i$$$, the optimization procedure perturbates the segments between two points $$$S_i=(p_i - p_{i-1})$$$. c) Example of a generated axonal configuration after the optimization procedure. Each segment has been perturbated individually to match the desired global metrics. d) Example of the post-processing used to add irregularities to the configuration.

Figure 2: Top panel shows a set of centre lines used to compute the tortuosity and micro/meso angular deviations. For more information refer to Mariam Andersson et. al., ISMRM-2019 submitted abstract: "Uncovering 3D Axonal Morphologies with Synchrotron Imaging: Impact on Microstructure Imaging with Diffusion MRI". The lower panel shows the obtained distribution of the computed metrics of the 3D-histological data. The resulting mean tortuosity, local angular deviation and global deviation are 1.07, 13.02 degrees, and 13.39 degrees respectively.

Figure 3: Top: Result of optimization with a total of 4952 axons. Each axon is parametrized using a set of control points to define the connected segments between them. The Below: Histograms of the tortuosity and angular deviations (blue for the mean global angular deviation, and red for the local) over the axons. The resulting mean tortuosity, local angular deviation and global deviation are 1.058, 13.68 degrees, and 14.68 degrees respectively.

Figure 4. View of the axon bundle which is bisected along the central line for visualization of the internal structure. In this view we can see a relatively homogeneous configuration; this is a result of the boundaries and the high packing density of 0.54 ICVF.

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