Exploring white matter dynamics and morphology through interactive numerical phantoms: the White Matter Generator

Brain white matter is a dynamic environment that continuously adapts and reorganizes in response to stimuli and pathological changes. Glial cells, especially, play a key role in tissue repair, inflammation modulation, and neural recovery. The movements of glial cells and changes in their concentrations can influence the surrounding axon morphology. We introduce the White Matter Generator (WMG) tool to enable the study of how axon morphology is influenced through such dynamical processes, and how this, in turn, influences the diffusion-weighted MRI signal. This is made possible by allowing interactive changes to the configuration of the phantom generation throughout the optimization process. The phantoms can consist of myelinated axons, unmyelinated axons, and cell clusters, separated by extra-cellular space. Due to morphological flexibility and computational advantages during the optimization, the tool uses ellipsoids as building blocks for all structures; chains of ellipsoids for axons, and individual ellipsoids for cell clusters. After optimization, the ellipsoid representation can be converted to a mesh representation which can be employed in Monte-Carlo diffusion simulations. This offers an effective method for evaluating tissue microstructure models for diffusion-weighted MRI in controlled bio-mimicking white matter environments. Hence, the WMG offers valuable insights into white matter's adaptive nature and implications for diffusion-weighted MRI microstructure models, and thereby holds the potential to advance clinical diagnosis, treatment, and rehabilitation strategies for various neurological disorders and injuries.

1 Introduction

The neuronal network of the brain is a dynamic structure, constantly adapting to internal and external stimuli. Within the white matter, the local environment of neurons is found to modulate their axon morphology (Andersson et al., 2020). Especially glial cells play a significant role in this modulation. While the glial cells do not directly participate in neuronal signaling, they carry out crucial supporting tasks: myelination (oligodendrocytes), maintaining an appropriate chemical environment for neuronal signaling (astrocytes), and removing cellular debris from sites of injury or normal cell turnover (microglial) (Purves et al., 2001). To carry out these tasks, the glial cells move through the tissue. Such dynamic behavior has been observed with microscopy (Davalos et al., 2005; Nimmerjahn et al., 2005; Tønnesen et al., 2018). Furthermore, it has been found that the morphology of axons adapts in response to these environmental changes (Andersson et al., 2020). Moreover, there is a relation between the morphology and function of the individual neurons. Morphological factors that affect the axonal conduction velocity include the axon diameter (Hursh, 1939; Skoven et al., 2023), the myelin sheath thickness (Sanders and Whitteridge, 1946), and the length and spacing of nodes of Ranvier (Arancibia-Cárcamo et al., 2017). Quantifying the morphology of axons and brain white matter in general is thus of great interest for enhancing our understanding of various brain processes, and can provide a potential biomarker of disease.

The most promising method for assessing axon morphology non-invasively, is diffusion-weighted MRI (dMRI). dMRI reflects the morphological properties of the underlying tissue by measuring diffusion properties of water molecules within the tissue (Basser et al., 1994; Beaulieu, 2002; Alexander et al., 2019; Novikov et al., 2019). To advance the method even further, we need access to extensive experimentation for scan sequence optimization and biophysical model validation informed by known ground truths. This is unfeasible to obtain with either preclinical or clinical MRI due to the limited spatial resolution.

It is therefore of great interest to supplement dMRI research with numerical simulations, because these come with a ground truth. dMRI can with advantage be simulated based on Monte Carlo (MC) diffusion simulations (Hall and Alexander, 2009; Rafael-Patino et al., 2020). The level of realism in simulations is crucial for the specificity and applicability obtainable by the models developed based on them (Nilsson et al., 2012; Andersson et al., 2020; Brabec et al., 2020; Lee et al., 2020).

Mimicking the highly complex tissue of white matter is a great challenge. White matter has commonly been represented as idealized straight, coaxial, infinitely long, and non-touching cylinders (Novikov et al., 2019). However, from recent advances in 3D histology, it has been validated that the white matter has a much more complex configuration (Tønnesen et al., 2018; Abdollahzadeh et al., 2019; Lee et al., 2019; Andersson et al., 2020)(Abdollahzadeh et al., 2021). While individual axons are found to express non-circular cross-sections, longitudinal diameter variations, and tortuosity, axons on the ensemble scale are found to express orientation dispersion and crossing fibers. Furthermore, recent studies show that such characteristics have a crucial impact on the diffusion signal and modeling (Nilsson et al., 2012; Andersson et al., 2020; Brabec et al., 2020; Lee et al., 2020; Winther et al., 2023), and should therefore be incorporated into our simulations to improve the realism.

Two primary approaches have been taken to generate realistic numerical white matter phantoms: segmentation (both semiautomatic and automatic) from 3D histology of tissue (Abdollahzadeh et al., 2019; Lee et al., 2019, 2021; Andersson et al., 2020), and numerical synthesis (Balls and Frank, 2009; Hall and Alexander, 2009; Budde and Frank, 2010; Landman et al., 2010; Mingasson et al., 2017; Ginsburger et al., 2019; Callaghan et al., 2020; Villarreal-Haro et al., 2023). Segmentation of white matter tissue provides a very high degree of realism. However, it can be very resource-consuming with respect to time for manual editing and the sacrifice of animal lives. Furthermore, this approach provides little flexibility in the morphological variation of the tissue and only provides static snapshots thereof. In contrast, numerical synthesis allows for much lower resource consumption and allows for a high degree of morphological flexibility. However, the outputs are based on assumptions of what the anatomy looks like (Dyrby et al., 2018).

Various tools have been developed for numerical synthesis of white matter phantoms with complex morphological properties. Recent tools include MEDUSA (Ginsburger et al., 2019), ConFiG (Callaghan et al., 2020), and CACTUS (Villarreal-Haro et al., 2023). MEDUSA (Ginsburger et al., 2019) enables the representation of the most diverse tissue elements by allowing multiple compartments including axons, astrocytes, oligodendrocytes, nodes of Ranvier, and myelinated axons. All structures are represented as spheres. While the sphere representation allows for high representational power of the longitudinal axon morphology characteristics, it does not allow for eccentric cross-sections documented by histology (Abdollahzadeh et al., 2019; Lee et al., 2019; Andersson et al., 2020). Meanwhile, both ConFiG (Callaghan et al., 2020) and CACTUS (Villarreal-Haro et al., 2023) allow only for the inclusion of myelinated axons. However, due to a higher refinement of the axonal cross-sections, these possess higher realism compared with the otherwise circular cross-sections. CACTUS (Villarreal-Haro et al., 2023) stands out for its superior computational efficiency, which in turn enables the creation of larger and more dense phantoms. However, none of these tools considers the dynamic aspects of white matter tissue which has a crucial influence on axon morphology (Andersson et al., 2020). In this work, we present the White Matter Generator (WMG); a new tool for generating interactive numerical phantoms for Monte Carlo dMRI simulations. The interactivity of the WMG tool enables the mimicking of brain white matter dynamics. The phantoms can consist of multiple compartments: unmyelinated axons, myelinated axons, (here defined as fibers), and cells, separated by extra-cellular space. The key assumption for the synthesis is that axon morphology is modulated by the local environment, as observed with XNH imaging in our previous work (Andersson et al., 2020, 2022). Due to morphological flexibility and computational advantages during the optimization, the tool uses ellipsoids as building blocks for all structures during the synthesis; chains of ellipsoids for fibers, and individual ellipsoids for cells. Thereby, fibers can obtain non-circular cross-sections, longitudinal diameter variations, and tortuosity as observed in 3D histology. To mimic a dynamic tissue environment, the tool allows the user to generate phantoms at consecutive time points by interactively changing parameters and tissue composition during the optimization process. The output is in the format of PLY surface meshes, which makes them directly compatible with existing Monte Carlo diffusion simulators such as the widely used MC-DC Simulator (Rafael-Patino et al., 2020) and Camino (Hall and Alexander, 2009). We demonstrate examples of various tissue configurations: varying degrees of fiber dispersion, types of bundle crossings, demyelination, inclusion of static cells, and cell dynamics. The biological relevance of the outputted axons is evaluated based on a set of morphological metrics and compared with those observed with XNH imaging of monkey brain white matter in our previous work (Andersson et al., 2020, 2022).

