Front. Neuroinform. Frontiers in Neuroinformatics Front. Neuroinform. 1662-5196 Frontiers Media S.A. 10.3389/fninf.2023.1208073 Neuroscience Methods CACTUS: a computational framework for generating realistic white matter microstructure substrates Villarreal-Haro Juan Luis 1 * Gardier Remy 1 Canales-Rodríguez Erick J. 1 Fischi-Gomez Elda 1 2 3 Girard Gabriel 1 2 3 4 Thiran Jean-Philippe 1 2 3 Rafael-Patiño Jonathan 1 3 1Signal Processing Laboratory (LTS5), École Polytechnique Frale de Lausanne (EPFL), Lausanne, Switzerland 2CIBM Center for Biomedical Imaging, Lausanne, Switzerland 3Radiology Department, Centre Hospitalier Universitaire Vaudois, University of Lausanne, Lausanne, Switzerland 4Department of Computer Science, University of Sherbrooke, Sherbrooke, QC, Canada

Edited by: Ludovico Minati, University of Electronic Science and Technology of China, China

Reviewed by: Viktor Vegh, The University of Queensland, Australia; Andrada Ianus, Champalimaud Foundation, Portugal

 *Correspondence: Juan Luis Villarreal-Haro juan.villarrealharo@epfl.ch 01 08 2023 2023 17 1208073 18 04 2023 13 07 2023 Copyright © 2023 Villarreal-Haro, Gardier, Canales-Rodríguez, Fischi-Gomez, Girard, Thiran and Rafael-Patiño. 2023 Villarreal-Haro, Gardier, Canales-Rodríguez, Fischi-Gomez, Girard, Thiran and Rafael-Patiño

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Monte-Carlo diffusion simulations are a powerful tool for validating tissue microstructure models by generating synthetic diffusion-weighted magnetic resonance images (DW-MRI) in controlled environments. This is fundamental for understanding the link between micrometre-scale tissue properties and DW-MRI signals measured at the millimetre-scale, optimizing acquisition protocols to target microstructure properties of interest, and exploring the robustness and accuracy of estimation methods. However, accurate simulations require substrates that reflect the main microstructural features of the studied tissue. To address this challenge, we introduce a novel computational workflow, CACTUS (Computational Axonal Configurator for Tailored and Ultradense Substrates), for generating synthetic white matter substrates. Our approach allows constructing substrates with higher packing density than existing methods, up to 95% intra-axonal volume fraction, and larger voxel sizes of up to 500μm3 with rich fibre complexity. CACTUS generates bundles with angular dispersion, bundle crossings, and variations along the fibres of their inner and outer radii and g-ratio. We achieve this by introducing a novel global cost function and a fibre radial growth approach that allows substrates to match predefined targeted characteristics and mirror those reported in histological studies. CACTUS improves the development of complex synthetic substrates, paving the way for future applications in microstructure imaging.

 microstructure imaging diffusion MRI brain imaging white matter Monte-Carlo simulations numerical phantom synthetic substrates high packing density 1. Introduction

Diffusion-weighted magnetic resonance imaging (DW-MRI) is a non-invasive technique used to study the microscopic structure of biological tissues in vivo. It is sensitive to the ensemble of water molecules (wherein each molecule follows a random motion pattern) as they interact with cellular surfaces (Simpson and Carr, 1958; Stejskal and Tanner, 1965; Bihan, 1995). This technique provides a valuable tool to study brain microstructure and its alterations following injury (Parizel et al., 2005; To et al., 2022) and neurological disease (van Gelderen et al., 1994; Budde and Frank, 2010; Narvaez-Delgado et al., 2019).

White matter is a crucial component of the brain, composed of highly organized axon bundles that interconnect cortical regions and subcortical regions (Brückner et al., 1996; Sporns, 2011). Various imaging techniques have been considered to characterise the white matter tissue microstructure in different species. For example, axon diameters have been measured in some white matter regions of the macaque monkey brain using histology and DW-MRI (Caminiti et al., 2013), and optical microscopy (Innocenti and Caminiti, 2017). These studies show that the estimated distribution of axon diameters is long-tailed, with a mean of around one micrometre. Recent studies that used high-resolution three-dimensional (3D) synchrotron X-ray nano-holotomography (Andersson et al., 2020) and 3D electron microscopy (Lee et al., 2019) found that axons are non-cylindrical and exhibit environment-dependent variations in diameter and trajectory. Alongside axon diameters, another relevant feature is the intracellular volume the axons occupy in a predetermined region. In histological postmortem data, the white matter intracellular space volume has been estimated as ranging between 60 and 85% of the brain volume for macaques (Stikov et al., 2015) and human adults (Syková and Nicholson, 2008). Interestingly, it goes as high as 70–95% in mice, as reported by light microscopy (Tønnesen et al., 2018), and cryo and chemical fixations (Korogod et al., 2015). It is speculated that this range might be influenced by the shrinkage of the extracellular space due to the fixation process (Dam, 1979; Bolduan et al., 2020).

Given the importance of studying white matter tissue microstructure in vivo, several DW-MRI models have been proposed (e.g., Murday and Cotts, 1968; Neuman, 1974; van Gelderen et al., 1994; Söderman and Jönsson, 1995; Stanisz et al., 1997; Assaf et al., 2004, 2008; Assaf and Basser, 2005; Alexander et al., 2010; Dyrby et al., 2011; Drobnjak et al., 2016; Jelescu and Budde, 2017; Kakkar et al., 2018; Novikov et al., 2018, 2019; Lee et al., 2020; Veraart et al., 2020, 2021; Harkins et al., 2021). However, validating these non-invasive techniques requires physical and numerical phantoms with a well-known microstructure (Campbell et al., 2005; Fieremans et al., 2008; Tournier et al., 2008; Fillard et al., 2011; Lavdas et al., 2013; Maier-Hein et al., 2017; Zhou et al., 2018; Schilling et al., 2019; Andersson et al., 2020; Lee et al., 2020; Rafael-Patino et al., 2020; Warner et al., 2023). Phantoms, in the context of this paper, are geometrical models of brain tissue structures that serve as a proxy or reference for evaluating the performance of imaging techniques. While physical phantoms have been widely used, they are often limited by their high costs and the impracticality of replicating axons' sizes and complex spatial arrangement. Therefore, numerical phantoms have emerged as the most popular validation technique for studying the complexities of diffusion phenomena in cases where analytical solutions are unavailable; because they only require a substrate that mimics the tissue of interest to simulate the displacements of water molecules and corresponding DW-MRI signal (Close et al., 2009; Côté et al., 2013; Neher et al., 2014). Nevertheless, the difficulty in Monte-Carlo simulations lies in accurately mimicking the geometry of white matter tissue (Hall and Alexander, 2009; Nilsson et al., 2012, 2017; Baxter and Frank, 2013; Plante and Cucinotta, 2013; Grussu et al., 2019; Truffet et al., 2020).

Various studies have proposed to generate numerical phantoms approaching the tissue's morphological complexity and density. For instance, two popular tools, MEDUSA (Ginsburger et al., 2019) and CONFIG (Callaghan et al., 2020), focus on generating specialized voxel-wise phantoms with microstructural geometries that replicate the properties of white matter. Recently, a tailored modification of Close et al. (2009) framework was used to build challenging substrates for the DiSCo challenge (Rafael-Patino et al., 2021), aimed to test fibre-tracking and connectivity methods on large-scale synthetic datasets from DW-MRI Monte-Carlo simulations. While these methods have provided valuable tools to characterise and simulate DW-MRI signals in numerical substrates, they still have important limitations regarding the maximum packing density and substrate size achieved. For instance, state-of-the-art frameworks can generate synthetic substrates with packing densities up to 75% (Ginsburger et al., 2019; Callaghan et al., 2020; Rafael-Patino et al., 2021), whereas the density found in histological data goes up to 95% in some regions (Korogod et al., 2015; Tønnesen et al., 2018). Moreover, they cannot sample substrate beyond 100μm3, which in turn restricts the sampling diversity achieved for morphological features (Romascano et al., 2018; Rafael-Patino et al., 2020). Therefore, the DW-MRI signals generated from these substrates may not accurately mimic the brain signals measured in white matter regions with higher packing densities.

