Simulation of Brownian motion of water molecules in WM tissue plays an important role in validating models that describe the link between the measured diffusion-weighted MR signal and the underlying tissue microstructure.^1-5 In such a simulation framework, it is crucial that specific tissue characteristics mimic those observed in real tissue as closely as possible. Failing to do so may induce a considerable bias that could confound subsequent analyses.

In this work, we present a novel approach for simulating the axon packing in nerve bundles, which optimizes tissue features in terms of fiber configuration. We demonstrate that in comparison with the conventional techniques,^3,6,7 our proposed method achieves higher packing densities and eliminates the bias in the radii distribution, while preserving the short-range disorder of the axon packing structure.

The distribution of axon radii within a nerve bundle, can be approximated with a gamma distribution.^8,9 In our proposed packing algorithm, axons are represented by circles with radii generated from a gamma distribution and are assumed to be confined in a square with width $$$W$$$. The circles do not overlap and are separated by interfaces with specific boundary thickness, $$$h=\eta r$$$, where $$$r$$$ is the axon radius and $$$\eta$$$ is the separation parameter.

Given a desired set of model parameters, the algorithm starts with generating a circle in the center of the square continuing in levels around the first circle. At each step, a new circle is placed such that its boundary is tangent to the boundaries of the two previously placed neighboring circles, i.e., $$\begin{cases}(x_i-x_{L_1})^2+(y_i-y_{L_1})^2=(r_{L_1}+h_{L_1}+r_i+h_i)^2\\(x_i-x_{L_2})^2+(y_i-y_{L_2})^2=(r_{L_2}+h_{L_2}+r_i+h_i)^2\end{cases},$$ where $$$(x_i,y_i)$$$ are the coordinates of the center of the $$$i$$$-th circle with radius $$$r_i$$$ and boundary thickness $$$h_i$$$, and $$$(x_{L_1},y_{L_1})$$$, $$$(x_{L_2},y_{L_2})$$$ are the coordinates of the center points of the first and second neighboring circles as depicted in Fig.1. The algorithm terminates when the number of placed circles reaches $$$N_A$$$. For a given packing density ($$$\phi$$$), $$$N_A$$$ can be calculated as $$N_A=\frac{W^2\phi}{(1+\eta)^{2}v_0},$$ where $$$v_0$$$ is the mean axon area. For radii generated from a gamma distribution, $$$Gamma(\alpha,\beta)$$$, where $$$\alpha$$$ and $$$\beta$$$ denote the shape and scale parameters, $$$v_0$$$ can be calculated as $$v_0=\pi\frac{\Gamma(\alpha+2)\beta^2}{\Gamma(\alpha)},$$ where $$$\Gamma(\alpha)$$$ is the gamma function. The flowchart of the algorithm is presented in Fig.2.

Monte-Carlo simulations (1000 repetitions) with model parameter settings as shown in Fig.3(a) were performed to evaluate performance in terms of bias in estimated mean radius ($$$\mu_r$$$), radius variance ($$$\sigma_r^2$$$) and packing density ($$$\phi$$$). In addition, the structure correlation function (SCF) in Fourier domain is used to describe the level of randomness of the axon packing results, i.e., $$C(k)=\iint\,e^{-i\mathbf{kr}}\rho(\mathbf{r}+\mathbf{r}_0)\rho(\mathbf{r}_0)d\mathbf{r}\frac{d\mathbf{r}_0}{W^{2}},$$ where $$$\rho(\mathbf{r})$$$ is the density of axons at point $$$\mathbf{r}$$$. $$$\rho(\mathbf{r})$$$ is 1 for points inside axons and 0 for outside. $$$\phi v_0$$$ and $$$\langle r_{ext}\rangle=\sqrt{v_0/\pi}$$$ (external object radius) are used for the SCF and $$$k$$$-axis normalization respectively.^10,11

Fig.4 shows the outcome of both techniques and the lattice packing. The related estimated radii PDF is presented in Fig.3(b), showing a bias in the estimated mean radius. The analysis of estimated mean radius and variance, depicted in Fig.5(a-b), illustrates that: (i) the conventional technique underestimates the model parameters and (ii) the uncertainty of the estimations is smaller for the proposed technique. The bias in the mean radius and variance for the conventional method is due to limited free space for new axons and rejection of larger axons as the procedure evolves.

The analysis of simulated packing densities presented in Fig.5(c) demonstrates that for a given separation parameter, the proposed method results in higher packing densities. However, the analysis of SCF, Fig.5(d), shows that the plateau in $$$C(k)$$$ as $$$k\to0$$$ for the proposed method is smaller than for the conventional (indicating stronger spatial correlation), but the method still preserves short-range disorder property (distinctive from completely ordered packing, lattice in Fig.4(c)).

A novel approach was proposed for simulating the placement of axons in nerve bundles. Comparison of the results to the conventional method demonstrated that our method (i) eliminates the bias in estimated radius mean and variance (ii) results in tissue models which are highly consistent with the desired characteristics of nerve bundles, (iii) can result in higher packing densities for a given set of model parameters, and (iv) while is more correlated at short scales, still preserves the short-range disorder property of the tissue. However, the exact degree of randomness in WM tissue and how nature packs the axons remains an open question.

The resultant axon packings can be used subsequently to study the Brownian motion of water molecules within nerve bundles and the effect of axon configuration properties on the measured diffusion MR signal.