2 Methods

Phantoms can consist of myelinated axons, axons, and cells, separated by extra-cellular space. To obtain morphological flexibility and computational advantages during the optimization, the tool uses ellipsoids as the building blocks for all structures during the synthesis; chains of ellipsoids for fibers, and individual ellipsoids for cells.

The optimization is based on a force-biased packing algorithm, first introduced by Altendorf and Jeulin (2011), where the phantoms are obtained as an equilibrium between repulsion forces (for avoiding overlap between individual axons and cells) and recover forces (for ensuring the structure of the individual axons) between the ellipsoidal building blocks.

Different brain regions and types of cell dynamics can be mimicked by adapting the configuration w.r.t. fiber diameter distributions, fiber volume fractions (FVF), cell volume fractions (CVF), global dispersion (ϵ), bundle crossings, and interactive changes. Changes to the configuration can be made interactively and at any time during the optimization.

An optimization starts with one configuration file (config-file) and ends with another updated config-file. Thereby, the updated config file can be easily adapted and given as input for another round of optimization according to these adaptations. A conceptual flow chart is seen in Figure 1, and each compartment is described in more detail later in this section.

Figure 1

Figure 1. The WMG tool can be divided into four compartments. Configuration: A config-file is generated based on a set of phantom parameters and optimization parameters. Within the config-file, all structures are represented by ellipsoids. Optimization: The optimization runs based on the config-file, and is carried out by iterating over the axonal ellipsoids. Cellular ellipsoids remain static. Maintaining consistency between the input and output format makes interactive adaptation of the configuration convenient. Re-configuration: After a round of optimization, the config-file can be adapted by changing phantom parameters (e.g., CVF and cell positions), and/or optimization parameters (e.g.,ellipsoidDensity, growSpeed and maxIterations). It can then be used as input for another round of optimization. Post-processing: After a round of optimization, the config-file can be post-processed to obtain a mesh-representation from the ellipsoid-representation. Furthermore, it is often beneficial to perform a Garland Heckbert simplification of these meshes to remove redundant vertices to reduce the computational load of e.g., Monte Carlo diffusion simulations.

The WMG tool is written in TypeScript, whereas supporting functions for configuration and post-processing are written in Python 3. A simplified version of the WMG tool is available as a graphical web interface at https://map-science.github.io/WhiteMatterGenerator. This version allows one to get a more intuitive feeling of the tool and its parameters while testing out different configurations. The full version is available as a command line interface which can be installed by following the instructions at https://map-science.github.io/WhiteMatterGenerator/help/cli. This version enables automation and large-scale production, and this is the version we will be focusing on here. The supporting functions for configuration and post-processing are available at https://github.com/MaP-science/WhiteMatterGenerator.

2.1 Configuration

There are two types of parameters: phantom parameters, which specify the content and limits of the phantom, and optimization parameters, which specify the optimization criteria of the phantom. As long as the config-file follows the correct formatting, the user can configure the phantom in any way they like. Below is a guide to the options provided with the WMG.

The phantom parameters are:

• Voxel:

- Outer (larger) voxel dimensions: Defines the spatial restriction of phantom content as a cuboid.

- Inner (smaller) voxel dimensions: Defines the cuboid volume wherein the FVF is optimized and represents the field-of-view of the phantom. This voxel is necessary because the discontinuation at the outer voxel boundary will lead to biased morphology at that boundary.

• Fibers:

- Initial base point (b) and direction (r): fibers are initialized as chains of ellipsoids arranged in straight lines defined with a direction and one fixed base point. The fixed points are uniformly distributed along the axis of primary fiber-orientation .

- Global fiber dispersion (ϵ): fiber directions are uniformly sampled on the surface of a spherical cone, where the height of the cone cap is scaled by ϵ. ϵ = 0.0 allows disperson of 0 deg and ϵ = 1.0 allows dispersion of 90 deg.

- Crossing bundles: Can be configured either as adjacent sheets (Kjer et al., 2023) or interwoven (see Figure 5). Specified by the number of bundles, their primary direction and the fraction of fibers they contain.

- fiber-ellipsoid                                                                 separation (mapFromMaxDiameterToEllipsoidSeparation): The ellipsoid separation within each fiber is defined from a user-specified mapping function which takes the maximum allowed fiber-diameter as input.

- Allowed diameter ranges for the individual fibers (mapFromMaxDiameterToMinDiameter):        To enable closer packing of the axons, it is possible to allow the diameter of each ellipsoid to vary. The range is defined by a user-specified mapping function which takes the maximum allowed fiber-diameter as input. We set the allowed diameter range of the individual axon as a margin around its targeted diameter. The targeted diameters come from an idealized packing of straight parallel cylinders with diameters that follow the targeted gamma-distribution. This initial packing is obtained with the CylinderGammaDistribution() function of the MC-DC Simulator (Rafael-Patino et al., 2020). By initializing all fibers with a diameter which is significantly lower than that of the idealized packing, the axons are able to move more freely during the initial iterations of the optimizations and thereby achieve a more dense packing.

- g-ratio: When the ellipsoid-representation is converted to the mesh-representation the fiber is compartmentalized into axon and myelin based on this g-ratio. That is, the myelin diameter matches that of the fiber, while the axon diameter is scaled in relation to the fiber's diameter through the g-ratio.

- Contract speed (contractSpeed): How much the fibers contract per step, i.e., how stiff the axons are. This number should be ≥0.

- Deformation                                                     factor (mapFromDiameterToDeformationFactor):     How much the fiber-ellipsoids should deform (change diameter and eccentricity) as opposed to change position when a collision occurs. The deformation is a number between 0 and 1, and is defined by a user-specified mapping function which takes the current fiber-diameter as input. A deformation factor of 0 means that the axonal ellipsoid cannot be deformed at all and hence remains as is. A deformation factor of 1 means that the ellipsoid will deform as much as possible rather than change the position of its center. By adjusting this map between stages of an optimization scheme, the diameters and cross-sections of axons can be fixed such that only the axon trajectories will adapt according to later cell mobility.

• Cells:

- Target CVF (targetCVF): Cells can be added based on ellipsoidal dimensions randomly sampled from normal distributions. One cell is sampled and placed at a time. Cells are added until the targetCVF is met. When placed, a cell may overlap with axons but not cells. In the following optimization, the axons will then adapt to avoid any overlap.

- Position: Defines the center of the cell.

- Shape: Specifies the 3 × 3 transformation matrix used when going from a unit sphere to an ellipsoid. This ellipsoid will be the shape of the cell.

• Output format:

- Radial resolution of outputted meshes: Determines how detailed the output mesh will be.

- Extend fibers around the voxel: Defines whether or not to extend the fibers around the voxel to create mirrored intra-axonal compartments.

The optimization parameters are:

• Target FVF (targetFVF): The targeted FVF.

• Grow speed (growSpeed): How much the fibers grow per optimization step. 0 means no growth, and 1 means that the axon will grow to 100% of its target size in 1 step.

• Maximum number of iterations (maxIterations): Even if the targetFVF is not reached at this point, the optimization will terminate here.

• Output interval (outputInterval): Interval with which an output with the current configuration is provided.

It is beneficial to apply optimization schemes that start out with higher growSpeed (bigger steps) and lower intra-fiber ellipsoid density (less computationally demanding), and go toward lower growSpeed (smaller steps) and higher intra-fiber ellipsoid resolution (more computationally demanding). An example is illustrated in Figure 2. Furthermore, it is crucial to choose realistic parameters to obtain a successful phantom optimization.

Figure 2

Figure 2. Example of an optimization scheme. At initialization, the phantom consists of an outer-voxel (for spatial restriction), an inner voxel (for the computing optimization metric FVF), straight fibers, and cells. Initially, smaller diameters are assigned to the fibers to allow for more flexibility as they grow into their converged morphology. (A) Inter-fiber ellipsoids grow step-wise toward their maximum allowed diameter while adapting to their surroundings as a consequence of collisions with these. (B) The fiber-ellipsoid density is increased. (C) Growth is repeated. (D) Steps (B) and (C) or variations thereof can be repeated. (E) Once optimized to the user's liking, the ellipsoid-representation can be converted to a mesh-representation and then applied in e.g., Monte Carlo diffusion simulations.

