Hong-Hsi Lee^{1}, Dmitry S Novikov^{1}, and Els Fieremans^{1}

We present a numerical framework for
validation of diffusion models applied to biological tissues. By performing
Monte Carlo simulations in a geometry created by automated segmentation of 3-dimensional
mouse brain electron-microscopy images, we study the time-dependence of the diffusion coefficient and kurtosis along individual white matter axons,
with or without orientation dispersion. Simulation results, together with the analysis
of spatial correlations of axonal cross-section variations along axons, point at the
short-range disorder in restrictions to diffusion, and agree with theoretical
predictions. Our results are consistent with diffusion time-dependence observed in *in vivo*
human brain data.

Diffusion MRI (dMRI)
is sensitive to a length scale of microns, and is an indispensable in vivo technique
for evaluating microstructure in biological tissues. For validating various
diffusional models developed for white matter (WM), previous numerical
simulations were performed either in simple geometries (circles/ellipses/cylinders/spheres/ellipsoids)$$$\,$$$[1-4], or in 2d realistic micro-geometry reconstructed from light or electron
microscopy$$$\,$$$[5-6]. However, simulations in a realistic *three-dimensional microstructure*, that reflects genuine tissue
properties, have been missing.

Here we present a framework for segmenting axons from sequential-slice election microscopy (SSEM) images in the corpus callosum (CC) of a mouse brain. By performing Monte Carlo (MC) simulations in the obtained realistic 3d intra-axonal space (IAS), we study the time-dependence of diffusion coefficient and kurtosis along axons with or without axonal orientation dispersion. Having full control of tissue microstructure allows us for the first time to validate theory and previously reported human brain data$$$\,$$$[7-8] .

$$$\bf{SSEM}$$$

A female 8-week-old C57BL/6 mouse was sacrificed and fixed by transcardiac-perfusion. The genu of CC was excised and analyzed with a scanning electron microscope (Zeiss Gemini-300), which acquired 401 consecutive images of$$$\,6000\times8000\,$$$pixels. Part of the data was discarded due to distortion, leading to a volume (Fig.1a) of$$$\,36\times48\times20\,\mu$$$m$$$^3\,$$$with $$$6\times6\times100\,$$$nm$$$^3$$$-resolution.

$$$\bf{IAS\,Segmentation}$$$

Data was downsampled to$$$\,24\times24\times100\,$$$nm$$$^3$$$-resolution before segmentation. To
segment long axons passing through all slices, we employed diffusion: We manually
seeded an initial position per axon within the central slice (451 seeds),
and filled the IAS by cubic-lattice-hopping diffusion trajectories of 1000
particles per axon, confined by a myelin mask obtained by
expectation-maximization algorithm. Segmented axons with imperfect myelin mask were
deleted (~135) by proofreading, resulting in an IAS mask of **216 fully segmented long axons**
($$$>20\,\mu$$$m in length, discarded ~100 short ones).
The IAS mask was downsampled into$$$\,(0.1\,\mu$$$m$$$)^3$$$-resolution, and transformed into 3d polyhedra by
tetrahedral-filling (Fig.1b). The effect of orientation dispersion was
controlled by subsequent realigning axons along the z-axis (Fig.1c).
To determine the structural universality class of
geometries along axons, we calculated the power-spectrum$$$\,\Gamma(k)=\tilde{\rho}(-k)\tilde{\rho}(k)/L$$$,
where$$$\,\tilde{\rho}(k)\,$$$is the Fourier transform of$$$\,\rho(z)$$$, which is
the normalized IAS cross section along axons (Fig.2b).

$$$\bf{Numerical\,Simulation}$$$

MC simulations of random walkers were implemented in Matlab in continuous space bounded by the boundary of polyhedral IAS, and validated in a polyhedron-approximated cylinder (Fig.3). Next, simulations were performed inside 216 IAS polyhedra (Fig.1c). $$$1.1\times10^5$$$ random walkers were distributed in proportion to axonal volume. Each particle diffused over $$$1.25\times10^5$$$ steps with a duration$$$\,8\times10^{-3}\,$$$ms and a step length$$$\,0.098\,\mu$$$m. Maximal diffusion time$$$\,t=100\,$$$ms. Axo-plasmic diffusivity $$$D_0=2\,\mu$$$m$$$^2$$$/ms. Total calculation time ~35hours on 200 CPU cores of a high-performance-computing cluster.

