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

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

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

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

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 c=12(pi,j-1+pi,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.

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

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 x(r)=SSTr||STr||. 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.

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)}(r)=rT(ei1,j1(r)-ei2,j2(-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 minro(i1,j1),(i2,j2)(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.

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 12o(i1,j1),(i2,j2)(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 s1=-δ(||xi1,j1(r)||)12o(i1,j1),(i2,j2)(r)||xi1,j1(r)||=-δ(||xi1,j1(r)||)o(i1,j1),(i2,j2)(r)2||xi1,j1(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: s2=μi-||xi1,j1(r)||||xi1,j1(r)||=μi||xi1,j1(r)||-1. So the resulting value of s is s = max(s₁, s₂). So the deformation matrix D₁ for ellipsoid j₁ is D1=max(s1,s2)rrT+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 pil,jl+(-1)l12o(i1,j1),(i2,j2)(r)r for l ∈ {1, 2}.

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

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.

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.

We quantify cross-sectional eccentricity based on elliptic parameterization by e=1-b2/a2 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.

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

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

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.

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.

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.

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