2.2 Optimization

2.2.1 Initialization of voxels, fibers, and cells

Each fiber is initialized with a base point, b, and a direction, r. A ray is traced along r and −r until the voxel boundary is hit. The two points hit by the ray define the endpoints of the fiber. The straight line connecting the two endpoints is then filled with ellipsoids according to the specified ellipsoid density and minimum fiber diameter. Each cell is represented by a single ellipsoid. They cannot move or deform, but otherwise act in the same way as fibers.

Each ellipsoid has two properties: a position vector and a shape matrix. With no deformation, the shape of an ellipsoid is a unit sphere S² = {q ∈ ℝ³| ||q|| = 1}. The deformation of the sphere is represented by a 3 × 3 matrix, S, such that the surface of the resulting ellipsoid is E = {Sq + p| ||q|| = 1} where p is the position.

Initially, the matrix S = μI where μ is given by the allowed minimum radius for the given fiber. For increased flexibility, μ should be very small relative to the voxel size.

Hence, each ellipsoid, j, on each fiber, i, is represented by a position vector, p_{i, j}, and a shape matrix, S_{i, j}, and the surface of the ellipsoid is thus defined by

$\begin{array}{ll}{E}_{i,j}=\left\{{\text{S}}_{i,j}\text{q}+{\text{p}}_{i,j}|||\text{q}||=1\right\}.& (1)\end{array}$

Ellipsoids are deformed by applying a deformation matrix D = srr^T +I . Multiplying (the points on) a surface by D thereby results in a scaling of 1 + s along r where s ∈ ℝ and r ∈ ℝ³, ||r|| = 1. This can be seen by taking an arbitrary point p and transforming it using D

$\begin{array}{ll}\text{D}\text{p}=(s\text{r}{\text{r}}^T+\text{I})\text{p}=s(\text{r}{\text{r}}^T)\text{p}+\text{p}=s({\text{r}}^T\text{p})\text{r}+\text{p}.& (2)\end{array}$

After n deformations of an ellipsoid the shape matrix will have the form S = D_n...D₂D₁.

2.2.2 Growth, contraction, and redistribution for fiber-ellipsoids

The growth step considers each ellipsoid j of each fiber i. The shape matrix of the updated ellipsoid S_new is set to a weighted average between the current shape S_old and the target shape (i.e., a sphere having the maximum radius r_max) by

$\begin{array}{ll}{\text{S}}_{\text{new}}=(1-a){\text{S}}_{\text{old}}+a{r}_{\text{max}}\text{I}& (3)\end{array}$

where a ∈ [0, 1] is the growSpeed.

A contraction step adjusts the curve on which the ellipsoids lie. The force along the fiber is dependent on the contraction coefficient κ. A value of κ = 0 will do nothing, and a value of κ = 1 will move the position p_{i, j} of the ellipsoid all the way over to $\text{c}=\frac{1}{2}({\text{p}}_{i,j-1}+{\text{p}}_{i,j+1})$ making this segment of the fiber a straight line. For any value of κ∈[0, 1] the position is updated to (1−κ)p_{i, j}+κc.

A redistribution step is then performed to prevent ellipsoids from clumping together and thereby creating gaps in the ellipsoid chain. This is done by updating the position of each ellipsoid in a chain such that they are evenly distributed along the trajectory of the fiber. Furthermore, the endpoints of the fibers are moved to the nearest side whenever they are inside the voxel to make sure that the fiber spans the entire voxel.

2.2.3 Collisions: fiber-voxel, fiber-fiber, and fiber-cell

Fiber-voxel collisions (outer voxel) are checked for each ellipsoid independently. Whenever an ellipsoid's center is outside of the voxel it is projected back onto the boundary. fiber-fiber collisions are checked by computing the potential overlap between each ellipsoid of a fiber and all other ellipsoids of all other fibers. fiber-cell collisions are checked by likewise computing the potential overlap between each ellipsoid of a fiber and all cells.

To compute overlaps between ellipsoids (see Figure 3), we utilize that finding the surface point of an ellipsoid furthest away along a vector r is the same as finding the surface point which has the normal r. When deforming an ellipsoid using the matrix D, the normals are transformed by D⁻¹. Since after n deformations D_i, i ∈ {1, ..., n} the shape has the form S = D_n...D₂D₁, the corresponding transformation matrix N for the normals is

$\begin{array}{r}\text{N}={\text{D}}_n^{-1}...{\text{D}}_2^{-1}{\text{D}}_1^{-1}={({({\text{D}}_n^T...{\text{D}}_2^T{\text{D}}_1^T)}^T)}^{-1}=\\ {({({\text{D}}_n...{\text{D}}_2{\text{D}}_1)}^T)}^{-1}={({\text{S}}^T)}^{-1}& (4)\end{array}$

since each D_i is symmetric. When going from a normal on an ellipsoid to the corresponding normal of the original sphere, the normal should be transformed by the inverse, i.e., N⁻¹ = S^T. After this transformation, the resulting vector can be normalized. Since we now have a unit normal on a unit sphere, the position on that sphere is given by the same vector. To get the corresponding position on the ellipsoid, one can simply multiply by S. Hence the point with the normal r on the ellipsoid with shape S and center at the origin is $\text{x}(\text{r})=\text{S}\frac{{\text{S}}^T\text{r}}{||{\text{S}}^T\text{r}||}$. Given the ellipsoid j of fiber i and a direction r the extremal point is thus e_{i, j}(r) = p_{i, j} + x_{i, j}(r). This extremum is used to compute the overlap between two ellipsoids.

Figure 3

Figure 3. 2D visualization of the concept applied for checking for potential overlap of two ellipsoids. In the example, we have two ellipsoids centered at p₁ and p₂, respectively. The objective is to check if there exists a set of hyperplanes, P₁(r) and P₂(−r), which separates the ellipsoids. This is done by finding the surface point e₁(r) for which the normal is parallel with r, and the surface point e₂(−r) for which the normal is parallel with −r. Now, if the dot product between the vector r^T [lying in the plane P₁(r)] and the vector e₁(r) − e₂(−r) is <0, the planes P₁(r) and P₂(−r) must separate the two ellipsoids. i.e., if there exists an r such that r^T · (e₁(r) − e₂(−r)) < 0, the two ellipsoids are not overlapping.

Define the overlap of two ellipsoids j₁ and j₂ belonging to the fibers i₁ and i₂ respectively along an axis r as o_{(i1, j1), (i2, j2)}$(\text{r})={\text{r}}^T({\text{e}}_{i_1,j_1}(\text{r})-{\text{e}}_{i_2,j_2}(-\text{r}))$. If one can find a separating axis (Gottschalk et al., 1997), r, such that o_{(i1, j1), (i2, j2)}(r) < 0, the ellipsoids do not overlap. If o_{(i1, j1), (i2, j2)}(r) > 0 for all values of r, the ellipsoids overlap, since they are convex. To check whether two ellipsoids overlap, one can find the minimum overlap and check if it is positive or negative. The ellipsoids overlap if and only if $\underset{\text{r}}{min}{o}_{(i_1,j_1),(i_2,j_2)}(\text{r})>0$. The minimum is computed numerically starting with the initial guess r = p_{i2, j2} − p_{i1, j1}.

When the axis r providing the least overlap o_{(i1, j1), (i2, j2)}(r) is found, it is checked if o_{(i1, j1), (i2, j2)}(r) < 0. If so, nothing should be done. In the case where o_{(i1, j1), (i2, j2)}(r) > 0 the ellipsoids should be updated such that the new value of o_{(i1, j1), (i2, j2)}(r) is 0. This is done by updating both the positions and shapes of the ellipsoids depending on the deformation factor. In practice, a small constant is added to o to make sure that there is a minimum distance between fibers in order to avoid overlaps in the final mesh-representation.