To overcome these limitations, we introduce a novel computational workflow, CACTUS (Computational Axonal Configurator for Tailored and Ultradense Substrates), to generate synthetic fibres with rich microstructure characteristics. Expanding on previous methods (Close et al., 2009; Ginsburger et al., 2019; Rafael-Patino et al., 2020), we develop a novel numerical phantom generator for white matter substrates. CACTUS solves the high-density packing problem and achieves up to 95% intracellular volume fractions while efficiently generating substrate sizes up to 500μm3. Furthermore, CACTUS is highly customisable, capable of generating synthetic substrates with a wide range of characteristics, such as single-bundle (Stikov et al., 2015), bundle crossings (Tuch, 2004; Tournier et al., 2007; Schilling et al., 2017; Canales-Rodr-guez et al., 2019), orientation dispersion (Zhang et al., 2012; Daducci et al., 2015), gamma-distributed axon radii (Assaf et al., 2008; Sepehrband et al., 2016), non-constant longitudinal fibre-radii (Andersson et al., 2020), substrates with non-cylindrical fibres and tortuous surfaces (Lee et al., 2019), and myelin compartments (Mackay et al., 1994; Stikov et al., 2015; Canales-Rodr-guez et al., 2021). Through these features, CACTUS expands on the capabilities of existing substrate generation methods, providing a flexible and versatile tool for studying white matter microstructure in controlled environments.

 2. Methods

CACTUS generates synthetic substrates in three steps (see Figure 1): a) substrate initialisation, b) joint fibre optimisation, c) fibre radial growth (FRG).

Example of the CACTUS method steps to create a synthetic substrate. (A) The substrate initialization orients fibres in a bundle to achieve a predefined mean angular dispersion (e.g., 15°). (B) The joint fibre optimization step removes fibre overlaps by adapting their trajectories and local radii. In the case of bundle crossings, the trajectories are trimmed to the centre of the crossing. (C) The fibre radial growth step further increases the fibre-packing density while keeping the predefined target radius distribution.

Firstly, in the substrate initialization step, synthetic straight cylindrical fibres are initialized and parameterized inside a cuboid. In CACTUS, a single fibre population (bundle) is a group of fibres arranged cohesively along one main orientation. A bundle has two main properties: the average global dispersion, which is the mean angle between the main orientation of each fibre and the bundle, and the target radii distribution, from which the fibre radii are sampled.

In the second step, the joint fibre optimization, CACTUS extends previously proposed frameworks (Close et al., 2009; Ginsburger et al., 2019) based on local optimization. In our case, we aim to minimise a cost function that penalises some essential fibre properties such as overlapping, high curvature, increase in length, and promote compactness. Moreover, CACTUS introduces a new fibre parameterization based on capsules, which reduces the number of parameters needed to characterise fibre trajectories and handles fibre overlapping more efficiently. The resulting optimization problem is solved via a gradient descent algorithm (Duchi et al., 2011). During optimization, CACTUS prioritises removing fibre overlapping, while the penalization of curvature, length and promotion of compactness maintains a coherent fibre structure at all time-points.

Finally, the fibre trajectories are used to mesh the fibre surfaces in the fibre radial growth (FRG) step. The FRG also increases the packing density while keeping the correspondent fibre's parameterization structure using a discrete grid to seed, to grow, and to rearrange the fibre into the final substrates. The grid discretization defines the fibres' isosurface needed to compute the final surfaces with a marching cube algorithm (Lewiner et al., 2003).

2.1. Substrate initialization

Our substrate initialization algorithm enhances the circle two-dimensional (2D) packing algorithm proposed by Hall and Alexander (2009) to create a 3D packing of bundles. The algorithm creates a single bundle by initializing the fibres inside a cuboid of dimensions L × L × H. The endpoints of the fibres are contained within the L × L squared faces, while the orientation of the cuboid's height H and the bundle are aligned to the Z-axis. The algorithm packs 2D circles in the opposite faces of the cuboid, sampling radii from a gamma Γ(α, β) distribution (Assaf et al., 2008; Sepehrband et al., 2016), until the target density is met. At the same time, the algorithm packs the two opposite 2D circles to create an initialization such that the bundle reach the specified mean angular dispersion η. In scenarios where the target density exceeds 75%, an adjustment is made by shrinking the radii of the distribution. This allows the algorithm to continue packing until reaching a density of 75%. It is important to note that the radii will subsequently grow non-uniformly back to their original size during the execution of the algorithm while simultaneously achieving the desired final target density. In order to create a substrate with two bundles crossing at an inter-bundle angle of θ, two different bundles are initialized in their respective cuboids and subsequently rotated and translated are applied. Figure 1 shows examples of a single bundle and a bundle crossing initialization. Finally, we parameterise each fibre's skeleton as the trajectory of its centre of mass. This trajectory is defined by several control points connecting the two endpoints sampled during the packing algorithm, where each point has a corresponding radius.

 2.2. Joint fibre optimization

Once the substrate is initiated as described in Section 2.1, fibres may overlap. CACTUS employs an optimization method to readjust the fibre trajectories and disentangle overlaps by defining several cost functions. These cost functions, inspired by Close et al. (2009) and Ginsburger et al. (2019), help to regularise and obtain coherent fibre structures with the specified target properties. Ordered by priority of penalization, these target properties are as follows: (i) fibre overlapping (see Section 2.2.1), (ii) high curvature, (iii) increased fibre length, (iv) changes in radii, and (v) compactness. The optimization algorithm alternative between two steps: first minimises the overlapping cost function, then the subsequent step aimed at minimizing the remaining cost functions. An algorithm requirement is to identify a solution that exhibits no overlaps. Once a solution without overlaps is achieved, the algorithm iterates further to reduce (when possible) the penalization associated with the remaining cost functions while maintaining the absence of overlaps.

In the following subsection, we introduce the novel parameterization and overlapping cost function based on capsules, which is a key contribution of our work. As the remaining cost functions are relatively straightforward and similar to those in previous studies, we have provided their definitions in Supplementary material (section joint fibre optimization).

2.2.1. Fibre capsule-parameterization

Fibres are parameterized as skeletons made of 3D control points. In the overlapping cost function, every pair of consecutive points in the skeleton forms a capsule, defined with the set of parameters [p0, p1, r0, r1], where p0, p1 ∈ ℝ3 are the initial/ending points of the capsule, and r0, r1 ∈ ℝ are their respective radius (see Figure 2A). In our scenario, the length of a capsule (distance form pi to pi+1) is not restricted, but we suggest the ranges between 12ri up to 2ri, and the change sampling frequency increases as the radii decrease.

(A) Capsule example, whose parameters are two points and a radius at each control point. (B) Example of fibre as a chain of capsules. Two adjacent capsules share a control point (position and radius parameters). (C) Illustrative example of the key components within the overlapping cost function. It showcases the intersection of two capsules, with two spheres representing the intersection region between the capsules.

In this framework, a fibre parameterization can be defined as a chain of capsules (see Figure 2B). The fibre Sa, with ma control points, is composed of the capsules determined by the subsequent point pairs as [xia,xi+1a,ria,ri+1a] with control points {x0a,x1a,…xma−1a}⊂ℝ3 and associated radius {r0a,r1a,…rma−1a}⊂ℝ.

 2.2.2. Overlapping cost function

The overlapping cost function handles the fibre collision by identifying overlaps from two capsules from two different fibres. In CACTUS, the detection step of capsule intersection is a generalization of the cylinder-to-cylinder collision detection (Van Verth and Bishop, 2015). We define the overlapping cost function between two capsules by computing the overlapping of the closest spheres centred in the capsules, as Figure 2C shows. Formally, the closest points between two given capsules [p0, p1, rp0, rp1], [q0, q1, rq0, rq1], with p0, p1, q0, q1 ∈ ℝ3 and , rpi, rqj ∈ R, are the points centred in the capsule ((1 − tp)p0 + tpp1) and, ((1 − tq)q0 + tqq1), where tp, tq are found by the following minimization problem:

g(ta,tb;p0,p1,q0,q1) =  ∥(1−ta)p0+tap1−(1−tb)q0−tbq1∥2 tp,tq: =  arg minta,tb g(ta,tb;  p0,p1,q0,q1),               s.t. 0≤ta,tb≤1

which has a closed-form solution.

After finding the values tp, tq that define the closest centre points between two capsules of different fibres, their overlapping cost function is defined as:

