CACTUS: a computational framework for generating realistic white matter microstructure substrates

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μm³ 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.

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μm³, 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μm³. 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).

FIGURE 1

Figure 1. 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 [p₀, p₁, r₀, r₁], where p₀, p₁ ∈ ℝ³ are the initial/ending points of the capsule, and r₀, r₁ ∈ ℝ are their respective radius (see Figure 2A). In our scenario, the length of a capsule (distance form p_i to p_i+1) is not restricted, but we suggest the ranges between $\frac{1}{2}{r}_i$ up to 2r_i, and the change sampling frequency increases as the radii decrease.

FIGURE 2

Figure 2. (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 ${\mathcal{S}}^a$, with m_a control points, is composed of the capsules determined by the subsequent point pairs as $\left[x_i^a,x_{i+1}^a,r_i^a,r_{i+1}^a\right]$ with control points $\left\{x_0{}^a,x_1{}^a,\dots x_{m_a-1}{}^a\right\}\subset {ℝ}^3$ and associated radius $\left\{r_0^a,r_1^a,\dots r_{m_a-1}^a\right\}\subset ℝ.$

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 [p₀, p₁, r_p0, r_p1], [q₀, q₁, r_q0, r_q1], with p₀, p₁, q₀, q₁ ∈ ℝ³ and ,r_pi, r_qj ∈ R, are the points centred in the capsule ((1 − t_p)p₀ + t_pp₁) and, ((1 − t_q)q₀ + t_qq₁), where t_p, t_q are found by the following minimization problem:

$\begin{array}{ll}g(t_a,t_b;p_0,p_1,q_0,q_1)=\parallel (1-t_a)p_0+t_a{p}_1-(1-t_b)q_0-t_b{q}_1{\parallel }^2& (1)\end{array}$

$\begin{array}{l}{t}_p,t_q:=\underset{t_a,t_b}{\text{ arg min}}g(t_a,t_b;p_0,p_1,q_0,q_1),\\ s.t.0\le t_a,t_b\le 1& (2)\end{array}$

which has a closed-form solution.

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

$\begin{array}{l}{f}_1\left(p_0,p_1,q_0,q_1,r_{p_0},r_{p_1},r_{q_0},r_{q_1};t_p,t_q\right)=\\ \{\begin{array}{l}{D}^2\parallel p_0-p_1\text{}\parallel \parallel q_0-q_1\parallel r_p{r}_q,\text{ if }D\ge 0\\ 0,\text{ if }D<0\end{array}& (3)\end{array}$

where,

$\begin{array}{ll}D:=\left(1-\frac{\parallel c_1-c_2\parallel }{r_p+r_q}\right),& (4)\end{array}$

$\begin{array}{ll}{c}_1:=(1-t_p)p_0+t_p{p}_1,& (5)\end{array}$

$\begin{array}{ll}{c}_2:=(1-t_q)q_0+t_q{q}_1,& (6)\end{array}$

$\begin{array}{ll}{r}_p\text{ : = }(1-t_p)r_{p_0}\text{+ }{t}_p{r}_{p1},& (7)\end{array}$

$\begin{array}{ll}{r}_q\text{ : = }(1-t_q)r_{q_0}\text{+ }{t}_q{r}_{q1},& (8)\end{array}$

and t_p, t_r 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 $\mathcal{O}(n^2)$ problem. To improve computational time, we implemented a fixed-radius-cell data-structure (Turau, 1991) for nearest neighbours queries, reducing the problem to $\mathcal{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 [−r_i/2, +r_i/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 (−r_i/2, r_i/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=\sqrt{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μm³ 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\times 10^{-3}\frac{{\mathrm{\text{mm}}}^2}{\mathrm{\text{s}}}$ (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) $\frac{\mathrm{\text{s}}}{{\mathrm{\text{mm}}}^2}$, and five directions (to calculate parallel and radial decay) with 20 b-values uniformly distributed from (500 to 10,000 ) $\frac{\mathrm{\text{s}}}{{\mathrm{\text{mm}}}^2}$. 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μm³. 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°.

FIGURE 3

Figure 3. (A–D) Mesh's renders of cross-sections with dimensions of (100μm)² 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.

TABLE 1

Table 1. 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.

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

FIGURE 4

Figure 4. (A–E) Mesh renders of cross-sections with dimensions of (100μm)² 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μm³, the meshes had around 46 million vertices and 92 million faces, and the file size is of 5.1 Gigabytes.

TABLE 2

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

FIGURE 5

Figure 5. Cross-sections of the synthetic substrates constructed to match the statistics of the histological values reported in Table 2. The substrate dimension is 300μm³. 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μm³, 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μm³. 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.

FIGURE 6

Figure 6. Four 3D substrates of varying sizes: (A) 50³, (B) 100³, (C) 200³, and (D) 500³ (μm)³, 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.

FIGURE 7

Figure 7. 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 $\frac{\mathrm{\text{s}}}{{\mathrm{\text{mm}}}^2}$. 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.

FIGURE 8

Figure 8. 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μm³, 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μm³, 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μm³ 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μm³, 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.

Acknowledgments

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., et al. (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. doi: 10.1002/nbm.4072

PubMed Abstract | CrossRef Full Text | Google Scholar

Alexander, D. C., Hubbard, P. L., Hall, M. G., Moore, E. A., Ptito, M., Parker, G. J., et al. (2010). Orientationally invariant indices of axon diameter and density from diffusion MRI. NeuroImage 52, 1374–1389. doi: 10.1016/j.neuroimage.2010.05.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Altendorf, H., and Jeulin, D. (2011). Random-walk-based stochastic modeling of three-dimensional fiber systems. Phys. Rev. E 83, e041804. doi: 10.1103/PhysRevE.83.041804

PubMed Abstract | CrossRef Full Text | Google Scholar

Amunts, K., and Axer, M. (2019). Dense fiber modeling for 3D-polarized light imaging simulations. Fut. Trends HPC Disrupt. Scenario 34, 240.

Google Scholar

Andersson, M., Kjer, H. M., Rafael-Patino, J., Pacureanu, A., Pakkenberg, B., Thiran, J.-P., et al. (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. doi: 10.1073/pnas.2012533117

PubMed Abstract | CrossRef Full Text | Google Scholar

Assaf, Y., and Basser, P. J. (2005). Composite hindered and restricted model of diffusion (charmed) MR imaging of the human brain. Neuroimage 27, 48–58. doi: 10.1016/j.neuroimage.2005.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Assaf, Y., Blumenfeld-Katzir, T., Yovel, Y., and Basser, P. J. (2008). Axcaliber: a method for measuring axon diameter distribution from diffusion mri. Magnet. Reson. Med. 59, 1347–1354. doi: 10.1002/mrm.21577

PubMed Abstract | CrossRef Full Text | Google Scholar

Assaf, Y., Freidlin, R. Z., Rohde, G. K., and 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. doi: 10.1002/mrm.20274

PubMed Abstract | CrossRef Full Text | Google Scholar

Baxter, G. T., and Frank, L. R. (2013). A computational model for diffusion weighted imaging of myelinated white matter. Neuroimage 75, 204–212. doi: 10.1016/j.neuroimage.2013.02.076

PubMed Abstract | CrossRef Full Text | Google Scholar

Berry, D. B., Regner, B., Galinsky, V., Ward, S. R., and Frank, L. R. (2018). Relationships between tissue microstructure and the diffusion tensor in simulated skeletal muscle. Magnet. Reson. Med. 80, 317–329. doi: 10.1002/mrm.26993

PubMed Abstract | CrossRef Full Text | Google Scholar

Bihan, D. L. (1995). Molecular diffusion, tissue microdynamics and microstructure. NMR Biomed. 8, 375–386. doi: 10.1002/nbm.1940080711

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolduan, F., Grosser, S., and Vida, I. (2020). Minimizing shrinkage of acute brain slices using metal spacers during histological embedding. Brain Struct. Funct. 225, 2577–2589. doi: 10.1007/s00429-020-02141-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Brückner, G., Härtig, W., Kacza, J., Seeger, J., Welt, K., and Brauer, K. (1996). Extracellular matrix organization in various regions of rat brain grey matter. J. Neurocytol. 25, 333–346. doi: 10.1007/BF02284806

PubMed Abstract | CrossRef Full Text | Google Scholar

Budde, M. D., and 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. doi: 10.1073/pnas.1004841107

PubMed Abstract | CrossRef Full Text | Google Scholar

Callaghan, R., Alexander, D. C., Palombo, M., and Zhang, H. (2020). ConFiG: Contextual Fibre Growth to generate realistic axonal packing for diffusion MRI simulation. NeuroImage 220, 117107. doi: 10.1016/j.neuroimage.2020.117107

PubMed Abstract | CrossRef Full Text | Google Scholar

Caminiti, R., Carducci, F., Piervincenzi, C., Battaglia-Mayer, A., Confalone, G., Visco-Comandini, F., et al. (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. doi: 10.1523/JNEUROSCI.0761-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, J. S. W., Siddiqi, K., Rymar, V. V., Sadikot, A. F., and 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. doi: 10.1016/j.neuroimage.2005.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Canales-Rodr-guez, E. J., Legarreta, J. H., Pizzolato, M., Rensonnet, G., Girard, G., Rafael-Patino, J., et al. (2019). Sparse wars: a survey and comparative study of spherical deconvolution algorithms for diffusion mri. NeuroImage 184, 140–160. doi: 10.1016/j.neuroimage.2018.08.071

PubMed Abstract | CrossRef Full Text | Google Scholar

Canales-Rodr-guez, E. J., Pizzolato, M., Piredda, G. F., Hilbert, T., Kunz, N., Pot, C., et al. (2021). Comparison of non-parametric T2 relaxometry methods for myelin water quantification. Med. Image Anal. 69, 101959. doi: 10.1016/j.media.2021.101959

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandra, R., Dagum, L., Kohr, D., Menon, R., Maydan, D., and McDonald, J. (2001). Parallel Programming in OpenMP. Burlington, MA: Morgan kaufmann.

Google Scholar

Close, T. G., Tournier, J. D., Calamante, F., Johnston, L. A., Mareels, I., and Connelly, A. (2009). A software tool to generate simulated white matter structures for the assessment of fibre-tracking algorithms. NeuroImage 47, 1288–1300. doi: 10.1016/j.neuroimage.2009.03.077

PubMed Abstract | CrossRef Full Text | Google Scholar

Côté, M.-A., Girard, G., Bor, A., Garyfallidis, E., Houde, J.-C., and Descoteaux, M. (2013). Tractometer: towards validation of tractography pipelines. Med. Image Anal. 17, 844–857. doi: 10.1016/j.media.2013.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Daducci, A., Canales-Rodríguez, E. J., Zhang, H., Dyrby, T. B., Alexander, D. C., and Thiran, J.-P. (2015). Accelerated microstructure imaging via convex optimization (amico) from diffusion MRI data. NeuroImage 105, 32–44. doi: 10.1016/j.neuroimage.2014.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Dam, A. M. (1979). Shrinkage of the brain during histological procedures with fixation in formaldehyde solutions of different concentrations. J. Hirnforschung 20, 115–119.

PubMed Abstract | Google Scholar

Drobnjak, I., Zhang, H., Ianuş, A., Kaden, E., and 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. doi: 10.1002/mrm.25631

PubMed Abstract | CrossRef Full Text | Google Scholar

Duchi, J., Hazan, E., and Singer, Y. (2011). Adaptive subgradient methods for online learning and stochastic optimization. J. Machine Learn. Res. 12.

Google Scholar

Dyer, E. L., Roncal, W. G., Prasad, J. A., Fernandes, H. L., Gürsoy, D., De Andrade, V., et al. (2017). Quantifying mesoscale neuroanatomy using x-ray microtomography. Eneuro 4, 2017. doi: 10.1523/ENEURO.0195-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Dyrby, T. B., Baar, W. F., Alexander, D. C., Jelsing, J., Garde, E., and 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. doi: 10.1002/hbm.21043

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, Q., Witzel, T., Nummenmaa, A., Van Dijk, K. R., Van Horn, J. D., Drews, M. K., et al. (2016). MGH–USC human connectome project datasets with ultra-high b-value diffusion MRI. Neuroimage 124, 1108–1114. doi: 10.1016/j.neuroimage.2015.08.075

PubMed Abstract | CrossRef Full Text | Google Scholar

Fieremans, E., De Deene, Y., Delputte, S., Özdemir, M. S., D'Asseler, Y., Vlassenbroeck, J., et al. (2008). Simulation and experimental verification of the diffusion in an anisotropic fiber phantom. J. Magnet. Reson. 190, 189–199. doi: 10.1016/j.jmr.2007.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Fillard, P., Descoteaux, M., Goh, A., Gouttard, S., Jeurissen, B., Malcolm, J., et al. (2011). Quantitative evaluation of 10 tractography algorithms on a realistic diffusion MR phantom. NeuroImage 56, 220–234. doi: 10.1016/j.neuroimage.2011.01.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Ginsburger, K., Matuschke, F., Poupon, F., Mangin, J. F., Axer, M., and Poupon, C. (2019). MEDUSA: a GPU-based tool to create realistic phantoms of the brain microstructure using tiny spheres. NeuroImage 193, 10–24. doi: 10.1016/j.neuroimage.2019.02.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Grussu, F., Ianuş, A., Tur, C., Prados, F., Schneider, T., Kaden, E., et al. (2019). Relevance of time-dependence for clinically viable diffusion imaging of the spinal cord. Magnet. Reson. Med. 81, 1247–1264. doi: 10.1002/mrm.27463

PubMed Abstract | CrossRef Full Text | Google Scholar

Hall, M. G., and Alexander, D. C. (2009). Convergence and parameter choice for monte-carlo simulations of diffusion MRI. IEEE Trans. Med. Imag. 28, 1354–1364. doi: 10.1109/TMI.2009.2015756

PubMed Abstract | CrossRef Full Text | Google Scholar

Harkins, K. D., Beaulieu, C., Xu, J., Gore, J. C., and Does, M. D. (2021). A simple estimate of axon size with diffusion MRI. Neuroimage 227, 117619. doi: 10.1016/j.neuroimage.2020.117619

PubMed Abstract | CrossRef Full Text | Google Scholar

Herrmann, L. R. (1976). Laplacian-isoparametric grid generation scheme. J. Eng. Mech. Div. 102, 749–756. doi: 10.1061/JMCEA3.0002158

CrossRef Full Text | Google Scholar

Hess, R. (2010). Blender Foundations: The Essential Guide to Learning Blender 2.6. Waltham, MA: Focal Press.

Google Scholar

Innocenti, G. M., and 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. doi: 10.1007/s00429-016-1266-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Jelescu, I. O., and Budde, M. D. (2017). Design and validation of diffusion mri models of white matter. Front. Phys. 28, 61. doi: 10.3389/fphy.2017.00061

PubMed Abstract | CrossRef Full Text | Google Scholar

Kakkar, L. S., Bennett, O. F., Siow, B., Richardson, S., Ianuş, A., Quick, T., et al. (2018). Low frequency oscillating gradient spin-echo sequences improve sensitivity to axon diameter: an experimental study in viable nerve tissue. Neuroimage 182, 314–328. doi: 10.1016/j.neuroimage.2017.07.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Korogod, N., Petersen, C. C., and Knott, G. W. (2015). Ultrastructural analysis of adult mouse neocortex comparing aldehyde perfusion with cryo fixation. Elife 4, e05793. doi: 10.7554/eLife.05793

PubMed Abstract | CrossRef Full Text | Google Scholar

Lam, S. K., Pitrou, A., and 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.

PubMed Abstract | Google Scholar

Lavdas, I., Behan, K. C., Papadaki, A., McRobbie, D. W., and Aboagye, E. O. (2013). A phantom for diffusion-weighted MRI (DW-MRI). J. Magnet. Reson. Imag. 38, 173–179. doi: 10.1002/jmri.23950

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H.-H., Jespersen, S. N., Fieremans, E., and Novikov, D. S. (2020). The impact of realistic axonal shape on axon diameter estimation using diffusion MRI. NeuroImage 223, 117228. doi: 10.1016/j.neuroimage.2020.117228

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H.-H., Yaros, K., Veraart, J., Pathan, J. L., Liang, F.-X., Kim, S. G., et al. (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.

PubMed Abstract | Google Scholar

Lewiner, T., Lopes, H., Vieira, A. W., and Tavares, G. (2003). Efficient implementation of marching cubes' cases with topological guarantees. J. Graph. Tools 8, 1–15. doi: 10.1080/10867651.2003.10487582

CrossRef Full Text | Google Scholar

Lin, C., Huang, Y., Quan, T., and Zhang, Y. (2018). Modelling brain-wide neuronal morphology via rooted cayley trees. Sci. Rep. 8, 1–10. doi: 10.1038/s41598-018-34050-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, A., Whittall, K., Adler, J., Li, D., Paty, D., and Graeb, D. (1994). In vivo visualization of myelin water in brain by magnetic resonance. Magnet. Reson. Med. 31, 673–677. doi: 10.1002/mrm.1910310614

PubMed Abstract | CrossRef Full Text | Google Scholar

Maier-Hein, K. H., Neher, P. F., Houde, J.-C., Côté, M.-A., Garyfallidis, E., Zhong, J., et al. (2017). The challenge of mapping the human connectome based on diffusion tractography. Nat. Commun. 8, 1349. doi: 10.1038/s41467-017-01285-x

PubMed Abstract | CrossRef Full Text | Google Scholar

McKerns, M. M., Strand, L., Sullivan, T., Fang, A., and Aivazis, M. A. (2012). Building a framework for predictive science. arXiv preprint arXiv:1202.1056. doi: 10.25080/Majora-ebaa42b7-00d

CrossRef Full Text | Google Scholar

Menzel, M., Michielsen, K., De Raedt, H., Reckfort, J., Amunts, K., and Axer, M. (2015). A jones matrix formalism for simulating three-dimensional polarized light imaging of brain tissue. J. Royal Soc. Interf. 12, 20150734. doi: 10.1098/rsif.2015.0734

PubMed Abstract | CrossRef Full Text | Google Scholar

Mingasson, T., Duval, T., Stikov, N., and Cohen-Adad, J. (2017). Axonpacking: an open-source software to simulate arrangements of axons in white matter. Front. Neuroinformat. 11, 5. doi: 10.3389/fninf.2017.00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, A., Gray, A., and Liu, T. (2003). New algorithms for efficient high dimensional non-parametric classification. Adv. Neural Inform. Process. Syst. 16.

Google Scholar

Murday, J. S., and Cotts, R. M. (1968). Self-diffusion coefficient of liquid lithium. J. Chem. Phys. 48, 4938–4945. doi: 10.1063/1.1668160

CrossRef Full Text | Google Scholar

Narvaez-Delgado, O., Rojas-Vite, G., Coronado-Leija, R., Ramírez-Manzanares, A., Marroquín, J. L., Noguez-Imm, R., et al. (2019). Histological and diffusion-weighted magnetic resonance imaging data from normal and degenerated optic nerve and chiasm of the rat. Data Brief 26, 104399. doi: 10.1016/j.dib.2019.104399

PubMed Abstract | CrossRef Full Text | Google Scholar

Neher, P. F., Laun, F. B., Stieltjes, B., and Maier-Hein, K. H. (2014). Fiberfox: facilitating the creation of realistic white matter software phantoms. Magnet. Reson. Med. 72, 1460–1470. doi: 10.1002/mrm.25045

PubMed Abstract | CrossRef Full Text | Google Scholar

Neuman, C. H. (1974). Spin echo of spins diffusing in a bounded medium. J. Chem. Phys. 60, 4508–4511. doi: 10.1063/1.1680931

CrossRef Full Text | Google Scholar

Nilsson, M., Lasič, S., Drobnjak, I., Topgaard, D., and 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. doi: 10.1002/nbm.3711

PubMed Abstract | CrossRef Full Text | Google Scholar

Nilsson, M., Lätt, J., Ståhlberg, F., van Westen, D., and Hagslätt, H. (2012). The importance of axonal undulation in diffusion mr measurements: a monte carlo simulation study. NMR Biomed. 25, 795–805. doi: 10.1002/nbm.1795

PubMed Abstract | CrossRef Full Text | Google Scholar

Novikov, D. S., Fieremans, E., Jespersen, S. N., and Kiselev, V. G. (2019). Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation. NMR Biomed. 32, e3998. doi: 10.1002/nbm.3998

PubMed Abstract | CrossRef Full Text | Google Scholar

Novikov, D. S., Kiselev, V. G., and Jespersen, S. N. (2018). On modeling. Magnet. Reson. Med. 79, 3172–3193. doi: 10.1002/mrm.27101

PubMed Abstract | CrossRef Full Text | Google Scholar

Parizel, P. M., Van Goethem, J., Özsarlak, Ö., Maes, M., and Phillips, C. (2005). New developments in the neuroradiological diagnosis of craniocerebral trauma. Eur. Radiol. 15, 569–581. doi: 10.1007/s00330-004-2558-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., et al. (2011). Scikit-learn: machine learning in Python. J. Machine Learn. Res. 12, 2825–2830.

Google Scholar

Plante, I., and 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. doi: 10.5772/53203

CrossRef Full Text | Google Scholar

Rafael-Patino, J., Girard, G., Truffet, R., Pizzolato, M., Thiran, J. P., and 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.

PubMed Abstract | Google Scholar

Rafael-Patino, J., Romascano, D., Ramirez-Manzanares, A., Canales-Rodríguez, E. J., Girard, G., and 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. doi: 10.3389/fninf.2020.00008

PubMed Abstract | CrossRef Full Text | Google Scholar

Romascano, D., Rafael-Patino, J., Jelescu, I., Barakovic, M., Tim, B., Jean-Philippe, D., et al. (2018). “Voxel size matters: big voxels are required to generate realistic extra-axonal DMRI signals from monte carlo simulations,” in ISMRM. 1–2.

Google Scholar

Schilling, K., Gao, Y., Janve, V., Stepniewska, I., Landman, B. A., and Anderson, A. W. (2017). Can increased spatial resolution solve the crossing fiber problem for diffusion MRI? NMR Biomed. 30, e3787. doi: 10.1002/nbm.3787

PubMed Abstract | CrossRef Full Text | Google Scholar

Schilling, K. G., Nath, V., Hansen, C., Parvathaneni, P., Blaber, J., Gao, Y., et al. (2019). Limits to anatomical accuracy of diffusion tractography using modern approaches. Neuroimage 185, 1–11. doi: 10.1016/j.neuroimage.2018.10.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneider-Mizell, C. M., Bodor, A. L., Collman, F., Brittain, D., Bleckert, A., Dorkenwald, S., et al. (2021). Structure and function of axo-axonic inhibition. Elife 10, e73783. doi: 10.7554/eLife.73783

PubMed Abstract | CrossRef Full Text | Google Scholar

Sepehrband, F., Alexander, D. C., Clark, K. A., Kurniawan, N. D., Yang, Z., and Reutens, D. C. (2016). Parametric probability distribution functions for axon diameters of corpus callosum. Front. Neuroanat. 10, 59. doi: 10.3389/fnana.2016.00059

PubMed Abstract | CrossRef Full Text | Google Scholar

Shekhar, R., Fayyad, E., Yagel, R., and 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.

Google Scholar

Simpson, J., and Carr, H. (1958). Diffusion and nuclear spin relaxation in water. Phys. Rev. 111, 1201. doi: 10.1103/PhysRev.111.1201

CrossRef Full Text | Google Scholar

Söderman, O., and Jönsson, B. (1995). Restricted diffusion in cylindrical geometry. J. Magnet. Reson. Ser. A 117, 94–97. doi: 10.1006/jmra.1995.0014

CrossRef Full Text | Google Scholar

Sorkine, O., Cohen-Or, D., Lipman, Y., Alexa, M., Rössl, C., and Seidel, H.-P. (2004). “Laplacian surface editing,” in Proceedings of the 2004 Eurographics/ACM SIGGRAPH Symposium on Geometry Processing (ACM), 175–184.

Google Scholar

Sporns, O. (2011). The non-random brain: efficiency, economy, and complex dynamics. Front. Comput. Neurosci. 5, 5. doi: 10.3389/fncom.2011.00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Stanisz, G. J., Wright, G. A., Henkelman, R. M., and Szafer, A. (1997). An analytical model of restricted diffusion in bovine optic nerve. Magnet. Reson. Med. 37, 103–111. doi: 10.1002/mrm.1910370115

PubMed Abstract | CrossRef Full Text | Google Scholar

Stejskal, E. O., and Tanner, J. E. (1965). Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J. Chem. Phys. 42, 288–292. doi: 10.1063/1.1695690

PubMed Abstract | CrossRef Full Text | Google Scholar

Stikov, N., Campbell, J. S., Stroh, T., Lavelée, M., Frey, S., Novek, J., et al. (2015). In vivo histology of the myelin g-ratio with magnetic resonance imaging. Neuroimage 118, 397–405. doi: 10.1016/j.neuroimage.2015.05.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Stroustrup, B. (1999). “An overview of the C++ programming language,” in Handbook of Object Technology, ed S. Zamir (Boca Raton, FL: CRC Press), 72.

Google Scholar

Sullivan, B., and Kaszynski, A. (2019). PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK). J. Open Sourc. Softw. 4, 1450. doi: 10.21105/joss.01450

CrossRef Full Text | Google Scholar

Syková, E., and Nicholson, C. (2008). Diffusion in brain extracellular space. Physiol. Rev. 88, 1277–1340. doi: 10.1152/physrev.00027.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

To, X. V., Mohamed, A. Z., Cumming, P., and 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. doi: 10.1016/j.bbi.2022.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Tønnesen, J., Inavalli, V. K., and Nägerl, U. V. (2018). Super-resolution imaging of the extracellular space in living brain tissue. Cell 172, 1108–1121. doi: 10.1016/j.cell.2018.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Tournier, J.-D., Calamante, F., and Connelly, A. (2007). Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution. NeuroImage 35, 1459–1472. doi: 10.1016/j.neuroimage.2007.02.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Tournier, J. D., Yeh, C. H., Calamante, F., Cho, K. H., Connelly, A., and Lin, C. P. (2008). Resolving crossing fibres using constrained spherical deconvolution: Validation using diffusion-weighted imaging phantom data. NeuroImage 42, 617–625. doi: 10.1016/j.neuroimage.2008.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Truffet, R., Rafael-Patino, J., Girard, G., Pizzolato, M., Barillot, C., Thiran, J.-P., et al. (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.

Google Scholar

Tuch, D. S. (2004). Q-ball imaging. Magnet. Reson. Med. 52, 1358–1372. doi: 10.1002/mrm.20279

PubMed Abstract | CrossRef Full Text | Google Scholar

Turau, V. (1991). Fixed-radius near neighbors search. Informat. Process. Lett. 39, 201–203. doi: 10.1016/0020-0190(91)90180-P

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Walt, S., Schönberger, J. L., Nunez-Iglesias, J., Boulogne, F., Warner, J. D., Yager, N., et al. (2014). Scikit-image: image processing in Python. PeerJ 2, e453. doi: 10.7287/peerj.preprints.336v2

CrossRef Full Text | Google Scholar

van Gelderen, P., de Vleeschouwer, M. H. M., DesPres, D., Pekar, J., van Zijl, P. C. M., and Moonen, C. T. W. (1994). Water diffusion and acute stroke. Magnet. Reson. Med. 31, 154–163. doi: 10.1002/mrm.1910310209

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Rossum, G., and Drake, F. L. (2009). Python 3 Reference Manual. Scotts Valley, CA: CreateSpace.

Google Scholar

Van Verth, J. M., and Bishop, L. M. (2015). Essential Mathematics for Games and Interactive Applications. Boca Raton, FL: CRC Press.

Google Scholar

Veraart, J., Nunes, D., Rudrapatna, U., Fieremans, E., Jones, D. K., Novikov, D. S., et al. (2020). Noninvasive quantification of axon radii using diffusion MRI. Elife 9, e49855. doi: 10.7554/eLife.49855

PubMed Abstract | CrossRef Full Text | Google Scholar

Veraart, J., Raven, E. P., Edwards, L. J., Weiskopf, N., and Jones, D. K. (2021). The variability of mr axon radii estimates in the human white matter. Hum. Brain Map. 42, 2201–2213. doi: 10.1002/hbm.25359

PubMed Abstract | CrossRef Full Text | Google Scholar

Warner, W., Palombo, M., Cruz, R., Callaghan, R., Shemesh, N., Jones, D. K., et al. (2023). Temporal diffusion ratio (TDR) for imaging restricted diffusion: optimisation and pre-clinical demonstration. NeuroImage 269, 119930. doi: 10.1016/j.neuroimage.2023.119930

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Schneider, T., Wheeler-Kingshott, C. A., and Alexander, D. C. (2012). NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain. Neuroimage 61, 1000–1016. doi: 10.1016/j.neuroimage.2012.03.072

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, F.-L., Li, Z., Gough, J. E., Cristinacce, P. L. H., and Parker, G. J. (2018). Axon mimicking hydrophilic hollow polycaprolactone microfibres for diffusion magnetic resonance imaging. Mater. Design 137, 394–403. doi: 10.1016/j.matdes.2017.10.047

PubMed Abstract | CrossRef Full Text | Google Scholar