The overlaps are computed at each iteration of the algorithm. In most cases, the overlap will be negative since most ellipsoids are far apart. To avoid the need for computing the overlap between all pairs of ellipsoids, a hierarchy of axis-aligned bounding boxes (AABB) is computed for each axon. If the AABBs of two axons do not overlap, then the overlap computation of all the pairs of ellipsoids between these two axons can be skipped. If the AABBs do overlap, the AABBs will be recursively split into smaller AABBs and checked for overlap until the bottom of the hierarchy is reached. In this case, the ellipsoid overlap is computed.

2.2.4 Updating shapes and positions of fiber-ellipsoids

The key assumption, that axon morphology is influenced by the local environment (Andersson et al., 2020, 2022), is implemented by requiring each axonal ellipsoid to adapt both shape and position according to its local environment in order to avoid overlapping structures. Each ellipsoid will have a specified deformation factor δ(||x_{i, j}(r)||) ∈ [0, 1] depending on the specified deformation factor map. Each fiber has a minimum allowed radius μ_i specifying that it should hold for any ellipsoid j that ||x_{i, j}(r)|| ≥ μ_i. The overlap o_{(i1, j1), (i2, j2)}(r) is divided in two equally big parts, and each part is handled by each of the two ellipsoids. So the deformation of an ellipsoid will result in a decrease in its size of at most $\frac{1}{2}{o}_{(i_1,j_1),(i_2,j_2)}(\text{r})$. This will ensure that the ellipsoids will not have a negative overlap after being deformed. Since multiplying by a deformation matrix corresponds to a scaling along r, the distance that the extremal point should be moved has to be divided by the current size along r to get the scaling factor s. So in the case with ellipsoid j₁ of fiber i₁ the value of s would be $s_1=\frac{-\delta (||{\text{x}}_{i_1,j_1}(\text{r})||)\frac{1}{2}{o}_{(i_1,j_1),(i_2,j_2)}(\text{r})}{||{\text{x}}_{i_1,j_1}(\text{r})||}=-\frac{\delta (||{\text{x}}_{i_1,j_1}(\text{r})||)o_{(i_1,j_1),(i_2,j_2)}(\text{r})}{2||{\text{x}}_{i_1,j_1}(\text{r})||}$. In the case where this deformation would make the ellipsoid's size smaller than the minimum μ_i, a less negative value of s is chosen such that the resulting size is exactly μ_i: $s_2=\frac{{\mu }_i-||{\text{x}}_{i_1,j_1}(\text{r})||}{||{\text{x}}_{i_1,j_1}(\text{r})||}=\frac{{\mu }_i}{||{\text{x}}_{i_1,j_1}(\text{r})||}-1$. So the resulting value of s is s = max(s₁, s₂). So the deformation matrix D₁ for ellipsoid j₁ is ${\text{D}}_1=max(s_1,s_2)\text{r}{\text{r}}^T+\text{I}$. Now the shape matrix can be updated by calculating D₁S_{i1, j1}. The equivalent computations are performed for ellipsoid j₂.

Secondly, the positions of the ellipsoids are updated. Each ellipsoid is treated as if they had equal “mass”, i.e., the center of their positions is conserved. The total distance to move the ellipsoids apart is equal to the overlap o_{(i1, j1), (i2, j2)}(r) after the shapes have been updated. By moving the two ellipsoids by the same amount, the new positions are updated to ${\text{p}}_{i_l,j_l}+{(-1)}^l\frac{1}{2}{o}_{(i_1,j_1),(i_2,j_2)}(\text{r})\text{r}$ for l ∈ {1, 2}.

2.2.5 Computational complexity

In the worst case, all ellipsoids are close to each other and the number of overlaps that need to be computed will be O(ellipsoidCount²). Since the number of ellipsoids is proportional to (1) the number of axons, (2) the voxel size, and (3) the ellipsoid density, this could also be written as O((axonCount·voxelSize·ellipsoidDensity)²).

2.3 Post-processing

A phantom is outputted in two formats: an ellipsoid-representation and a mesh-representation. The ellipsoid-representation is the updated config-file where the shape and position of each ellipsoid is specified. This config-file can be edited and used as input for further optimization. The mesh-representation is acquired from the ellipsoid-representation by generating tube-like mesh-segments to connect all ellipsoidal cross-sections within a fiber. The user specifies the number of radial segments that the tubes should have (radial resolution). The number of longitudinal segments is equal to the number of ellipsoids on the fiber. The cells are exported as they are—i.e., meshes having the shape of the corresponding ellipsoids. The mesh-representation can be outputted as individual files for each fiber and cell, or as one combined file containing all elements.

2.3.1 Optimization of meshes for Monte Carlo diffusion simulations

The higher the ellipsoid-density, the smaller deviation is reached between the ellipsoid-representation and the mesh-representation due to gaps between ellipsoids. Hence, a high ellipsoid density is required to avoid overlapping fibers after conversion to the mesh-representation. This means that the resulting meshes have a likewise high longitudinal resolution. Especially if the meshes are intended for Monte Carlo simulations of dMRI, this is unfavorable since each mesh element comes with a computational cost due to the extensive collision-detection involved. Meanwhile, the high mesh resolution is often redundant.

By performing Garland Heckbert simplification (Garland and Heckbert, 1997) of each myelin and axon mesh, we reduce the number of faces and vertices significantly. The reduction depends on the initial resolution and the allowed deviation between the initial and the simplified mesh. We allow a deviation of minDistance/2 − 0.0005 μm, where minDistance is the minimum distance allowed between ellipsoids of different fibers. For a radial resolution of 16 segments and ellipsoidDensityScaler of 0.20, the Garland Heckbert simplification decreases the number of faces by around 50%.

Before employing the phantoms in Monte Carlo diffusion simulations, the ends of the meshes are sealed to obtain “water tightness”.

The intra-axonal compartment can be extended by making two extra copies of the ellipsoids mirrored through each of the fiber's endpoints without the requirement of further optimization.

2.4 Morphological analysis

To demonstrate the biological relevance of the morphological features expressed in the numerical phantoms generated with the WMG, we compare the morphology of a set of numerical phantoms with real axons. The real axons originate from a monkey brain, and have been quantified by segmentation from XNH-volumes in previous work (Andersson et al., 2020). A total of 58 axons with lengths between 120 and 304 μm were segmented. The meshes are available at www.drcmr.dk/powder-averaging-dataset. The morphological metrics quantify axon diameter variations, cross-sectional eccentricity, and tortuosity. All metrics are sampled with ~375 nm spacing along axon trajectories (some variation due to tortuosity).

2.4.1 Axon diameters

We quantify axon diameter (AD) by the equivalent diameter as in Abdollahzadeh et al. (2019). i.e., axon diameters reported here are the diameter of a circle with an area equal to that of the axonal cross-section perpendicular to its local trajectory. For the real axons, the metric is extracted from the segmentation. This method is described in more detail in Andersson et al. (2020). For the WMG-generated axons, the metric is extracted from the ellipsoid-representation.

For individual axons, mean AD [mean(AD)] is calculated over all measurements along an axon, and the standard deviation of the AD [std(AD)] quantifies the variation of these values.

2.4.2 Cross-sectional eccentricity

We quantify cross-sectional eccentricity based on elliptic parameterization by $e=\sqrt{1-b^2/a^2}$ where a is the major axis and b is the minor axis (see Figure 10). For the real axons, the circumference of the segmentation is sampled by 24 points. The length of the major and minor axes of a corresponding ellipse is then extracted by fitting a Principal Component Analysis to the sampled points. The results depend on the resolution and smoothing of the images. For the WMG-generated axons, the lengths of the major and minor axes are extracted directly from the ellipsoids.

For individual axons, mean eccentricity [mean (eccentricity)] is calculated over all measurements along an axon, while the standard deviation of the eccentricity [std (eccentricity)] quantifies the variation of these values.

2.4.3 Tortuosity

We quantify tortuosity based on the tortuosity factor and the maximum deviation. The tortuosity factor is given as the ratio between the geodesic length of the centerline of an axon, and the Euclidean length between the end points of the centerline (see Figure 11). The maximum deviation is given by the maximum Euclidean distance between the centerline and a straight line spanning the endpoints of the centerline when measured perpendicular to the straight line (see Figure 11).

