Simulation of diffusion in axons with harmonic and stochastic trajectories
Jan Brabec1, Samo Lasic2, and Markus Nilsson1

1Clinical Sciences Lund, Radiology, Lund University, Lund, Sweden, 2Random Walk Imaging, Lund, Sweden


We investigated effects on the diffusion MRI signal from harmonic and stochastic variations in the trajectories of thin axon-like structures. Trajectories with different types of variations exhibited characteristically different diffusion spectra (time-dependence) that, for low frequencies, bare similarities with those from restricted diffusion environments. Non-straight trajectories can thus bias axon diameter estimation. Results also indicated that observable effects of structural disorder may not be specific for extracellular diffusion but may also be found for intra-axonal spaces.


Time-dependent diffusion allows the probing of morphological properties of biological tissues.1 The key question is which properties that are preserved in the dMRI signal. Such questions are often answered by using simulations that assume axons to be straight cylinders2 or, for undulating axons, periodic sinusoidal waves.3 However, microscopy of human white matter shows the presence of stochastic variations of axonal trajectories.4 Here, we investigated the impact of both harmonic and stochastic trajectories of axon-like structures on the time-dependence of diffusion. Trajectories were characterized by descriptive parameters, and related to properties of the diffusion spectra. Results highlight impacts of an often-neglected aspect of white matter microstructure: non-straight axonal trajectories.


Generation of axonal trajectories: To generate a toy model that mimics biologically plausible axonal trajectories, axons were represented as thin wires with a position given by

$${\rm y(\mathit{x} )=a\cdot sinā”(\mathit{\phi})=a\cdot sin(2\pi\cdot \mathit{x}/\lambda+r(\mathit{x}))}$$

with amplitudes of $$${\rm a = [2;10]\, \mu m}$$$, wavelengths $$${\rm \lambda=[25;150] \, \mu m}$$$ and a stochastic part, $$${\rm r(\mathit{x})}$$$, generated by cumulative sum of random numbers with variable degree of correlation. For harmonic sine waves $$${\rm r(\mathit{x})=0}$$$. The intrinsic diffusivity within the impermeable paths was set to $$${\rm D_0=1.7\, \mu m^2/ms}$$$.

Quantification strategies: Three parameters describing the axonal trajectories were defined from the trajectory, $$${\rm y(\mathit{x})}$$$. First, we defined the microscopic orientation dispersion, $$${\rm \mu OD}$$$, as:

$${\rm \mu OD=\bigg\langle\left(tan^{-1}\left(\frac{d}{d\mathit{x}}y(\mathit{x})\right)\right)^2\bigg\rangle }$$

Second, the dispersion-weighted wavelength was defined by

$${\rm \lambda_{\sigma}=\Bigg(\frac{\langle \mu OD(\mathit{x})\cdot\lambda^{-1}(\mathit{x})\rangle }{\langle\mu OD(\mathit{x})\rangle}\Bigg)^{-1} }$$

where the ‘local’ wavenumber is $$${\rm \lambda^{-1} (\mathit{x})=d/d\mathit{x} \, \phi(\mathit{x}) }$$$. Third, effects of correlation lengths, $$${\rm l_d}$$$, were also studied (Figure 3). This was defined as the decay rate of the envelope of square of autocorrelation function of $$${\rm y(\mathit{x})}$$$.

Simulating diffusion spectra: The diffusion spectra, $$${\rm D(\omega)}$$$, is the spectrum of the velocity auto-correlation function5,

$${\rm \chi (\mathit{t})= \langle v(\mathit{t})v(0) \rangle}$$

And it has a close relationship with the mean square displacement, according to