Diffusivity and kurtosis were estimated by cumulants ($$$\langle\,x^2\rangle$$$,$$$\,\langle\,x^4\rangle$$$). To simulate dispersion, we calculated total diffusivity and kurtosis when axons oriented on a cone with a polar angle$$$\,\theta=[0^\circ,15^\circ,30^\circ,45^\circ]\,$$$(no dispersion to high dispersion). We then compared simulation results with theory (described below) and previously published human studies using mono-polar PGSE$$$\,$$$[8] and STEAM$$$\,$$$[7] (Table 1b-c).

$$$\bf{Theory}$$$

The non-Gaussian diffusion along WM tracts at long
diffusion time$$$\,t$$$ probes the geometry and length scale of restrictions
along axons$$$\,$$$[9]. Randomly-positioned restrictions along axons (short-range
disorder in 1d) yield the following functional
form$$$\,$$$[9] for the axial diffusivity
$$D^\parallel(t)\simeq\,D_\infty+c\cdot\,t^{-\vartheta}\,,\quad\quad(1)$$
and axial kurtosis$$$\,$$$[10]:
$$K^\parallel(t)\simeq\,K_\infty+\frac{2c}{D_\infty}\cdot\,t^{-\vartheta}\,,\quad\quad(2)$$
characterized by the dynamical exponent $$$\vartheta=1/2$$$.
Here$$$\,D_\infty\,$$$is the bulk diffusivity
at$$$\,t\to\infty$$$,$$$\,c\,$$$is the strength of restrictions, and$$$\,K_\infty\,$$$is the kurtosis at $$$t\to\infty$$$, resulting from variances of$$$\,D_\infty^{(i)}\,$$$between different axons. **The time-dependence of$$$\,K^\parallel(t)\,$$$is completely determined by the effective-medium parameters** ($$$D_\infty,\,c$$$). We tested this
prediction using simulations.

Fig.2a shows the diameter histogram of long axons, and Fig.2b shows the 1-dimensional power spectrum along axons, where the red-dashed line indicates that$$$\,\Gamma(k)\propto\,k^p\,$$$approaches a plateau at$$$\,k\to0\,$$$(structural exponent$$$\,p=0\,$$$equivalent to short-range disorder), predicting the dynamical-exponent$$$\,\vartheta=(p+d)/2=1/2\,$$$in dimension$$$\,d=1$$$, Eqs.(1-2)$$$\,$$$[9].

Fig.3 validates our MC simulations in a polyhedron-approximated cylinder, reproducing known short-$$$t\,$$$[11] and long-$$$t\,$$$[12] theoretical results.

Fig.4a-d shows that, in realistic IAS polyhedra, simulated$$$\,D^\parallel(t)\,$$$is consistent with the time-dependence in Eq.(1) (fitted red-dashed line) even for highly dispersed axons ($$$\theta=45^\circ$$$). Table 1a shows corresponding fitted parameters ($$$D_\infty$$$,$$$\,c$$$), consistent with the human data in Table 1b.

To confirm that the micro-geometry along axons
falls into the short-range-disorder universality class$$$\,$$$[9], we predicted the
time-dependence of simulated $$$K^\parallel(t)$$$ based on effective medium
parameters from$$$\,D^\parallel(t)\,$$$in Table 1a. The theoretical prediction (slope$$$\,2c/D_\infty\,$$$of red-dashed line with respect to$$$\,1/\sqrt{t}$$$), based
on Eq.(2), is shown in Fig.4e-h, agreeing with the time-dependence at long$$$\,t\,$$$from MC, when dispersion angle$$$\,\theta\leq30^\circ$$$.

Simulation results are consistent with in vivo measurements and theoretical predictions: diffusion along axons, either with or without dispersion, is characterized by 1d short-range disorder and the$$$\,\vartheta=1/2\,$$$dynamical exponent, Eqs(1)-(2).

The proposed framework of axon segmentation and MC simulation in 3d realistic SSEM microstructure suggests a new validation method that could be extended for various other MRI contrast mechanisms, e.g. magnetization transfer, susceptibility, relaxation, exchange, etc.

[1] Li, J. R., et al., JMR, 248, 54-65 (2014)

[2] Harkins, K. D. & Does, M. D., Phys Med & Bio, 61, 4729-4745 (2016)

[3] Lin, M., et al., MRM, 76, 290-300 (2016)

[4] Palombo, et al., Proc ISMRM, 25, 1087 (2017)