2.5 Phantoms presented here

All phantoms included in this paper, are generated based on the optimization scheme shown in Table 1.

Table 1

Table 1. Optimization scheme used to generate the phantoms presented in the Section 3.

The growth in the first stage is rougher (i.e., lower ellipsoidDensityScaler and higher growSpeed), and then gets finer through the stages (i.e., higher ellipsoidDensityScaler, and lower growSpeed). Different optimization schemes may yield different FVFs and generation times. The following parameters are kept constant and equal for all phantoms:

• targetFVF = 0.8

• mapFromMaxDiameterToMinDiameter =

    {'from' : [0.2, 0.5, 1.25], 'to'  :

[0.2, 0.2, 0.5]}

Some phantoms include cells with CVF = 0.05. The size and shape of the cell clusters are determined by randomly sampling the ellipsoid axes from normal distributions. Each primary axis l₁ is sampled from the distribution l₁ ~ N(μ = 13 μm, σ = 2μm), while the secondary l₂ and tertiary l₃ axes are set to be equal, and sampled from the distribution l₂ = l₃ ~ N(μ = 5 μm, σ = 1 μm). This results in a mean fractional anisotropy of 0.54. The dispersion angles of the cell clusters (i.e., orientation relative to the primary axis of axons) are sampled from a uniform distribution over the interval (−23, 23) deg. This is in accordance with observations of CVF, anisotropy and fiber dispersion in Andersson et al. (2020), where cell clusters from healthy monkey white matter were approximated as ellipsoids (or tensors).

3 Results

To demonstrate the applications of the interactive feature of the WMG, we focus on two scenarios: general cell mobility and inflammatory response. Firstly, we demonstrate the biological relevance of the morphological features expressed in the numerical phantoms generated with the WMG, by comparing the morphology of a set of numerical phantoms with real axons. The real axons originate from a monkey brain, and have been quantified by segmentation from XNH-volumes in previous work (Andersson et al., 2020). Then, by interactively changing the cell configurations over time, we show the morphology of the surrounding fibers changes consequently. We quantify the changes by analyzing longitudinal axon diameter variations, longitudinal variations of cross-sectional eccentricity, and tortuosity of fibers.

3.1 The WMG mimicking different tissue compositions

The global fiber dispersion is controlled by the parameter ϵ. Examples of varying degrees of dispersion are shown both with and without cells in Figure 4. When increasing ϵ, the fibers are forced to bend and deform more to achieve a higher FVF. Similarly, the presence of the static cells forces the fibers to adapt around them.

Figure 4

Figure 4. Demonstration of the global fiber dispersion parameter ϵ and the presence of cells (blue ellipsoids). The unit spheres on the left show the global dispersion associated with each ϵ. A cut is made at 1/3 of the voxel's height to enhance the visualization of individual fiber morphology. Above this height, the meshes are pruned such that only 7% are left. The black voxel marks the boundary of the ellipsoid centers, while the gray voxel marks the volume for which the FVF is optimized. The two voxels are used to avoid boundary effects within the optimized volume. Cross-sections and individual axons are visualized in Supplementary Figure S1.

Crossing fiber bundles can be generated in two configurations—sheets (Kjer et al., 2023) and interwoven. Examples are shown in Figure 5.

Figure 5

Figure 5. Demonstration of two crossing bundle configurations. The black voxel marks the boundary spatial limit of the ellipsoid centers, while the gray voxel marks the volume for which the FVF is optimized. In both examples, the number of fibers is split equally between the two perpendicular bundles. The split and fiber directions can be varied, and the number of crossing bundles can be increased. (A) Sheet configuration where the fibers of each bundle are separated into sheets. (B) Interwoven configuration where the fibers of each bundle are interwoven with each other.

Demyelination can be mimicked by selectively stripping the myelin from axons and adding cells to mimic an inflammatory response as shown in Figure 6. This can be done either between optimization steps by decreasing the allowed diameter range of a fiber to that of the corresponding axon, or by removing the myelin mesh after optimization.

Figure 6

Figure 6. Demonstration of demyelination showing a demyelinated version of the phantom shown in Figure 4 (ϵ = 0.2, with cells). The black voxel marks the boundary for the ellipsoid centers, while the gray voxel marks the volume for which the FVF is optimized.

3.2 Performance of the WMG

Figure 7 shows the obtained FVFs and the processing time for the parameters and optimization scheme described in Table 1. We achieve mean (FVF) ≥ 0.72 for phantoms of all degrees of ϵ without cells and mean (FVF) ≥ 0.66 with cell clusters (CVF = 0.05). It is seen that the processing time increases with ϵ. This is a consequence of the increased fiber dispersion leading to a higher frequency of collisions of the chain of ellipsoids comprising the axons. Each collision necessitates correction, leading to a longer optimization time. Likewise, the increased dispersion complicates dense packing and results in lower FVF. Similarly, the presence of cells (CVF = 0.05) results in a much increased processing time compared to phantoms which have otherwise identical parameter configurations. This is due to the additional volume taken up by the static cell clusters further complicating the packing and resulting in lower FVF.

Figure 7

Figure 7. Performance of the WMG in terms of achieved axonal volume fraction (FVF) and processing time. Each point is based on five samples/repetitions. Mean (FVF) > 0.66 is achieved for all configurations. The complexity of the phantom configuration is increased by increasing dispersion (ϵ), adding cell clusters, and introducing fiber crossings. The higher complexity makes it more challenging to achieve higher FVFs.

3.3 Morphological variation of WMG-generated phantoms

To demonstrate that the morphological features expressed in the WMG-generated phantoms are relevant in relation to real 3D morphology, we compare the morphology of the axons with real axons quantified from XNH-volumes. The comparisons are based on longitudinal axon diameter variations, cross-sectional eccentricity, and tortuosity— all metrics which are reflected in the dMRI signal.

3.3.1 Longitudinal axon diameters

Figure 8 shows an example of target and output distributions of mean axon diameters for individual axons. Much flexibility in morphology parameters is required to obtain the high FVFs. Hence, some deviation from the target is to be expected. For the parameters and optimization scheme described in Table 1, the mean μ values of the distributions are within 0.6 μm of that of the target distribution with μ = 1.8μm.

Figure 8

Figure 8. Target vs. outputted distributions of mean axon diameters for individual axons mean(AD). Gamma distributions are fitted to the histogram bins of the outputted mean axon diameters. The desired diameter of each axon is sampled from a gamma distribution (the desired distribution). Due to the packing procedure, some deviation is expected for the outputted distributions. Here, deviations of mean values are within 0.6 μm.

Figure 9 shows a positive correlation between the std (AD) and mean (AD). We compare with the XNH-samples and see an overlap of the metrics and a likewise positive correlation. We see that std (AD) can be regulated both by varying ϵ and CVF. The higher ϵ forces the individual axons to deform more to pack to the dispersing environment—and vice versa. Similarly for increased CVF (see “ϵ = 0.0, CVF=0.00” vs. “ϵ = 0.0, CVF=0.05”), where the presence of cells force additional ellipsoid deformation in order to obtain the dense packing.

Figure 9

Figure 9. Longitudinal axon diameter variation. Each marker represents one axon from a given sample, while the lines are the linear fits of all axons in a sample/numerical phantom. The standard deviation of the mean axon diameter per axon [std (AD)] correlates positively with the mean axon diameter per axon [mean (AD)]—both for the XNH and WMG axons. Note that the distribution of axon diameters of the XNH axons is not representative of the tissue as a whole, but biased toward the larger axons due to those axons being more clear on the XNH images.

3.3.2 Axonal cross-sectional eccentricity

