Jonathan Rafael-Patino^{1}, Thomas Yu^{1}, Mariam Andersson^{2,3}, Hans Martin Kjer^{2,3}, Vedrana Andersen Dahl^{3}, Alexandra Pacureanu^{4}, Anders Bjorholm Dahl^{3}, Tim B. Dyrby^{2,3}, and Jean-Philippe Thiran^{1,5}

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.

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).

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.