f1(p0,p1,q0,q1,rp0,rp1,rq0,rq1;tp,tq)={D2∥p0−p1​∥∥q0−q1∥rprq,  if D≥00,                                                               if D<0

where,

D: =(1−∥c1−c2∥rp+rq), c1: =(1−tp)p0+tpp1, c2: =(1−tq)q0+tqq1, rp : = (1−tp) rp0+ tprp1, rq : = (1−tq) rq0+ tqrq1,

and tp, tr are the minimal values from the function in Equation 2.

Consequently, the total overlapping cost function in a substrate is computed by adding the evaluated cost of all possible pairwise capsule combinations. If capsules overlap, a penalization is added; otherwise, it is set to zero.

 2.2.3. Implementation details

At last, we mention the technical implementation details of the joint-fibre optimization algorithm, including strategies for reducing computational complexity and the use of specific data structures. Firstly, in the total overlapping cost function, the capsule-to-capsule comparison is a O(n2) problem. To improve computational time, we implemented a fixed-radius-cell data-structure (Turau, 1991) for nearest neighbours queries, reducing the problem to O(n). Since all the cost functions are analytical, we calculated their analytical derivatives for the gradient descent algorithm (see Supplementary material, section joint fibre optimization). We used the adaptative gradient Adagrad (Duchi et al., 2011), iterating until there were no overlapping fibres. All the cost functions, queries, and gradients calculations were implemented in C++ (Stroustrup, 1999) and parallelized with OpenMP (Chandra et al., 2001). To handle bundle crossings, we trim the optimized fibre trajectories to keep only a subregion with fibres that truly belong to the crossing, as shown in Figure 1. This step eliminates boundary fibres that may not fully represent the crossing characteristics.

 2.3. Fibre radial growth 2.3.1. FRG description

After completing the substrate initialization and joint fibre optimization steps, it follows to compute the fibre mesh. Previous studies have managed to achieve a fibre density up to 75% (Altendorf and Jeulin, 2011; Mingasson et al., 2017) with cylindrical-shaped fibres and gamma-distributed diameter, and up 75% with non-cylindrical shaped fibres (Callaghan et al., 2020). In this study, we propose a new method, called Fibre Radial Growth (FRG), to obtain higher packing density and complex axon morphologies beyond the cylindrical shape. The FRG algorithm discretises the 3D space that the fibres occupy to define individual masks for each fibre in it. The FRG algorithm begins to generate the fibre masks by randomly placing seed points within all capsule fibres. These seed points grow iteratively by adding neighbouring points to the fibre mask, employing a breadth-first-search approach through the grid. The seeds grow for a fixed number of iterations as long as they do not interfere with other fibres' boundaries. The propagation through random initializations avoids uniform growth and adds irregularities to the fibre shape, allowing tortuous surface reconstructions in the fibre surfaces. Since the seeding is done inside capsules, the final axon radius in the mesh is related to the radii used in the capsules. We employ two distinct seeding strategies to manage the radius variation effect in our study. The first one depends on the strategy of seeding in the FRG, which depends on how we seed points within the capsule and then grow the seeding points. We can achieve radii variations in intervals like [−ri/2, +ri/2] or (−∈, +∈) depending on whether we decide to seed more randomly or uniformly within the capsule.

In the second case, when we aim to increase the radii variation further in the (−ri/2, ri/2) range, we modify the fibre initialization step. We can define specific patterns in the radii of the fibres' capsules. For example, to incorporate a radii periodicity oscillation, we can set the radii at the start of the capsule to be 1μm and increase the end radii to 2μm. Then, in the subsequent capsule, we can choose to maintain the radii at 2μm or revert them to 1μm, based on the desired frequency of change specified by the user.

Once the FRG step is completed, the fibre density of the particular configuration inputted is maximized. We generate the fibre's outer surface mesh using the marching cubes algorithm (Lewiner et al., 2003; Pedregosa et al., 2011) applied to the fibre mask. This algorithm produces a mesh object consisting of vertices and triangles. Then, we applied a Laplacian smoothing (Herrmann, 1976; Sorkine et al., 2004; Sullivan and Kaszynski, 2019) to remove sharp angles, and finally decimate the mesh to reduce the number of triangles without affecting the morphology of the substrates (Shekhar et al., 1996). Subsequently, we generate a new mesh representing the fibre's inner surface by eroding the previously estimated outer grid and following the same procedure for the meshing. The space between these two surfaces defines the myelin volume.

 2.3.2. Implementation details

Finally, we would like to elaborate on the technical implementation details of the FRG algorithm to mention the specific design choices we made to ensure its computational efficiency. FRG is implemented in Python (Van Rossum and Drake, 2009), parallelized with its multiprocessing ibraries (McKerns et al., 2012), and compiled with Numba (Lam et al., 2015). Image 3D processing and meshing are done using van der Walt et al. (2014), Sullivan and Kaszynski (2019), and Hess (2010). Moreover, the FRG is designed to run a ball-tree structure (Moore et al., 2003) from the Sklearn library (Pedregosa et al., 2011) as a preprocessing to store fibres and their interactions. The fine-tuned FRG algorithm's design allows for the independent execution of fibre growth and meshing on multiple computers in a distributed manner, eliminating the need for multi-thread or computer synchronization.

 3. Experiments

To evaluate the performance of CACTUS, we designed a comprehensive set of substrates with specific geometries. Each experiment below involves several metrics essential for quantifying the microstructure properties of the brain white matter. The metrics include the axon volume fractions, the radius distribution per substrate, the radii change along the fibres, the myelin volume, the g-ratio, the orientation dispersion and bundle crossings. Finally, we conducted testing on the generated substrates and performed Monte-Carlo diffusion simulations to assess their usability and explore the signal decay characteristics associated with these substrates.

3.1. Maximum fibre volume fraction

In our first experiment, we aim to explore the macro-structural parameters of substrates, such as substrate size (i.e., the voxel size in MRI experiments), fibre dispersion, two bundle crossings, and the ability to create high-density packing substrates. We assess the maximum fibre volume fraction that CACTUS achieves in two scenarios: a single bundle and two bundles. In the single bundle case, we generated six substrates with mean angle dispersions of 0, 5, 10, 15, 20, and 25°, respectively. In the two bundles case, we generated five crossing substrates with inter-bundle angles of 30, 45, 60, 75, and 90°, and the fibres of each bundle were initialized with a mean angle dispersion of 5° around the main bundle orientation.

 3.2. Substrates targeting predefined microstructure features

The following two paragraphs describe experiments conducted to explore the ability of CACTUS to replicate desired microstructural parameters into its synthetic substrates. These parameters include the axon volume fraction (AVF), myelin volume fraction (MFV), g-ratio, and radii distribution. We compare the reference values taken from previous histological studies and those achieved by CACTUS.

In the second experiment, we created a series of synthetic substrates that emulate the histological values reported by Stikov et al. (2015) in various white matter regions. Specifically, the target characteristics are the fibre volume fraction, myelin volume fraction, and aggregated g-ratio, g=1-MVF/FVF (Stikov et al., 2015). In our scenario, the axon volume fraction (AVF) is the volume of the fibre inner surface. The myelin volume fraction (MVF) represents the volume of the space between the inner and outer fibre surfaces. The fibre volume fraction (FVF) is the sum of AVF and MVF.

In the third experiment, we investigated the effect of substrate size on radii distribution. To measure the radii distribution for each fibre, we cut the mesh skeleton in an orthogonal plane at regular 1μm intervals and calculated the cross-sectional area of the polygon defined by the plane. The equivalent fibre radius is defined as the radius of a circle with the same area as the polygon (Lee et al., 2019). The global radii distribution per substrate was computed using the mean radius for each fibre.

 3.3. CACTUS substrates usage for Monte-Carlo simulations

The final experiment aims to evaluate the usability of CACTUS substrates in Monte-Carlo diffusion simulations. Synthetic DW-MRI data is generated using meshes obtained from previous experiments. This experiment aims to assess the feasibility and reliability of utilizing CACTUS meshes in Monte-Carlo simulations and examine the resulting DW-MRI signals of such substrates.

To simulate diffusion within non-permeable tissue, we utilized the MCDC Simulator (Rafael-Patino et al., 2020). In this context, the diffusion process within distinct biological structures was assumed to contribute independently to the DW-MRI signal. As a result, the intracellular and extracellular signals were generated separately and combined to generate the overall signal.

The four substrates simulated are composed of ~8,500 fibres. The fibres' outer diameter ranges from 0.5 to 4 um, sampled from a gamma distribution with parameters θ = 1.1, κ = 0.5. Each fibre's inner diameter is calculated from the following log-curve found in Lee et al. (2019).

The simulation substrates have a 300μm3 volume, split in an image size of (42 × 42 × 42) voxels of 7.14μm resolution. Within each voxel, the signal is simulating using random particle sampling with a density of one particle per cubic micrometre Rafael-Patino et al. (2020, 2021) and Romascano et al. (2018) showed that it is a sufficient number of particles to obtain a robust estimation of the diffusion signal in complex fibre geometries.

Particles initiated within the inner diameter of the fibres and outside the outer diameter of the fibres were used to generate the DW-MRI signal. The particles initiated between the outer and inner diameters were discarded because, in this case, we are not simulating the diffusion in the myelin compartment. This was done as previously in the DiSCo Challenge (Rafael-Patino et al., 2021), where no diffusion contrast is assumed in the myelin compartment, however T2 effects could be considered if necessary. The diffusion coefficient, which is user-defined, was fixed to D=0.6×10-3mm2s (corresponding to an ex-vivo diffusivity), for intra and extracellular compartments.

The DW-MRI protocol is based on the protocol of HCP (Fan et al., 2016). It contains four shells of 50 directions with b-values of (1,000, 2,000, 3,000, and 4,000) smm2, and five directions (to calculate parallel and radial decay) with 20 b-values uniformly distributed from (500 to 10,000 ) smm2. The protocol values are fixed for TE = 0.057 s, Δ = 21.8 ms, and δ = 12.9 ms.

 4. Results 4.1. Maximum fibre volume fraction

Figure 3 shows the internal morphology of four substrates consisting of a single bundle with a dispersion of 0, 5, 10, and 20°, respectively. All the substrates were generated with dimensions of 500μm3. Table 1 (top panel) reports the substrate characteristics, including the number of fibres, the obtained fibre volume fraction, and the dispersion parameters. We note that the maximum fibre volume fraction decreased from 94.7 to 90.8% as the dispersion increased from 0 to 25°.

(A–D) Mesh's renders of cross-sections with dimensions of (100μm)2 to visualize the internal morphology of four substrates with a single bundle, each with different mean angular dispersion. For all cases, the outer surface volume is colored black. The inner surface volume is superimposed over the outer volume and colored gray. White represents the extracellular space, i.e., the volume not occupied by any fibre. All bundles are vertically aligned, and the substrates were built to have a mean angular dispersion of (A) 0°, (b) 5°, (C) 10°, and (D) 20°, respectively.

Substrate characteristics, including the number of bundles, mean dispersion angle, mean inter-bundle crossing angle, number of fibres per substrate, and fibre volume fraction (in per cent), respectively.

Nbr of Bundles Bundle dispersion (η) Crossing angle (θ) Nbr of fibres FVF
1 0° - 31,954 94.7%
1 5° - 31,023 93.4%
1 10° - 30,241 92.6%
1 15° - 30,412 92.2%
1 20° - 31,161 91.8%
1 25° - 31,863 90.8%
2 5° 30° 30,026 93.9%
2 5° 45° 30,712 93.3%
2 5° 60° 31,023 93.5%
2 5° 75° 31,152 92.3%
2 5° 90° 30,245 92.2%

The top and bottom panels correspond to the substrates with a single bundle and two bundles. Each row represents a different substrate.

Results from the experiment generating bundle crossings with different inter-bundle angles are depicted in Figure 4 and Table 1 (bottom panel). Figure 4 displays a cross-section of the substrates, where each bundle has a distinctive color for visualization purposes. Although local perturbations in fibre trajectories (on the order of 5°) may occur in the substrates due to the high fibre packing, the average bundle orientation is sustained. The bottom panel of Table 1 reports the fibre volume fraction of these bundle crossing substrates with inter-bundle angles of 30, 45, 60, 75, and 90°. For all the evaluated substrates, the fibre volume fraction remains nearly constant at ~93% (92.2 − 93.9%).

(A–E) Mesh renders of cross-sections with dimensions of (100μm)2 portraying the internal morphology of five substrates consisting of two bundles with different inter-bundle angles. In all cases, the outer surface volume of bundle 1 and bundle 2 is displayed in black. The inner surface volume is superimposed over the outer volume and colored light gray for bundle 1 and dark gray for bundle 2. The extracellular space, representing the volume not occupied by any fibre, is colored white. The dark gray fibres are aligned parallel to the X axis, and the light gray fibres crossed at angles of (A) 30°, (B) 45°, (C) 60°, (D) 75°, and (E) 90°. The outer volume, which is defined by the outer surfaces minus the inner volume, is colored in black for both bundles. Extra axonal space is colored in white.

 4.2. Substrates targeting predefined microstructure features 4.2.1. Axon volume fraction, myelin volume fraction, and g-ratio

We simulated various substrates of a single bundle to mimic microstructure properties previously reported in Stikov et al. (2015). The histological values used as a reference are the myelin volume fraction (MVF), fibre volume fraction (FVF), axonal volume fraction (AVF=FVF-MVF), and g-ratio. The values achieved by CACTUS are shown in Table 2. The difference between the target and obtained substrate properties was lower than 2% in all cases. Examples of the generated substrates and histology data are shown in Figure 5. Electron microscopy images were generously provided by Prof. Nikola Stikov and Dr. Jennifer Campbell, and are used to highlight the geometric similarities of synthetic fibre shapes. On average, for these substrates of 300μm3, the meshes had around 46 million vertices and 92 million faces, and the file size is of 5.1 Gigabytes.

Target microstructure histological properties (left) reported in Stikov et al. (2015), and corresponding properties of the substrates generated by CACTUS (right).

Target Achieved
Substrate AVF MVF FVF g-ratio AVF MVF FVF g-ratio
(a) 25 35 60 64.5 26.0 36.0 62 64.7
(b) 25 43 68 60.6 26.3 43.6 69.9 61.3
(c) 31 44 75 64.2 32.2 43.8 76.07 65
(d) 39 37 76 71.6 41.2 35.0 76.0 73.5

The axon volume fraction (AVF) is the volume of the inner axon surface. The myelin volume fraction (MVF) represents the volume of the space between the inner and outer axon surfaces. The fibre volume fraction (FVF) is the sum of AVF and MVF. The aggregated g-ratio, g=1-MVF/FVF is equal to the mean inner and outer axon radius ratio for all the fibres in the substrate.

Cross-sections of the synthetic substrates constructed to match the statistics of the histological values reported in Table 2. The substrate dimension is 300μm3. For visualization purposes, the axonal space is colored blue, and the myelin is red, and extra axonal space is colored white. (a–d) Correspond to the same substrates shown in Table 2. The bottom panel shows the representative histological images courtesy of Prof. Nikola Stikov and Dr. Jennifer Campbell. The EM images are used to show the similarities in fibre shape and packing.

 4.2.2. Radii distribution and substrate size

Figure 6 show the CACTUS substrates with different sizes, ranging from 30 to 500μm3, and the target and empirical radius distributions obtained for each substrate. The empirical radius distributions closely replicated the targeted ones for substrates equal to or bigger than 200μm3. The optimization algorithm step ran for ~4 h for the largest substrate (right panel) on a node with 64 cores (2.4 GHz) and 400 Mb of RAM. The reconstruction time of the FRG algorithm was ~1 min per fibre, using one core with 500 Mbs of memory per core.

Four 3D substrates of varying sizes: (A) 503, (B) 1003, (C) 2003, and (D) 5003 (μm)3, with 341, 1,316, 4,859, and 33,478 fibres, respectively. The empirical and target radii distributions are displayed on the bottom of each substrate. The empirical distribution better approximates the target distribution as substrate size increases.

We extracted three representative fibre segments from the substrates shown in Figure 5 and displayed them in Figure 7. The top panel of the figure exhibits the cross-sections of the outer and inner surfaces of the fibre, along with the cross-sections of their diameters. The bottom panel shows the diameter distribution of each axon. We observed that, regardless of the tortuosity of the fibre trajectory, the diameter distribution of both the inner and outer diameters of all three cases was centred around the target diameter.

Representative group of three fibres extracted from the substrate shown in Figure 5. In the top panel, we display the fibres with varying diameters and tortuous trajectories. The straightest fibre is presented on the left, while the most tortuous one is displayed on the right. We display each fibre's outer surface in red and its inner surface in blue. We also show the fibre's skeleton in black and cross-sections orthogonal to the fibre's skeleton. The diameter is measured every 1μm along its trajectory. The cross-section cut of the outer surface is shown in red, and the cross-section cut of the inner surface is shown in blue. The bottom panel presents the violin plots of the outer and inner diameters measured. The red (blue) dotted line represents the target outer (inner) diameter of the three fibres.

 4.2.3. Monte-Carlo diffusion simulations

The generated signals for the four substrates from Table 2 depicted in Figure 8. In terms of the parallel diffusivity signal decay, Figures 8A, B demonstrate that the four substrates exhibit similar behavior, making it challenging to distinguish them even at b-values around 8,000 smm2. However, when examining the radial diffusivity signal decay in Figures 8C, D, distinct curves are observed for each substrate. Notably, the logarithmic plot reveals a characteristic tail, indicating the non-Gaussian nature of the diffusion process occurring in the plane perpendicular to the fibre orientations.

Synthetic DW-MRI signals of substrates generated in Table 2 and shown in Figure 8. The X-axis represents the b-value used for measurement, while the Y-axis represents the normalized signal. (A) Parallel signal decay measured along the direction parallel to the substrate fibers (Z-axis). (B) Parallel signal decay plot with a logarithmic scale on the Y-axis. (C) Radial signal decay averaged over four different diffusion directions orthogonal to the orientation of the fibers. (D) Radial signal decay plot with a logarithmic scale on the Y-axis.

 5. Discussion

Over the last 20 years, Monte-Carlo diffusion simulations have been used to optimise DW-MRI data acquisition protocols and validate microstructure models. Nevertheless, doubts have been raised regarding the accuracy of the simple geometries used to construct the diffusion substrates.

Various tools have been developed to address the challenges associated with substrate complexity, such as MEDUSA (Ginsburger et al., 2019) and CONFIG (Callaghan et al., 2020), each offering distinct approaches to substrate generation. MEDUSA primarily focuses on creating substrates with multiple compartments, including axons, oligodendrocytes, and astrocytes, and also models their interactions within the substrate. The axon compartment encompasses features such as axon diameter distribution, bundle dispersion, local tortuosity, myelin presence, Ranvier nodes, and beading. The oligodendrocytes and astrocytes compartments incorporate parameters such as total diameter distribution, body diameter distribution, number of branchings (processes), and a balancing factor. On the other hand, CONFIG employs biologically motivated rules to model the intricate interactions among axons during growth. Parameters within CONFIG include axon mean radius, standard deviation of radius, bundle dispersion, and packing density. Additionally, it encompasses parameters related to axon growth, such as chemoattraction, fiber collapse, cell-adhesion, and fasciculation.

In this work, we introduced CACTUS, a novel framework to produce numerical substrates mimicking white matter tissue with high volume packings, rich microstructural features and geometries that closely matching the desired input parameters. Among the controllable parameters in CACTUS we include the target distribution for the fibre radii, radii variation per fibre, a myelin compartment, target g-ratio, bundle dispersion, bundle crossings, fibre tortuosity, and packing density. The high versatility of CACTUS is founded on its efficient computational implementation and its mathematical formulation divided into three algorithmic steps (substrate initialization, joint fibre optimization, and fibre radial growth) composed of various competing terms controlling different substrate parameters.

To generate the substrates with CACTUS, we introduced a new algorithm to initialise fibre bundles with a target mean degree of orientation dispersion. Moreover, we introduced a novel capsule-based parametrization for optimizing fibre structures. Compared to circle parametrizations (Close et al., 2009; Ginsburger et al., 2019), the capsule parameterization requires fewer parameters, reducing the complexity of the optimization problem. We adapted the cost functions inspired by Close et al. (2009) and Ginsburger et al. (2019) for capsules and provided analytical derivatives, making the optimization faster and computationally more efficient. Finally, we proposed the fibre radial growth algorithm, which increases the fibre packing density in white matter substrates.

CACTUS was able to enhance the complexity of the fibre microstructure. In particular, our results showed CACTUS can produce substrate with fibre volume fraction beyond the 75% previously achieved. CACTUS reached high fibre volume fractions, up to 95% in its substrates (Table 1). Moreover, it consistently reached fibre volume fractions superior to 90% at all the various levels of bundle dispersion and crossing angles (Table 1, Figure 4).

In the single bundle case, the fibre volume fraction was the highest at 94.7% when fibres were aligned and decreased to 90.8% with increasing mean angular dispersion. Conversely, the fibre volume fraction remained consistently around ~93% in the two-bundle cases, regardless of the crossing angle. However, we note that the packing complexity of substrates with a single bundle and two bundles crossing differs. The former mimics the spatial arrangement of thousands of fibres with different crossing angles, which may produce more empty pockets between fibres and less densely packed substrates.

Another important feature of CACTUS is that it can create substrates with statistical characteristics informed by histological data. Indeed, we can closely adhere to the target statistics of axon volume fraction, myelin volume fraction, and g-ratio reported in histological studies (Stikov et al., 2015; see Figure 5). In all cases, the difference between the target and obtained substrate properties was lower than 2% (see Table 1). Notably, CACTUS is the first tool incorporating the g-ratio as a target characteristic and successfully matching it for large-scale substrates.

Also, CACTUS has the capability to generate substrates with a targeted radii distribution. In our experiments, the approximation of the target distribution improves as substrate size increases, as illustrated in Figure 6, underscoring the importance of generating large substrates. Furthermore, we have the availability to measure fibre geometry accurately. For instance, as seen in Figure 7, the generated fibres have a non-constant longitudinal radius and non-circular cross-sections. Despite the tortuous trajectories of the fibres, the diameter distribution remains centred around the target mean outer (inner) diameter of 1.5μm (1.1μm). Additionally, the diameter distribution presented replicates the diameter variations observed in 3D synchrotron images (Andersson et al., 2020), including longitudinal changes and a lack of skewness.

Finally, while previous works were able to achieve substrate sizes between 30 and 100μm3, CACTUS demonstrated a substantial improvement in the generation of larger substrates (Ginsburger et al., 2019; Callaghan et al., 2020). As shown in Figure 6, CACTUS generated substrate sizes ranging from 50 to 500μm3, all with up to a 95% fibre volume fraction. Our tool's ability to generate larger substrate sizes is advantageous for Monte-Carlo diffusion simulations in DW-MRI as it has been shown in previous studies (Rafael-Patino et al., 2020), that substrate sizes larger than 200μm3 can reduce the sampling bias caused by smaller substrate sizes, potentially leading to more accurate DW-MRI numerical simulations (Romascano et al., 2018; Rafael-Patino et al., 2020). In addition, the ability to generate large substrate sizes is advantageous as DW-MRI modelling is moving toward incorporating more microstructure features such as somas, astroglia, and vascularity (Dyer et al., 2017; Lin et al., 2018; Schneider-Mizell et al., 2021). This makes the generation of large substrates essential for capturing these additional features and moving toward more accurate and comprehensive microstructure imaging.

5.1. Limitations and future work

Although CACTUS incorporates complex microstructural features required to mimic some of the most relevant white matter geometrical properties, it still requires fibre-modelling assumptions to reduce the computational burden. Also, CACTUS generates substrates with characteristics resembling those from healthy white matter, but generating pathological tissue requires additional work, which we reserve for future studies.

Additionally, CACTUS focuses solely on generating white matter fibre structures. However, its capacity to generate large substrate sizes expands the potential for including other tissue components in future studies, such as astrocytes, oligodendrocytes, microglia, and capillaries.

Finally, although CACTUS output substrates are suitable for simulators like the MCDC (Rafael-Patino et al., 2020), a thorough analysis is necessary to comprehend the influence of mesh quality, like the number of triangles, on the DW-MRI signals generated by Monte-Carlo simulation. Such analysis is crucial for developing computationally viable simulations.

 5.2. Applications beyond diffusion MR

The applications of CACTUS are not limited to studying white matter microstructure using DW-MRI. For instance, it can be applied in DW-MRI studies outside the brain (Adelnia et al., 2019), where muscle fibres are organized into fascicles. The microscopic arrangement of muscle fibres can vary between different muscle groups, regions of the same muscle, and multiple pathological conditions (Berry et al., 2018). Moreover, the fibre meshes generated by CACTUS could be used in other applications, like Polarized Light Imaging (PLI; Menzel et al., 2015; Amunts and Axer, 2019, a technique used to infer the local fibre orientation in histological brain sections based on the birefringent properties of the myelin sheaths. The limitations of the birefringence PLI model were investigated in Menzel et al. (2015) by generating synthetic PLI data from a hexagonal bundle of straight parallel cylindrical fibres. Although a more general fibre constructor was recently proposed for validating 3D-PLI techniques (Amunts and Axer, 2019), the white matter substrates generated in our study could provide more realistic geometries for conducting similar studies.

 6. Conclusion

The generation of realistic substrates is critical for validating DW-MRI models, as it allows researchers to simulate and analyse the effect of microstructural changes on the DW-MRI signal.

In this work, we introduced CACTUS, a novel framework for generating axonal-like substrates with predefined geometrical features of interest. Our experiments show that CACTUS can generate white matter substrates with the desired spatial dimensions, fibre radii, g-ratio, non-circular cross-sections, tortuous trajectories, smooth surfaces, predefined inter-fibre angles and fibre dispersion. Notably, the generated fibre substrates reached up to 95% fibre volume fraction, the highest density reported in the literature to date, in agreement with previous histology studies. We also generated the large substrates/voxels of up to 500μm3, with dimensions similar to or higher than those used in preclinical MRI scanners, reducing the gap between numerical and real voxel sizes.

In conclusion, the CACTUS substrate generator tool presented in this study has the potential to advance white matter microstructure modelling. It provides a versatile and customisable platform for generating fibre substrates with quantifiable geometrical characteristics. It is open-source and accessible to the broader research community at: http://cactus.epfl.ch, facilitating the validation and comparison of current and future DW-MRI models.

 Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: http://cactus.epfl.ch.

 Author contributions

JV-H contributed to the methodology, coding, simulations, analysis, writing, and visualization. RG contributed to the discussion about substrate generation, the choice of simulation parameters, and writing and reviewing. EC-R and GG contributed to the discussion, methodology, experimental design, and writing review and editing. EF-G contributed to the discussion, methodology, and writing review and editing. J-PT supervised the project, provided funding, contributed to the writing review, and participated in discussions. JR-P contributed to the methodology, experimental design, discussions about simulations, actively contributed to the analysis of results, provided supervision, and writing review and editing. All authors contributed to the article and approved the submitted version.

 Funding

This work was supported by the Swiss National Science Foundation under grants 205320_175974 and 205320_204097. EC-R was supported by the Swiss National Science Foundation (Ambizione grant: PZ00P2_185814). Open access was funded by the École Polytechnique Fédérale de Lausanne.

We acknowledge access to the facilities and expertise of the CIBM Center for Biomedical Imaging, a Swiss research centre of excellence founded and supported by Lausanne University Hospital (CHUV), University of Lausanne (UNIL), Ecole Polytechnique Federale de Lausanne (EPFL), University of Geneva (UNIGE), and Geneva University Hospitals (HUG). We want to express our gratitude to Prof. Nikola Stikov and Dr. Jennifer Campbell for generously donating the histology images from their previous work Stikov et al. (2015). We thank Thomas Yu for his help proofreading this manuscript.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fninf.2023.1208073/full#supplementary-material

 References Adelnia F. Shardell M. Bergeron C. M. Fishbein K. W. Spencer R. G. Ferrucci L. . (2019). Diffusion-weighted MRI with intravoxel incoherent motion modeling for assessment of muscle perfusion in the thigh during post-exercise hyperemia in younger and older adults. NMR Biomed. 32, e4072. 10.1002/nbm.407230861224 Alexander D. C. Hubbard P. L. Hall M. G. Moore E. A. Ptito M. Parker G. J. . (2010). Orientationally invariant indices of axon diameter and density from diffusion MRI. NeuroImage 52, 1374–1389. 10.1016/j.neuroimage.2010.05.04320580932 Altendorf H. Jeulin D. (2011). Random-walk-based stochastic modeling of three-dimensional fiber systems. Phys. Rev. E 83, e041804. 10.1103/PhysRevE.83.04180421599195 Amunts K. Axer M. (2019). Dense fiber modeling for 3D-polarized light imaging simulations. Fut. Trends HPC Disrupt. Scenario 34, 240. Andersson M. Kjer H. M. Rafael-Patino J. Pacureanu A. Pakkenberg B. Thiran J.-P. . (2020). Axon morphology is modulated by the local environment and impacts the noninvasive investigation of its structure–function relationship. Proc. Natl. Acad. Sci. U. S. A. 117, 33649–33659. 10.1073/pnas.201253311733376224 Assaf Y. Basser P. J. (2005). Composite hindered and restricted model of diffusion (charmed) MR imaging of the human brain. Neuroimage 27, 48–58. 10.1016/j.neuroimage.2005.03.04215979342 Assaf Y. Blumenfeld-Katzir T. Yovel Y. Basser P. J. (2008). Axcaliber: a method for measuring axon diameter distribution from diffusion mri. Magnet. Reson. Med. 59, 1347–1354. 10.1002/mrm.2157718506799 Assaf Y. Freidlin R. Z. Rohde G. K. Basser P. J. (2004). New modeling and experimental framework to characterize hindered and restricted water diffusion in brain white matter. Magnet. Reson. Med. 52, 965–978. 10.1002/mrm.2027415508168 Baxter G. T. Frank L. R. (2013). A computational model for diffusion weighted imaging of myelinated white matter. Neuroimage 75, 204–212. 10.1016/j.neuroimage.2013.02.07623507381 Berry D. B. Regner B. Galinsky V. Ward S. R. Frank L. R. (2018). Relationships between tissue microstructure and the diffusion tensor in simulated skeletal muscle. Magnet. Reson. Med. 80, 317–329. 10.1002/mrm.2699329090480 Bihan D. L. (1995). Molecular diffusion, tissue microdynamics and microstructure. NMR Biomed. 8, 375–386. 10.1002/nbm.19400807118739274 Bolduan F. Grosser S. Vida I. (2020). Minimizing shrinkage of acute brain slices using metal spacers during histological embedding. Brain Struct. Funct. 225, 2577–2589. 10.1007/s00429-020-02141-332918613 Brückner G. Härtig W. Kacza J. Seeger J. Welt K. Brauer K. (1996). Extracellular matrix organization in various regions of rat brain grey matter. J. Neurocytol. 25, 333–346. 10.1007/BF022848068818977 Budde M. D. Frank J. A. (2010). Neurite beading is sufficient to decrease the apparent diffusion coefficient after ischemic stroke. Proc. Natl. Acad. Sci. U. S. A. 107, 14472–14477. 10.1073/pnas.100484110720660718 Callaghan R. Alexander D. C. Palombo M. Zhang H. (2020). ConFiG: Contextual Fibre Growth to generate realistic axonal packing for diffusion MRI simulation. NeuroImage 220, 117107. 10.1016/j.neuroimage.2020.11710732622984 Caminiti R. Carducci F. Piervincenzi C. Battaglia-Mayer A. Confalone G. Visco-Comandini F. . (2013). Diameter, length, speed, and conduction delay of callosal axons in macaque monkeys and humans: comparing data from histology and magnetic resonance imaging diffusion tractography. J. Neurosci. 33, 14501–14511. 10.1523/JNEUROSCI.0761-13.201324005301 Campbell J. S. W. Siddiqi K. Rymar V. V. Sadikot A. F. Pike G. B. (2005). Flow-based fiber tracking with diffusion tensor and Q-ball data: validation and comparison to principal diffusion direction techniques. NeuroImage 27, 725–736. 10.1016/j.neuroimage.2005.05.01416111897 Canales-Rodr-guez E. J. Legarreta J. H. Pizzolato M. Rensonnet G. Girard G. Rafael-Patino J. . (2019). Sparse wars: a survey and comparative study of spherical deconvolution algorithms for diffusion mri. NeuroImage 184, 140–160. 10.1016/j.neuroimage.2018.08.07130193974 Canales-Rodr-guez E. J. Pizzolato M. Piredda G. F. Hilbert T. Kunz N. Pot C. . (2021). Comparison of non-parametric T2 relaxometry methods for myelin water quantification. Med. Image Anal. 69, 101959. 10.1016/j.media.2021.10195933581618 Chandra R. Dagum L. Kohr D. Menon R. Maydan D. McDonald J. (2001). Parallel Programming in OpenMP. Burlington, MA: Morgan kaufmann. Close T. G. Tournier J. D. Calamante F. Johnston L. A. Mareels I. Connelly A. (2009). A software tool to generate simulated white matter structures for the assessment of fibre-tracking algorithms. NeuroImage 47, 1288–1300. 10.1016/j.neuroimage.2009.03.07719361565 Côté M.-A. Girard G. Bor A. Garyfallidis E. Houde J.-C. Descoteaux M. (2013). Tractometer: towards validation of tractography pipelines. Med. Image Anal. 17, 844–857. 10.1016/j.media.2013.03.00923706753 Daducci A. Canales-Rodríguez E. J. Zhang H. Dyrby T. B. Alexander D. C. Thiran J.-P. (2015). Accelerated microstructure imaging via convex optimization (amico) from diffusion MRI data. NeuroImage 105, 32–44. 10.1016/j.neuroimage.2014.10.02625462697 Dam A. M. (1979). Shrinkage of the brain during histological procedures with fixation in formaldehyde solutions of different concentrations. J. Hirnforschung 20, 115–119.556570 Drobnjak I. Zhang H. Ianuş A. Kaden E. Alexander D. C. (2016). Pgse, ogse, and sensitivity to axon diameter in diffusion MRI: insight from a simulation study. Magnet. Reson. Med. 75, 688–700. 10.1002/mrm.2563125809657 Duchi J. Hazan E. Singer Y. (2011). Adaptive subgradient methods for online learning and stochastic optimization. J. Machine Learn. Res. 12. Dyer E. L. Roncal W. G. Prasad J. A. Fernandes H. L. Gürsoy D. De Andrade V. . (2017). Quantifying mesoscale neuroanatomy using x-ray microtomography. Eneuro 4, 2017. 10.1523/ENEURO.0195-17.201729085899 Dyrby T. B. Baar W. F. Alexander D. C. Jelsing J. Garde E. Sgaard L. V. (2011). An ex vivo imaging pipeline for producing high-quality and high-resolution diffusion-weighted imaging datasets. Hum. Brain Map. 32, 544–563. 10.1002/hbm.2104320945352 Fan Q. Witzel T. Nummenmaa A. Van Dijk K. R. Van Horn J. D. Drews M. K. . (2016). MGH–USC human connectome project datasets with ultra-high b-value diffusion MRI. Neuroimage 124, 1108–1114. 10.1016/j.neuroimage.2015.08.07526364861 Fieremans E. De Deene Y. Delputte S. Özdemir M. S. D'Asseler Y. Vlassenbroeck J. . (2008). Simulation and experimental verification of the diffusion in an anisotropic fiber phantom. J. Magnet. Reson. 190, 189–199. 10.1016/j.jmr.2007.10.01418023218 Fillard P. Descoteaux M. Goh A. Gouttard S. Jeurissen B. Malcolm J. . (2011). Quantitative evaluation of 10 tractography algorithms on a realistic diffusion MR phantom. NeuroImage 56, 220–234. 10.1016/j.neuroimage.2011.01.03221256221 Ginsburger K. Matuschke F. Poupon F. Mangin J. F. Axer M. Poupon C. (2019). MEDUSA: a GPU-based tool to create realistic phantoms of the brain microstructure using tiny spheres. NeuroImage 193, 10–24. 10.1016/j.neuroimage.2019.02.05530849528 Grussu F. Ianuş A. Tur C. Prados F. Schneider T. Kaden E. . (2019). Relevance of time-dependence for clinically viable diffusion imaging of the spinal cord. Magnet. Reson. Med. 81, 1247–1264. 10.1002/mrm.2746330229564 Hall M. G. Alexander D. C. (2009). Convergence and parameter choice for monte-carlo simulations of diffusion MRI. IEEE Trans. Med. Imag. 28, 1354–1364. 10.1109/TMI.2009.201575619273001 Harkins K. D. Beaulieu C. Xu J. Gore J. C. Does M. D. (2021). A simple estimate of axon size with diffusion MRI. Neuroimage 227, 117619. 10.1016/j.neuroimage.2020.11761933301942 Herrmann L. R. (1976). Laplacian-isoparametric grid generation scheme. J. Eng. Mech. Div. 102, 749–756. 10.1061/JMCEA3.0002158 Hess R. (2010). Blender Foundations: The Essential Guide to Learning Blender 2.6. Waltham, MA: Focal Press. Innocenti G. M. Caminiti R. (2017). Axon diameter relates to synaptic bouton size: structural properties define computationally different types of cortical connections in primates. Brain Struct. Funct. 222, 1169–1177. 10.1007/s00429-016-1266-127372337 Jelescu I. O. Budde M. D. (2017). Design and validation of diffusion mri models of white matter. Front. Phys. 28, 61. 10.3389/fphy.2017.0006129755979 Kakkar L. S. Bennett O. F. Siow B. Richardson S. Ianuş A. Quick T. . (2018). Low frequency oscillating gradient spin-echo sequences improve sensitivity to axon diameter: an experimental study in viable nerve tissue. Neuroimage 182, 314–328. 10.1016/j.neuroimage.2017.07.06028774648 Korogod N. Petersen C. C. Knott G. W. (2015). Ultrastructural analysis of adult mouse neocortex comparing aldehyde perfusion with cryo fixation. Elife 4, e05793. 10.7554/eLife.0579326259873 Lam S. K. Pitrou A. Seibert S. (2015). “Numba: A llvm-based python jit compiler,” in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC. 1–6.35500032 Lavdas I. Behan K. C. Papadaki A. McRobbie D. W. Aboagye E. O. (2013). A phantom for diffusion-weighted MRI (DW-MRI). J. Magnet. Reson. Imag. 38, 173–179. 10.1002/jmri.2395023576443 Lee H.-H. Jespersen S. N. Fieremans E. Novikov D. S. (2020). The impact of realistic axonal shape on axon diameter estimation using diffusion MRI. NeuroImage 223, 117228. 10.1016/j.neuroimage.2020.11722832798676 Lee H.-H. Yaros K. Veraart J. Pathan J. L. Liang F.-X. Kim S. G. . (2019). Along-axon diameter variation and axonal orientation dispersion revealed with 3D electron microscopy: implications for quantifying brain white matter microstructure with histology and diffusion MRI. Brain Struct. Funct. 223, 117228.30790073 Lewiner T. Lopes H. Vieira A. W. Tavares G. (2003). Efficient implementation of marching cubes' cases with topological guarantees. J. Graph. Tools 8, 1–15. 10.1080/10867651.2003.10487582 Lin C. Huang Y. Quan T. Zhang Y. (2018). Modelling brain-wide neuronal morphology via rooted cayley trees. Sci. Rep. 8, 1–10. 10.1038/s41598-018-34050-130353025 Mackay A. Whittall K. Adler J. Li D. Paty D. Graeb D. (1994). In vivo visualization of myelin water in brain by magnetic resonance. Magnet. Reson. Med. 31, 673–677. 10.1002/mrm.19103106148057820 Maier-Hein K. H. Neher P. F. Houde J.-C. Côté M.-A. Garyfallidis E. Zhong J. . (2017). The challenge of mapping the human connectome based on diffusion tractography. Nat. Commun. 8, 1349. 10.1038/s41467-017-01285-x31685826 McKerns M. M. Strand L. Sullivan T. Fang A. Aivazis M. A. (2012). Building a framework for predictive science. arXiv preprint arXiv:1202.1056. 10.25080/Majora-ebaa42b7-00d Menzel M. Michielsen K. De Raedt H. Reckfort J. Amunts K. Axer M. (2015). A jones matrix formalism for simulating three-dimensional polarized light imaging of brain tissue. J. Royal Soc. Interf. 12, 20150734. 10.1098/rsif.2015.073426446561 Mingasson T. Duval T. Stikov N. Cohen-Adad J. (2017). Axonpacking: an open-source software to simulate arrangements of axons in white matter. Front. Neuroinformat. 11, 5. 10.3389/fninf.2017.0000528197091 Moore A. Gray A. Liu T. (2003). New algorithms for efficient high dimensional non-parametric classification. Adv. Neural Inform. Process. Syst. 16. Murday J. S. Cotts R. M. (1968). Self-diffusion coefficient of liquid lithium. J. Chem. Phys. 48, 4938–4945. 10.1063/1.1668160 Narvaez-Delgado O. Rojas-Vite G. Coronado-Leija R. Ramírez-Manzanares A. Marroquín J. L. Noguez-Imm R. . (2019). Histological and diffusion-weighted magnetic resonance imaging data from normal and degenerated optic nerve and chiasm of the rat. Data Brief 26, 104399. 10.1016/j.dib.2019.10439931516943 Neher P. F. Laun F. B. Stieltjes B. Maier-Hein K. H. (2014). Fiberfox: facilitating the creation of realistic white matter software phantoms. Magnet. Reson. Med. 72, 1460–1470. 10.1002/mrm.2504524323973 Neuman C. H. (1974). Spin echo of spins diffusing in a bounded medium. J. Chem. Phys. 60, 4508–4511. 10.1063/1.1680931 Nilsson M. Lasič S. Drobnjak I. Topgaard D. Westin C.-F. (2017). Resolution limit of cylinder diameter estimation by diffusion MRI: The impact of gradient waveform and orientation dispersion. NMR Biomed. 30, e3711. 10.1002/nbm.371128318071 Nilsson M. Lätt J. Ståhlberg F. van Westen D. Hagslätt H. (2012). The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study. NMR Biomed. 25, 795–805. 10.1002/nbm.179522020832 Novikov D. S. Fieremans E. Jespersen S. N. Kiselev V. G. (2019). Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation. NMR Biomed. 32, e3998. 10.1002/nbm.399830321478 Novikov D. S. Kiselev V. G. Jespersen S. N. (2018). On modeling. Magnet. Reson. Med. 79, 3172–3193. 10.1002/mrm.2710129493816 Parizel P. M. Van Goethem J. Özsarlak Ö. Maes M. Phillips C. (2005). New developments in the neuroradiological diagnosis of craniocerebral trauma. Eur. Radiol. 15, 569–581. 10.1007/s00330-004-2558-z15696294 Pedregosa F. Varoquaux G. Gramfort A. Michel V. Thirion B. Grisel O. . (2011). Scikit-learn: machine learning in Python. J. Machine Learn. Res. 12, 2825–2830. Plante I. Cucinotta F. A. (2013). Monte-carlo simulation of particle diffusion in various geometries and application to chemistry and biology. Theory Appl. Monte Carlo Simul. 2013, 193–225. 10.5772/53203 Rafael-Patino J. Girard G. Truffet R. Pizzolato M. Thiran J. P. Caruyer E. (2021). “The microstructural features of the diffusion-simulated connectivity (disco) dataset,” in Computational Diffusion MRI: 12th International Workshop, CDMRI 2021, Held in Conjunction with MICCAI 2021, Strasbourg, France, October 1, 2021, Proceedings 12 (Berlin: Springer), 159–170.34632021 Rafael-Patino J. Romascano D. Ramirez-Manzanares A. Canales-Rodríguez E. J. Girard G. Thiran J. P. (2020). Robust monte-carlo simulations in diffusion-MRI: effect of the substrate complexity and parameter choice on the reproducibility of results. Front. Neuroinformat. 14, 8. 10.3389/fninf.2020.0000832210781 Romascano D. Rafael-Patino J. Jelescu I. Barakovic M. Tim B. Jean-Philippe D. . (2018). “Voxel size matters: big voxels are required to generate realistic extra-axonal DMRI signals from monte carlo simulations,” in ISMRM. 1–2. Schilling K. Gao Y. Janve V. Stepniewska I. Landman B. A. Anderson A. W. (2017). Can increased spatial resolution solve the crossing fiber problem for diffusion MRI? NMR Biomed. 30, e3787. 10.1002/nbm.378728915311 Schilling K. G. Nath V. Hansen C. Parvathaneni P. Blaber J. Gao Y. . (2019). Limits to anatomical accuracy of diffusion tractography using modern approaches. Neuroimage 185, 1–11. 10.1016/j.neuroimage.2018.10.02930317017 Schneider-Mizell C. M. Bodor A. L. Collman F. Brittain D. Bleckert A. Dorkenwald S. . (2021). Structure and function of axo-axonic inhibition. Elife 10, e73783. 10.7554/eLife.7378334851292 Sepehrband F. Alexander D. C. Clark K. A. Kurniawan N. D. Yang Z. Reutens D. C. (2016). Parametric probability distribution functions for axon diameters of corpus callosum. Front. Neuroanat. 10, 59. 10.3389/fnana.2016.0005927303273 Shekhar R. Fayyad E. Yagel R. Cornhill J. F. (1996). “Octree-based decimation of marching cubes surfaces,” in Proceedings of Seventh Annual IEEE Visualization'96 (San Francisco, CA: IEEE), 335–342. Simpson J. Carr H. (1958). Diffusion and nuclear spin relaxation in water. Phys. Rev. 111, 1201. 10.1103/PhysRev.111.1201 Söderman O. Jönsson B. (1995). Restricted diffusion in cylindrical geometry. J. Magnet. Reson. Ser. A 117, 94–97. 10.1006/jmra.1995.0014 Sorkine O. Cohen-Or D. Lipman Y. Alexa M. Rössl C. Seidel H.-P. (2004). “Laplacian surface editing,” in Proceedings of the 2004 Eurographics/ACM SIGGRAPH Symposium on Geometry Processing (ACM), 175–184. Sporns O. (2011). The non-random brain: efficiency, economy, and complex dynamics. Front. Comput. Neurosci. 5, 5. 10.3389/fncom.2011.0000521369354 Stanisz G. J. Wright G. A. Henkelman R. M. Szafer A. (1997). An analytical model of restricted diffusion in bovine optic nerve. Magnet. Reson. Med. 37, 103–111. 10.1002/mrm.19103701158978638 Stejskal E. O. Tanner J. E. (1965). Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J. Chem. Phys. 42, 288–292. 10.1063/1.169569034864390 Stikov N. Campbell J. S. Stroh T. Lavelée M. Frey S. Novek J. . (2015). In vivo histology of the myelin g-ratio with magnetic resonance imaging. Neuroimage 118, 397–405. 10.1016/j.neuroimage.2015.05.02326004502 Stroustrup B. (1999). “An overview of the C++ programming language,” in Handbook of Object Technology, ed S. Zamir (Boca Raton, FL: CRC Press), 72. Sullivan B. Kaszynski A. (2019). PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK). J. Open Sourc. Softw. 4, 1450. 10.21105/joss.01450 Syková E. Nicholson C. (2008). Diffusion in brain extracellular space. Physiol. Rev. 88, 1277–1340. 10.1152/physrev.00027.200718923183 To X. V. Mohamed A. Z. Cumming P. Nasrallah F. A. (2022). Subacute cytokine changes after a traumatic brain injury predict chronic brain microstructural alterations on advanced diffusion imaging in the male rat. Brain Behav. Immunity 102, 137–150. 10.1016/j.bbi.2022.02.01735183698 Tønnesen J. Inavalli V. K. Nägerl U. V. (2018). Super-resolution imaging of the extracellular space in living brain tissue. Cell 172, 1108–1121. 10.1016/j.cell.2018.02.00729474910 Tournier J.-D. Calamante F. Connelly A. (2007). Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution. NeuroImage 35, 1459–1472. 10.1016/j.neuroimage.2007.02.01617379540 Tournier J. D. Yeh C. H. Calamante F. Cho K. H. Connelly A. Lin C. P. (2008). Resolving crossing fibres using constrained spherical deconvolution: Validation using diffusion-weighted imaging phantom data. NeuroImage 42, 617–625. 10.1016/j.neuroimage.2008.05.00218583153 Truffet R. Rafael-Patino J. Girard G. Pizzolato M. Barillot C. Thiran J.-P. . (2020). “An evolutionary framework for microstructure-sensitive generalized diffusion gradient waveforms,” in International Conference on Medical Image Computing and Computer-Assisted Intervention (Berlin: Springer), 94–103. Tuch D. S. (2004). Q-ball imaging. Magnet. Reson. Med. 52, 1358–1372. 10.1002/mrm.2027915562495 Turau V. (1991). Fixed-radius near neighbors search. Informat. Process. Lett. 39, 201–203. 10.1016/0020-0190(91)90180-P34520382 van der Walt S. Schönberger J. L. Nunez-Iglesias J. Boulogne F. Warner J. D. Yager N. . (2014). Scikit-image: image processing in Python. PeerJ 2, e453. 10.7287/peerj.preprints.336v2 van Gelderen P. de Vleeschouwer M. H. M. DesPres D. Pekar J. van Zijl P. C. M. Moonen C. T. W. (1994). Water diffusion and acute stroke. Magnet. Reson. Med. 31, 154–163. 10.1002/mrm.19103102098133751 Van Rossum G. Drake F. L. (2009). Python 3 Reference Manual. Scotts Valley, CA: CreateSpace. Van Verth J. M. Bishop L. M. (2015). Essential Mathematics for Games and Interactive Applications. Boca Raton, FL: CRC Press. Veraart J. Nunes D. Rudrapatna U. Fieremans E. Jones D. K. Novikov D. S. . (2020). Noninvasive quantification of axon radii using diffusion MRI. Elife 9, e49855. 10.7554/eLife.4985532048987 Veraart J. Raven E. P. Edwards L. J. Weiskopf N. Jones D. K. (2021). The variability of mr axon radii estimates in the human white matter. Hum. Brain Map. 42, 2201–2213. 10.1002/hbm.2535933576105 Warner W. Palombo M. Cruz R. Callaghan R. Shemesh N. Jones D. K. . (2023). Temporal diffusion ratio (TDR) for imaging restricted diffusion: optimisation and pre-clinical demonstration. NeuroImage 269, 119930. 10.1016/j.neuroimage.2023.11993036750150 Zhang H. Schneider T. Wheeler-Kingshott C. A. Alexander D. C. (2012). NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain. Neuroimage 61, 1000–1016. 10.1016/j.neuroimage.2012.03.07222484410 Zhou F.-L. Li Z. Gough J. E. Cristinacce P. L. H. Parker G. J. (2018). Axon mimicking hydrophilic hollow polycaprolactone microfibres for diffusion magnetic resonance imaging. Mater. Design 137, 394–403. 10.1016/j.matdes.2017.10.04729307950