Figure 10 shows the variation in cross-sectional eccentricity along WMG-generated axons as a function of the mean cross-sectional eccentricity along the axons. Both the mean (eccentricity) and the std (eccentricity) can be regulated by varying ϵ and CVF. Again, because the higher ϵ forces the individual axons to deform more to pack to the dispersing environment—and vice versa. Similarly for increased CVF (see “ϵ = 0.0, CVF=0.00” vs. “ϵ = 0.0, CVF=0.05”), where the presence of cells force additional ellipsoid deformation to obtain the dense packing. The types of eccentricities in the XNH-segmented axons and the WMG-generated axons are different in nature. While the cross sections of the XNH axons are asymmetric and squiggly by nature, the cross sections of the synthetic WMG axons are limited to being elliptic due to the axon being comprised of ellipsoids. Hence, while the eccentricity measure does provide a rough estimation of the dominating cross-sectional eccentricity in both cases, it is not directly comparable across the two.

Figure 10

Figure 10. Cross-sectional eccentricity. We achieve an overlapping degree of eccentricity for the WMG-generated axons as we observed for the XNH axons. Meanwhile, the std (eccentricity) shows a higher degree of variation than the XNH samples. We see a higher degree of eccentricity in the crossing fiber (CF) compared to the corpus callosum (CC) quantified from the XNH images.

3.3.3 Axonal tortuosity

A positive correlation between maximum deviation and tortuosity both for axons from the WMG and XNH is observed. This is shown in Figure 11. The tortuosity can be regulated by the dispersion parameter ϵ and by adding cells. The higher dispersion and CVF mean that the axons have to bend more around each other to not overlap and vice versa.

Figure 11

Figure 11. Tortuosity. We quantify the tortuosity of the axons based on the tortuosity factor and the maximum deviation (both metrics are explained in the figure) to compare with the real axons segmented from XNH-imaging of corpus callosum (CC) and crossing fiber (CF) from a monkey. It is possible to reach overlapping tortuosity factors for the synthetic WMG axons as is observed for the real XNH axons from the corpus callosum (CC) region. However, we do not match the high tortuosity factors and maximum deviations observed for the crossing fiber (CF) region. a, b, c, and d show examples of 2D projections of the axonal centerlines at different ends of the spectrum.

3.4 Interactive configuration allows the mimicking of dynamic environments

Dynamic cells are mimicked by interactively changing the cell configuration during phantom optimization and letting the fibers adapt accordingly. Figure 12 shows how a dynamic cell environment can influence the morphological metrics of the surrounding axons. In the right column, we see how moving cells are causing an increase in axon diameter variations, eccentricity, and tortuosity over time. In Figure 13, we see a similar trend when the CVF is increased over time—except that in this example, the axons are generally straight and have non-eccentric cross-sections at the first time point due to CVF = 0.00.

Figure 12

Figure 12. General cell mobility: The CVF is kept constant while cells are moved around over time. (A) Visualization of three temporal snapshots of a phantom. A cut is made at 1/3 of the voxel's height to enhance the visualization of individual axon morphology. Above this height, the meshes are pruned such that only 7% are left. The black voxel marks the boundary of the ellipsoid centers, while the gray voxel marks the volume for which the FVF is optimized. (B) Morphological metrics computed for the phantom snapshots seen in (A). When starting from otherwise straight axons, the moving cells cause increasing diameter variations (top), eccentricity (middle), and tortuosity (bottom) over time. For an explanation of metrics for eccentricity and tortuosity, see Figures 10, 11, respectively.

Figure 13

Figure 13. Cell increase: the CVF is increased in steps of 0.035 from 0.00 to 0.14 over time. (A) Visualization of three temporal snapshots of a phantom. A cut is made at 1/3 of the voxel's height to enhance the visualization of individual axon morphology. Above this height, the meshes are pruned such that only 7% are left. The black voxel marks the boundary of the ellipsoid centers, while the gray voxel marks the volume for which the FVF is optimized. (B) Morphological metrics computed for the phantom snapshots seen in (A). When starting from otherwise straight axons, the increasing cell fraction causes increasing diameter variations (top), eccentricity (middle), and tortuosity (bottom) over time. For an explanation of metrics for eccentricity and tortuosity, see Figures 10, 11, respectively.

4 Discussion

With the development of the WMG tool, we enable the generation of interactive numerical white matter phantoms composed of the compartments of axons, myelin, cells and extra-cellular space. This allows us to mimic the dynamic nature of white matter by interactively changing the configuration of the cell clusters within a given phantom while the surrounding axon morphology adapts accordingly. With these phantoms, we can thereby analyze how specific cell dynamics [such as mobility and inflammation (Davalos et al., 2005; Nimmerjahn et al., 2005; Tønnesen et al., 2018)] influence the morphology of the surrounding axons. Furthermore, the very same phantoms can be employed as mesh inputs for Monte Carlo diffusion simulations. This enables us to study how these finer detailed morphological changes affect the dMRI signal and to evaluate microstructure models.

When comparing morphological metrics between axons generated with the WMG tool and real axons segmented from XNH-imaging of a vervet monkey brain (corpus callosum and crossing fiber), we see an overlap with the metrics of axon morphologies from the corpus callosum. However, the WMG tool does not offer the degree of tortuosity observed for the crossing fiber region.

Based on the assumption that real white matter morphology carries a history of dynamic events, the WMG tool provides an improved understanding of the dynamics of white matter morphology—both for healthy and diseased tissue.

4.1 WMG-generated axons possess histology-resembling morphological metrics

We analyzed the morphology of the individual axons from phantoms by the WMG tool. Based on the metrics for longitudinal diameter variations (Figure 9), longitudinal variations of cross-sectional eccentricity (Figure 10), and tortuosity (Figure 11), we show that the WMG tool is capable of generating axons with morphology closely resembling that observed for axons from the corpus callosum of a vervet monkey brain with XNH imaging in Andersson et al. (2020, 2022). These morphological metrics have been found to influence the dMRI signal from individual axons (Nilsson et al., 2012; Andersson et al., 2020, 2022; Brabec et al., 2020; Lee et al., 2020; Winther et al., 2023), and are hence essential when generating numerical white matter phantoms with the purpose of evaluating microstructure models for dMRI.

The longitudinal axon diameter variation (Figure 9, y-axis) and the longitudinal variations of cross-sectional eccentricity (Figure 10), and the tortuosity factor (Figure 11, x-axis) reached by the WMG tool, are in agreement with that of the XNH-imaged corpus callosum. The maximum deviations (Figure 11, y-axis), however, are generally much higher in the XNH-imaged corpus callosum. Meanwhile, the XNH-imaged crossing fiber is more challenging to match. Here, only the degree of eccentricity (Figure 10) is in agreement. To obtain a higher degree of tortuosity, and especially the maximum deviation, the ellipsoid chains that make up the foundation of the axons would have to be initialized with an initial tortuosity rather than in a straight line. While this is not a standard configuration of the WMG tool, it can be obtained through manual configuration. No other tools for the numerical synthesis of white matter phantoms are known to output this degree of tortuosity.

Because the low signal-to-noise ratio of the XNH volumes challenged a robust manual segmentation of axons with mean diameters smaller than 2 μm, the segmented axons are not representative of the underlying distribution w.r.t. diameter (Andersson et al., 2020). This quantification does, however, provide valuable insight into the realistic range of morphological metrics for the larger axons.

For extra-cellular space, we targeted a volume fraction of 0.2 (i.e., targetFVF = 0.8) in agreement with what has been reported for normal adult brain tissue-which is between 0.15 and 0.30, and typically 0.20 (Syková and Nicholson, 2008). Figure 7 shows that all WMG phantoms without cells, reach well within the targeted range for all degrees of fiber dispersion. More precisely, we obtain FVFs between 0.83 ± 0.00 and 0.80 ± 0.01 for single fiber bundle phantoms with dispersion between 0 and 18 deg (ϵ of 0.0 and 0.2, respectively), and FVF of 0.72 ± 0.01 for dispersion 90 deg (ϵ of 1.0). Meanwhile, the MEDUSA tool (Ginsburger et al., 2019) reaches FVFs between 0.72 and 0.62 for single-fiber-bundle phantoms with dispersion between 0 and 20 deg, and the ConFiG tool (Callaghan et al., 2020) reaches FVFs between 0.75 and 0.71 for single-fiber-bundle phantoms with dispersion between 7.75 ± 4.23 deg and 17.46 ± 9.92 deg. Although the CACTUS tool (Villarreal-Haro et al., 2023) reaches impressive FVFs between 0.91 and 0.95 for single-fiber-bundle phantoms with dispersion between 0 and 25 deg, such high volume fractions are not necessarily relevant.