$$\langle \Delta x^2 (\mathit{t})\rangle=\int_0^{\mathit{t'}}\int_0^{\mathit{t''}}\chi(\mathit{t'})\,d\mathit{t'}d\mathit{t''}=\frac{4}{\pi}\int_0^{\infty}\frac{D(\mathit{\omega})}{\mathit{\omega^2}}[1-cos(\mathit{\omega t})]\,d\mathit{\omega} $$

Diffusion spectra were numerically obtained by Gaussian sampling of the mean square displacement perpendicularly to the direction of axonal propagation (Figure 2 and 3). Gaussian sampling relies on simulating and averaging over different diffusion particle trajectories weighted by a Gaussian function.

Characterizing diffusion spectra: The obtained diffusion spectra were characterized by three parameters: spectral width, spectral height, and low frequency power law behavior (Figure 1A). The height was estimated at high frequencies where the diffusion spectrum remained constant, the width was estimated as the frequency at the half of the height, and the power law behavior at low frequencies was estimated by fitting the spectra in a loglog plot constrained up to half of the width.The descriptive parameters were also predicted from the trajectories themselves. The spectral width was predicted as

$${\rm\omega_{\, predicted}=\bigg(\frac{\langle\mu OD(\mathit{x})\cdot\lambda^{-1}(\mathit{x})\rangle}{\langle\mu OD(\mathit{x}) \rangle}\bigg)^{-2}\cdot D_0}$$

Spectral height was predicted as

$${\rm D_{\infty; \, predicted}=\langle\mu OD(\mathit{x})\rangle}\cdot D_0$$

Power law behavior was not predicted, but we note that the diffusion spectra at low frequencies scales as $$${\rm \omega^p}$$$, with $$${\rm p=2}$$$ for restricted environments and $$${\rm p<2}$$$ for disordered systems.6

Results and Discussion

Characteristics of the trajectories had a clear effect on the diffusion spectra (Fig. 2 and 3). Interestingly, varying the level of randomness in the stochastic trajectories had little impact on the overall shape of the spectra (Fig. 3). The estimated width and height of the diffusion spectra agreed well with the values predicted from the trajectories (Fig. 4). The height reflects that, at short diffusion times, the random walks of diffusing particles do not sense the curvature of the trajectory but behave as if they diffused in straight but orientationally dispersed subsegments. Moreover, the width was determined by the average wavelength of the undulating trajectories, weighted by the local orientation dispersion. The power law exponent was found to vary between 1 and 2 for different trajectories, but no substantial difference was found between harmonic and stochastic undulations (Fig 4C).


Relevant characteristics of the diffusion spectra can be predicted by a small number of structural properties of the trajectories. These properties do not describe well-defined aspects of the trajectories, such as their amplitude and wavelengths, but rather parameters such as microscopic orientation dispersion and a dispersion-weighted wavelength. Custom gradient waveforms, and in particular short diffusion times, may be used to estimate such features7 (Fig. 5). Finally, we found a type of power law behavior in these axon-like structures otherwise assigned specifically to the extra-axonal space. This indicates that several emerging methods for axon diameter estimation8 may be confounded by axonal trajectories.


No acknowledgement found.


1. Novikov, Dmitry S., et al. "Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation." arXiv preprint arXiv:1612.02059 (2016).

2. Hall, Matt G., and Daniel C. Alexander. "Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI." IEEE transactions on medical imaging 28.9 (2009): 1354-1364.

3. Nilsson, Markus, et al. "The importance of axonal undulation in diffusion MR measurements: a Monte Carlo simulation study." NMR in Biomedicine 25.5 (2012): 795-805.

4. Schilling, Kurt G., et al. "Histological validation of diffusion MRI fiber orientation distributions and dispersion." NeuroImage (2017).

5. Stepisnik, Janez. "Time-dependent self-diffusion by NMR spin-echo." Physica B: Condensed Matter 183.4 (1993): 343-350.

6. Novikov, Dmitry S., et al. "Revealing mesoscopic structural universality with diffusion." Proceedings of the National Academy of Sciences 111.14 (2014): 5088-5093.

7. Aslund, Ingrid, et al. "Diffusion NMR for determining the homogeneous length-scale in lamellar phases." The Journal of Physical Chemistry B 112.10 (2008): 2782-2794.

8. Burcaw, Lauren M., Els Fieremans, and Dmitry S. Novikov. "Mesoscopic structure of neuronal tracts from time-dependent diffusion." NeuroImage 114 (2015): 18-37.


Figure 1: Parameters used to characterize the diffusion spectra (A) and the trajectories (B). A: The spectral width ($$${\rm \omega}$$$); spectral height at high frequencies ($$${\rm D_{\infty}}$$$) that represents diffusivity at short diffusion times, and power law behavior at low frequencies (power $$${\rm p}$$$ in $$${\rm \omega^p}$$$). B: Illustration of parameters used to characterize the stochastic trajectories.

Figure 2: Diffusion spectra in harmonic sine trajectories, with various wavelengths (columns) and amplitudes (rows). Varying the amplitude was reflected in the height of the spectra ($$${\rm D_{\infty}}$$$). Varying wavelengths impacted both the height and the width ($$${\rm D_{\infty}}$$$ and $$${\rm \omega}$$$, respectively).

Figure 3. Diffusion spectra in stochastic trajectories. Columns show variable microscopic orientation dispersions, and rows show variable correlation lengths. Shorter correlation lengths indicate a more disordered trajectory, but contrary to our intuition, it did not lead to significant changes in the width of the diffusion spectra. However, disorder did change the microscopic orientation dispersion leading to higher diffusivities at short diffusion times ($$$\rm{D_{\infty}}$$$). Importantly, the spectra from the different trajectories within each column are indistinguishable, meaning that diffusion experiments cannot be used to tell them apart.

Figure 4. Comparison between predicted and estimated parameters of the diffusion spectra. A: spectral width ($$${\rm\omega}$$$), B: spectral height ($$${\rm D_{\infty}}$$$) and C: power law behavior ($$${\rm p}$$$). Black points correspond to harmonic trajectories, brown to stochastic. Overall, the estimated parameters agree with the predicted ones. The power behavior diffusion spectra at low frequencies showed rather a distribution of power coefficients, not a single power.

Figure 5: How to design an experiment. When using gradients that probe long diffusion times, observable differences, being the inner product between dephasing and diffusion spectra, becomes almost indistinguishable between cylinders, stochastic and harmonic trajectories. However, gradients that probe rather short diffusion times could be used to distinguish between diverse types of restriction.

Proc. Intl. Soc. Mag. Reson. Med. 26 (2018)