[5] Chin, C. L., et al., MRM, 47, 455-460 (2002)

[6] Xu, T., et al., MRM (2017)

[7] Fieremans, E., et al., NeuroImage, 129, 414-427 (2016)

[8] Papaioannou, A., et al., Proc ISMRM, 25, 0723 (2017)

[9] Novikov, D. S., et al., PNAS, 111(14), 5088-5093 (2014)

[10] Novikov, D. S. & Kiselev, V. G., NMR Biomed, 23, 682-697 (2010)

[11] Mitra, P. P., Phys Rev B, 47(14), 8565-8574 (1993)

[12] Callaghan, P. T., JMR A, 113, 53-59 (1995)

[13] Burcaw, L. M., et al., NeuroImage, 114, 18-37 (2015)

[14] Lee, H. H., et al., Proc ISMRM, 25, 0839 (2017)

Figure 1. (a) The EM image with segmented IAS of all
axons passing through the central slice. (b) 3d polyhedrons of IAS generated
from the axon segmentations in (a) (from a different view angle). (c) Examples
of axon polyhedrons aligned to the z-axis. Only axons passing through all
slices (> 20 micron) were chosen, yielding 216 long axons in total. Each
axon is then cropped into 18 micron in length.

Figure 2. (a) Diameter
histogram of the 216 long axons. (b) Power spectrum $$$\Gamma(k)$$$ of
normalized cross sections along axons shows a plateau at $$$k\to0$$$ (red
dashed line), indicating a $$$\vartheta=1/2$$$ dynamical exponent, corresponding
to short-range disorder in 1d (randomly positioned restrictions along axons),
leading to the diffusion time-dependence in Eqs.(1-2).

Figure 3. Simulation
results (data points) in a cylinder approximated by a polyhedron of radius = 2 micron.
Dashed lines are theoretical predictions. $$$t_c=r^2/D_0=2$$$ ms (a) At short
time, axial diffusivity (AD) (blue) shows no time-dependence, and radial
diffusivity (RD) (red) shows a time-dependence in surface-to-volume limit [11]. (b) At long
time, AD shows no time-dependence, and RD shows a $$$1/t$$$ dependence. (c) At
short time, axial kurtosis (AK) (blue) is zero, and radial kurtosis (RK) (red)
shows a time-dependence derived from the exact solution of the signal [12].
(d) At long time, AK=0, and RK=-1/2, consistent with theory [13].

Figure 4. Simulation
results (data points) in 216 3d realistic IAS polyhedra shown in Fig.1c.
Each axon's axial direction lies on a cone with a dispersion polar angle
$$$\theta=\,$$$(a)$$$\,0^\circ$$$ (no dispersion), (b)$$$\,15^\circ$$$,
(c)$$$\,30^\circ$$$, and (d)$$$\,45^\circ\,$$$(high dispersion), for
$$$D^\parallel(t)$$$ time-dependence. The red-dashed lines in (a-d) are fits to
Eq.(1), with parameters shown in Table 1a. To further validate the
time-dependence along axons (short-range disorder in 1d), we predict the
$$$K^\parallel(t)$$$ time-dependence at (clinically) long time
($$$t=20-100\,$$$ms) with Eq.(2) and parameters in Table 1a with
$$$\theta=\,$$$(a)$$$\,0^\circ\,$$$(no dispersion), (b)$$$\,15^\circ$$$, (c)$$$\,30^\circ$$$,
and (d)$$$\,45^\circ\,$$$(high dispersion). The slope of red-dashed lines in
(e-h) are predictions, with no tunable parameters.

Table 1. Fitting parameters of $$$D^\parallel(t)$$$ in (a) our
simulations in 3d realistic IAS polyhedrons reconstructed from EM of mouse
brain genu of CC, (b) human brain data measured using monopolar PGSE [8], and
(c) human brain data measured using monopolar STEAM [7]. Simulation results for
dispersion angle$$$\,\theta=15^\circ$$$ are consistent to the human brain PGSE
data in genu. Results in human brain STEAM data have a much larger$$$\,c$$$
than other data, since potentially STEAM is sensitive to water exchange between
myelin water and intra-/extra-axonal water [14], magnifying changes of
$$$D^\parallel(t)$$$ with $$$t$$$. (ACR/SCR/PCR=anterior/superior/posterior
corona radiate, ALIC/PLIC=anterior/posterior limb of internal capsule,
genu/midbody/splenium of CC)