4.2 Axon morphology is modulated by cell dynamics

By exploiting the interactive component of the WMG tool to mimic cell dynamics in the form of general cell mobility and cell increase, we show how axon morphology can be modulated over time by these expected dynamics (Figures 12, 13). Under the assumption that real white matter follows similar principles, we can thereby use the WMG tool to study expected tissue response in various physiological or pathological scenarios involving dynamic cell behavior. No other tool offers this feature.

Cell dynamics are continuously occurring in healthy tissue in the form of general glial cell mobility (Davalos et al., 2005; Nimmerjahn et al., 2005; Tønnesen et al., 2018). However, because these changes are homogeneously distributed across the tissue, their influence on morphological metrics should stabilize over time and not cause net changes to the dMRI signal. Here, we mimic cell mobility by changing the configuration of cells over time while keeping the CVF constant (Figure 12). By starting from straight axons, we see how the continuous cell mobility induces an increase in longitudinal morphological variation as the axons adapt to the changing environment.

Cell dynamics are crucial during an inflammatory response (Purves et al., 2001; Davalos et al., 2005; Nimmerjahn et al., 2005; Tønnesen et al., 2018). With the WMG tool, inflammation can be mimicked by the increase of cells (Figure 13). We show how this dynamic behavior induces an increase in longitudinal morphological variation as the axons adapt to the changing environment—similar to the observations for general cell mobility (Figure 12). Inflammation can also involve e.g., edema due to the accumulation of fluid in the extra-axonal space (Nehring et al., 2022). This will in turn cause an expansion of the extra-axonal space and affect the compartmental volume fractions. Such changes were not modeled here.

4.3 Synthesis method and tissue configurations

The high morphological flexibility enabled by the ellipsoidal building blocks of the WMG tool allows for a large variety of different tissue configurations. Here, we show examples of varying dispersion (Figure 4), fiber crossings (Figure 5), demyelination (Figure 6), inclusion of static cells (Figure 4), general cell mobility (Figure 12), and cell increase (Figure 13). A key assumption for the WMG tool is that axon morphology is modulated to the local environment, and thereby adapting according to surrounding axons and cells.

The WMG tool is the first tool developed for generating numerical white matter phantoms which includes an interactive component to allow the mimicking of cell dynamics. However, the tools MEDUSA (Ginsburger et al., 2019), ConFiG (Callaghan et al., 2020), and CACTUS (Villarreal-Haro et al., 2023) do excel when it comes to computational efficiency and more complex morphological features.

Our concept of ellipsoidal building blocks is very similar to that of MEDUSA (Ginsburger et al., 2019) where spherical building blocks are used. While the sphere representation allows for a high representational power of the longitudinal axon morphology characteristics, it does not allow for eccentric cross-sections as documented by histology (Abdollahzadeh et al., 2019; Lee et al., 2019; Andersson et al., 2020). On the other hand, the ellipsoid geometry allows for more degrees of freedom which, in turn, enables eccentric cross-sections—although, still lacking cross-sectional squigglyness to truly resemble histology. Both tools are based on variations of a force-biased packing algorithm, first introduced by Altendorf and Jeulin (2011), where the phantoms are obtained as an equilibrium between “repulsion forces” (for avoiding overlapping axons and cells) and “recover forces” (for ensuring the structure of the individual axons). Similar mechanics are used to shape the axons in CACTUS (Villarreal-Haro et al., 2023), where each axon is represented by capsular building blocks during the optimization of axon trajectories. However, the CACTUS tool contains an additional step of radial optimization, which can increase packing densities and increase the realism of the axonal cross-sections by introducing the desired natural squigglyness. Meanwhile, the ConFiG tool (Callaghan et al., 2020) focuses on the initial growth and maturation of white matter by mimicking how axons are guided by chemical cues and adapt to the available space as they grow and extend from one end to the other. Hence, rather than initializing all axons along straight lines that extend the entire length of the voxel as done in the other tools, each axon is grown stepwise from one end to the other. This is followed by a radial optimization step similar to that of CACTUS (Villarreal-Haro et al., 2023). A similar finishing radial optimization could beneficially be applied in the WMG tool to obtain likewise more realistic cross-sections.

4.4 Limitations and future work

4.4.1 Tissue representativity

Numerical phantoms from the WMG tool can include axons, myelin, and glial cells. These are the most prominent structures in brain white matter; both in regards to volume fractions and signal contrast contributions in dMRI (Jelescu and Budde, 2017)..

However, additional structures present in the brain white matter might affect water diffusion and hence the dMRI signal. Such structures include microtubules, mitochondria, nodes of Ranvier, neuron cell bodies, dendrites, glial cell processes, and morphological changes related to pathology.. These will be implemented in a future version of the WMG tool.

With the WMG tool, the g-ratio is modeled as constant throughout the axon length. However, it has been found that this ratio can vary across myelin segments (Andersson et al., 2020). Demyelination in the WMG tool is mimicked by randomly stripping myelinated axons from their myelin sheets and by increasing CVF. The intermediate steps of the demyelination process are not modeled.

The size of the phantoms is crucial for the representational power of tissue features. For the WMG tool, the processing time is the limiting factor. It has been shown that phantom sizes larger than (200 μm)³ can reduce the sampling bias (Rafael-Patino et al., 2020) and hence improve the representational power. Moreover, while efforts are made to minimize structural boundary effects during synthesis, complete avoidance is challenging.

The WMG tool does offer a way to extend the intra-axonal compartment. All axon and myelin meshes can be mirrored around voxel boundaries without the requirement of further optimization. However, the added amounts of meshing do come with a computational cost when applied in Monte Carlo diffusion simulations—by 3-fold.

4.4.2 Improving the computational efficiency

Creating phantoms is currently a time-consuming process with 28 ± 3 h spent on generating phantoms containing a single fiber direction with ϵ = 0.2 and 55 ± 8 axons in a voxel size of 43.3±1.0 μm using one core (Figure 7). In comparison, for generating a phantom of comparable configuration, the parallelizable CACTUS tool (Villarreal-Haro et al., 2023) requires 4 h for generating a phantom of 33,478 axons in a voxel size of (500 μm)³ using 64 cores. However, since the optimization procedure for the WMG tool relies on checking overlaps between ellipsoids, and this is formulated with linear algebra, there presents an excellent opportunity for leveraging GPUs to enhance efficiency. GPUs are optimized for basic numerical linear algebra operations, and will therefore in many cases outperform CPUs for tasks involving such operations. Therefore, future efforts will focus on implementing the WMG tool in a GPU-compatible manner to accelerate the phantom generation.

While it is already possible to generate large and numerous phantoms with the WMG tool, addressing the computational efficiency will significantly expand the capacity. This enhancement will be crucial for facilitating the generation of large-scale and machine learning-friendly datasets.

4.4.3 Mimicking compression and stress of fibrous tissues

The interactive config-files and the force-biased packing algorithm of the WMG tool could further be applied to mimic the compression and stress of tissue. Such events can be mimicked by applying selective stretch or compression to a voxel and its content through interactive changes to the config-file. Meanwhile, cell elements can be added to mimic the accompanying necrotic debris and cellular immune response. Applications lie within traumatic brain injuries where tissue compression often arises due to external forces, and/or swelling in one area causing compression in adjacent regions (Knight and Kreitzer, 2020).

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding authors.

Author contributions

SW: Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. OP: Conceptualization, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. MA: Conceptualization, Methodology, Resources, Writing – review & editing. HK: Formal analysis, Methodology, Writing – review & editing. JB: Conceptualization, Methodology, Writing – original draft, Writing – review & editing. TD: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. SW has been supported by a DTU Alliance Stipend (project number 102682). SW, MA, and TD have received funding from the European Research Council (ERC) under the European Union's Horizon Europe research and innovation program (grant agreement No. 101044180) (Principal Investigator: TD). OP, MA, and HK have been supported by Capital Region Research Foundation Grant (grant number A5657; principal investigator: TD).

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.2024.1354708/full#supplementary-material

References

Abdollahzadeh, A., Belevich, I., Jokitalo, E., Sierra, A., and Tohka, J. (2021). DeepACSON automated segmentation of white matter in 3D electron microscopy. Commun. Biol. 4:179. doi: 10.1038/s42003-021-01699-w

PubMed Abstract | Crossref Full Text | Google Scholar

Abdollahzadeh, A., Belevich, I., Jokitalo, E., Tohka, J., and Sierra, A. (2019). Automated 3D axonal morphometry of white matter. Sci. Rep. 9:6084. doi: 10.1038/s41598-019-42648-2

PubMed Abstract | Crossref Full Text | Google Scholar

Alexander, D. C., Dyrby, T. B., Nilsson, M., and Zhang, H. (2019). Imaging brain microstructure with diffusion MRI: practicality and applications. NMR Biomed. 32:e3841. doi: 10.1002/nbm.3841

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 Stat. Nonlin. Soft Matter Phys. 83 (4 Pt 1):041804. doi: 10.1103/PhysRevE.83.041804

PubMed Abstract | Crossref Full Text | 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. 117, 33649–33659. doi: 10.1073/pnas.2012533117

PubMed Abstract | Crossref Full Text | Google Scholar

Andersson, M., Pizzolato, M., Kjer, H. M., Skodborg, K. F., Lundell, H., and Dyrby, T. B. (2022). Does powder averaging remove dispersion bias in diffusion MRI diameter estimates within real 3D axonal architectures? Neuroimage 248:118718. doi: 10.1016/j.neuroimage.2021.118718

PubMed Abstract | Crossref Full Text | Google Scholar

Arancibia-Cárcamo, I. L., Ford, M. C., Cossell, L., Ishida, K., Tohyama, K., and Attwell, D. (2017). Node of ranvier length as a potential regulator of myelinated axon conduction speed. Elife 6:e23329. doi: 10.7554/eLife.23329

PubMed Abstract | Crossref Full Text | Google Scholar

Balls, G. T., and Frank, L. R. (2009). A simulation environment for diffusion weighted MR experiments in complex media. Magn. Reson. Med. 62, 771–778. doi: 10.1002/mrm.22033

PubMed Abstract | Crossref Full Text | Google Scholar

Basser, P. J., Mattiello, J., and LeBihan, D. (1994). MR diffusion tensor spectroscopy and imaging. Biophys. J. 66, 259–267. doi: 10.1016/S0006-3495(94)80775-1

PubMed Abstract | Crossref Full Text | Google Scholar

Beaulieu, C. (2002). The basis of anisotropic water diffusion in the nervous system - a technical review. NMR Biomed. 15, 435–455. doi: 10.1002/nbm.782

PubMed Abstract | Crossref Full Text | Google Scholar

Brabec, J., Lasič, S., and Nilsson, M. (2020). Time-dependent diffusion in undulating thin fibers: impact on axon diameter estimation. NMR Biomed. 33:e4187. doi: 10.1002/nbm.4187

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

Davalos, D., Grutzendler, J., Yang, G., Kim, J. V., Zuo, Y., Jung, S., et al. (2005). ATP mediates rapid microglial response to local brain injury in vivo. Nat. Neurosci. 8, 752–758. doi: 10.1038/nn1472

PubMed Abstract | Crossref Full Text | Google Scholar

Dyrby, T. B., Innocenti, G. M., Bech, M., and Lundell, H. (2018). Validation strategies for the interpretation of microstructure imaging using diffusion MRI. Neuroimage 182, 62–79. doi: 10.1016/j.neuroimage.2018.06.049

PubMed Abstract | Crossref Full Text | Google Scholar

Garland, M., and Heckbert, P. S. (1997). “Surface simplification using quadric error metrics,” in Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH '97 (ACM Press/Addison-Wesley Publishing Co.), 209–216. doi: 10.1145/258734.258849

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

Gottschalk, S., Lin, M., and Manocha, D. (1997). OBBTree: a hierarchical structure for rapid interference detection. Comput. Graph 30:237244. doi: 10.1145/237170.237244

Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Hursh, J. B. (1939). Conduction velocity and diameter of nerve fibers. Am. J. Physiol. Legacy Content 127, 131–139. doi: 10.1152/ajplegacy.1939.127.1.131

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

Kjer, H. M., Andersson, M., He, Y., Pacureanu, A., Daducci, A., Pizzolato, M., et al. (2023). Bridging the 3D geometrical organisation of white matter pathways across anatomical length scales and species. eLife 562488. doi: 10.7554/eLife.94917.1

Crossref Full Text | Google Scholar

Knight, W. A., and Kreitzer, N. P. (2020). “Traumatic brain injury,” in Emergency Department Critical Care, eds. J. R. Shiber, and S. D. Weingart (Cham: Springer International Publishing), 393–407.

Google Scholar

Landman, B. A., Farrell, J. A. D., Smith, S. A., Reich, D. S., Calabresi, P. A., and van Zijl, P. C. M. (2010). Complex geometric models of diffusion and relaxation in healthy and damaged white matter. NMR Biomed. 23, 152–162. doi: 10.1002/nbm.1437

PubMed Abstract | Crossref Full Text | Google Scholar

Lee, H.-H., Fieremans, E., and Novikov, D. S. (2021). Realistic microstructure simulator (RMS): Monte carlo simulations of diffusion in three-dimensional cell segmentations of microscopy images. J. Neurosci. Methods 350:109018. doi: 10.1016/j.jneumeth.2020.109018

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. 224, 1469–1488. doi: 10.1007/s00429-019-01844-6

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. Neuroinform. 11:5. doi: 10.3389/fninf.2017.00005

PubMed Abstract | Crossref Full Text | Google Scholar

Nehring, S. M., Tadi, P., and Tenny, S. (2022). Cerebral Edema. StatPearls Publishing. Available online at: https://www.ncbi.nlm.nih.gov/books/NBK537272/

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

Nimmerjahn, A., Kirchhoff, F., and Helmchen, F. (2005). Resting microglial cells are highly dynamic surveillants of brain parenchyma in vivo. Science 308, 1314–1318. doi: 10.1126/science.1110647

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

Purves, D., Augustine, G. J., Fitzpatrick, D., Katz, L. C., LaMantia, A.-S., McNamara, J. O., et al. (2001). Neuroglial Cells. Sinauer Associates. Available online at: https://www.ncbi.nlm.nih.gov/books/NBK10869/

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. Neuroinform. 14:8. doi: 10.3389/fninf.2020.00008

PubMed Abstract | Crossref Full Text | Google Scholar

Sanders, F. K., and Whitteridge, D. (1946). Conduction velocity and myelin thickness in regenerating nerve fibres. J. Physiol. 105, 152–174. doi: 10.1113/jphysiol.1946.sp004160

Crossref Full Text | Google Scholar

Skoven, C. S., Andersson, M., Pizzolato, M., Siebner, H. R., and Dyrby, T. B. (2023). Mapping axon diameters and conduction velocity in the rat brain-different methods tell different stories of the structure-function relationship. arXiv [preprint]. doi: 10.1101/2023.10.20.558833

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

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

PubMed Abstract | Crossref Full Text | Google Scholar

Villarreal-Haro, J. L., Gardier, R., Canales-Rodríguez, E. J., Fischi-Gomez, E., Girard, G., Thiran, J.-P., et al. (2023). CACTUS: a computational framework for generating realistic white matter microstructure substrates. Front. Neuroinform. 17:1208073. doi: 10.3389/fninf.2023.1208073

PubMed Abstract | Crossref Full Text | Google Scholar

Winther, S., Lundell, H., Rafael-Patiño, J., Andersson, M., Thiran, J.-P., and Dyrby, T. B. (2023). Susceptibility-induced internal gradients reveal axon morphology and cause anisotropic effects in the dMRI signal. arXiv [preprint]. doi: 10.1101/2023.05.01.538981

Crossref Full Text | Google Scholar