# Realistic Microstructure Simulator (RMS): Monte Carlo simulations of
diffusion in three-dimensional cell segmentations of microscopy images

Hong-Hsi Lee  Els Fieremans  Dmitry S. Novikov  Center for Biomedical Imaging,
Department of Radiology, New York University School of Medicine, New York, NY,
USA  Center for Advanced Imaging Innovation and Research (CAI2R), New York
University School of Medicine, New York, NY, USA

###### Abstract

Background: Monte Carlo simulations of diffusion are commonly used as a model
validation tool as they are especially suitable for generating the diffusion
MRI signal in complicated tissue microgeometries.

New method: Here we describe the details of implementing Monte Carlo
simulations in three-dimensional (3d) voxelized segmentations of cells in
microscopy images. Using the concept of the corner reflector, we largely
reduce the computational load of simulating diffusion within and exchange
between multiple cells. Precision is further achieved by GPU-based parallel
computations.

Results: Our simulation of diffusion in white matter axons segmented from a
mouse brain demonstrates its value in validating biophysical models.
Furthermore, we provide the theoretical background for implementing a
discretized diffusion process, and consider the finite-step effects of the
particle-membrane reflection and permeation events, needed for efficient
simulation of interactions with irregular boundaries, spatially variable
diffusion coefficient, and exchange.

Comparison with existing methods: To our knowledge, this is the first Monte
Carlo pipeline for MR signal simulations in a substrate composed of numerous
realistic cells, accounting for their permeable and irregularly-shaped
membranes.

Conclusions: The proposed RMS pipeline makes it possible to achieve fast and
accurate simulations of diffusion in realistic tissue microgeometry, as well
as the interplay with other MR contrasts. Presently, RMS focuses on
simulations of diffusion, exchange, and $T_{1}$ and $T_{2}$ NMR relaxation in
static tissues, with a possibility to straightforwardly account for
susceptibility-induced $T_{2}^{*}$ effects and flow.

###### keywords:

Monte Carlo simulation , Realistic microstructure , Numerical validation ,
Diffusion MR , Biophysical modeling

††journal: Journal of Neuroscience Methods

##  1 Introduction

The MRI measurements of self-diffusion of water molecules in biological
tissues provide the sensitivity to the diffusion length scales ranging from
microns to tens of microns at clinically feasible diffusion times. As the
feasible range of diffusion lengths is commensurate with the sizes of cells,
diffusion MRI allows one to evaluate pathological changes in tissue
microstructure  in vivo. To balance between accuracy and precision in
estimation of tissue parameters through biophysical modeling of diffusion MR
signal, assumptions are inevitably made to simplify tissue microgeometry
(Grebenkov, 2007; Jones, 2010; Kiselev, 2017; Jelescu and Budde, 2017;
Alexander et al., 2018; Novikov et al., 2019). It is necessary to validate the
assumptions of models before use, either through experiments in physical
phantoms (Fieremans and Lee, 2018), or testing the model functional forms in
animals and human subjects (Novikov et al., 2018), or numerical simulations
(Fieremans and Lee, 2018).

So far, numerical simulation is the most flexible and economic choice among
all kinds of validation. Benefiting from the recent advances in microscopy,
realistic cell geometries for simulations have been directly reconstructed
from the microscopy data of neuronal tissues in 2 dimension (2d) (Chin et al.,
2002; Xu et al., 2018) and 3d (Nguyen et al., 2018; Palombo et al., 2019; Lee
et al., 2020b, a), as shown in Fig. 1. The emerging need for simulations in
realistic substrates prompts the development of open-source software congenial
to physicists, biologists and clinicians.

Here, we describe our implementation of Monte Carlo (MC)-based diffusion
simulations: the Realistic Microstructure Simulator (RMS), which entails a
fast and accurate model validation pipeline. While our pipeline has been
recently announced and applied to simulate diffusion MRI in axonal
microstructure (Lee et al., 2020b, a, d), these publications are mainly
focused on the physics of diffusion and model validation. In this work, we
describe the methodology in detail, building on the algorithms introduced by
our team over the past decade (Fieremans et al., 2008, 2010; Novikov et al.,
2011, 2014; Burcaw et al., 2015; Fieremans and Lee, 2018; Lee et al., 2020c),
and in particular, derive the finite MC-step effects relevant for the
interactions (reflection and permeation) of random walkers and membranes.

Fig. 1: Benefiting from advances in microscopy techniques and from increases
in computational power, the research focus of diffusion simulations has
gradually switched to 3-dimensional tissue microgeometries. Realistic white
matter axons in (Lee et al., 2020b) are examples of voxelized simulation
substrates used in RMS. The figure is adapted from (Chin et al., 2002; Xu et
al., 2018; Nguyen et al., 2018; Palombo et al., 2019; Lee et al., 2020b) with
permission from Wiley, Elsevier, and Springer Nature.

RMS is introduced as follows. In Section 2, we provide an overview of RMS
implementation. Theoretical results and implementation details of particle-
membrane collisions and exchange are presented in A and B; the first order
correction of membrane permeability due to a discretized diffusion process is
derived in C. In Section 3, we demonstrate the application of RMS to diffusion
simulations in realistic axonal shapes reconstructed from electron microscopy
data of a mouse brain (Lee et al., 2019). Simulated diffusion MR signals are
shown to be closely related to features of cell shape, facilitating the
interpretation of diffusion measurements in biological tissues. Finally, in
Section 4 we provide an outlook for microstructure simulation tools in
general, and RMS in particular, as a platform for MR-relevant simulations of
diffusion and relaxation in microscopy-based realistic geometries.

##  2 RMS Implementation

###  2.1 Realistic Microstructure Simulator: An overview

The goal of our RMS implementation is to provide a universal platform of MC
simulations of diffusion in any realistic microgeometry based on microscopy
data. Therefore, the RMS has the following properties:

  1. 1.

The simulation is performed in 3d continuous space with voxelized
microgeometry. We will introduce this main feature of RMS in Section 2.1.1.

  2. 2.

The particle-membrane interaction of impermeable membrane is modeled as
classic elastic collision. The reason for this choice will be explained in
Section 2.1.2, A and B.

  3. 3.

The models of permeable and absorbing membranes are implemented based on the
novel theoretical results tailored for diffusion simulations of high accuracy,
with details in Section 2.1.3 and C.

  4. 4.

The interplay of diffusion and other MR contrasts, such as $T_{1}$ and $T_{2}$
relaxation, surface relaxation and magnetization transfer effects, are
incorporated in MC simulated Brownian paths and MR signal generation (Section
2.1.4).

  5. 5.

The simulation kernel is accelerated by parallel computation on the GPU
(Section 2.1.5).

####  2.1.1 Substrates and masks

To generate the 3d substrates based on microscopy data for diffusion
simulations, we can translate the voxelized cell segmentation into either
smoothed meshes or binary masks (Nguyen et al., 2018). Each approach has its
own pros and cons.

For the smoothed-mesh approach, the generated cell model has smooth surface,
potentially having surface-to-volume ratio similar to the real cells. However,
it is non-trivial to decide on the degree of smoothing while generating the
cell model. In addition, in simulations, the problem of floating-point
precision may arise, especially for determining whether a random walker
encounters a membrane.

For the binary mask approach, it is fast and simple to translate the discrete
microscopy data into the 2d pixelated or 3d voxelized cell geometries.  In
this study, we only focus on the voxelized geometry of the isotropic voxel
size in 3d. Further, the corresponding simulation kernel is easy to implement
and maintain, and has low computational complexity with minimal problem of
floating-point precision. On the other hand, the generated cell model has
unrealistic surface-to-volume ratio due to its “boxy” cell surface.

To achieve fast and accurate MC simulations of diffusion in microscopy-based
cell geometries, we chose to use the binary mask approach, where the random
walker has at most three interactions with membranes (elastic collision or
membrane permeation) within each step, if the step size is smaller than the
voxel size of the geometry (Fig. 2).  We define the voxel size as the length
of the side of each 3$d$ cube. Applying a small step size smaller than the
voxel size within the voxelized geometry, we only need to take the integer
part of the walker’s position (in a unit of voxel size) to determine where the
walker resides after a particle-membrane interaction, minimizing the precision
problem in numerical calculations.

Here we clarify that the six faces of each 3$d$ voxel do not necessarily
represent the interface (i.e., cell membrane) of compartments; the membranes
always coincide with faces in the voxelized geometry, but not all voxel faces
are therefore part of membranes. When a random walker encounters a face
between two voxels belonging to the same compartment, the random walker will
permeate though the face as it does not exist; in contrast, when a random
walker encounters a face between two voxels, each belonging to a different
compartment, this face is effectively part of the cell membrane, and the check
of permeation probability and particle-membrane interaction will be triggered
(Sections 2.1.2 and 2.1.3).

A proper voxel size of the voxelized geometry should be (1) smaller than the
length scale of cell shapes for accurate simulations, and (2) larger than the
hopping step size in simulations to ensure at most three particle-membrane
interactions in each step. For example, the voxel size of axon geometry should
be much smaller than the axon diameter to capture the fine structures of
cells, such as axon caliber variation and axonal undulation. However, choosing
a small voxel size leads to an even smaller step size and a subsequently short
time for each step (Section 2.1.3), considerably increasing the number of
steps and calculation time. To speed up simulations without losing the
accuracy, it is recommended to choose the voxel size based on simulations in
cell-mimicking simple geometries (e.g., cylinders for axons, spheres for cell
bodies) with known analytical solutions of diffusion signals or metrics.

Finally, for simulations of diffusion in a substrate consisting of multiple
cells, the most computationally expensive calculation is to identify which
compartment a random walker resides in. A way to solve this problem is to
build a lookup table (Yeh et al., 2013; Fieremans and Lee, 2018).  By
partitioning cell geometries into many small voxels, the lookup table records
compartment labels in each voxel. When a random walker hops across few voxels
in a step, we only need to check compartments recorded in these voxels.
Interestingly, each voxel in voxelized geometry records only one compartment
label. In other words, the voxelized geometry serves as the lookup table
itself and thus saves the computational load and memory usage, which could be
non-trivial for the GPU parallelization (Nguyen et al., 2018).

Fig. 2: In the voxelized geometry, a random walker has at most three
interactions with membranes in each step. The simulation of diffusion in 3d
voxelized cell geometry can be simplified as at most three interactions
between a random walker and membranes in each step, when the step size is
smaller than the voxel size; the interactions include either an elastic
collision with or a permeation through the membrane. This simplification is
similar to the concept of the corner reflector (three consecutive coordinate
reflections send vector ${\bf x}\to-{\bf x}$, after which no wall can be
encountered within the corner, representing a given voxel), which has been
applied to radar sets, glass prisms, and bicycle reflectors.

####  2.1.2 Particle-membrane interaction: Why elastic collision is preferred
to other approaches?

In MC simulations of diffusion, the diffusion process is discretized into
multiple steps of constant length and random orientation in 3d, during a
constant time duration for each step. To model the particle-membrane
interaction, the most commonly used implementation for diffusion simulations
is elastic collision (Szafer et al., 1995; Hall and Alexander, 2009), which
properly equilibrates the homogeneous particle density around impermeable,
permeable and absorbing membranes (Fieremans and Lee, 2018). In RMS, we adopt
this approach for simulations of high accuracy.

To reduce the computational load of simulations in complicated geometries,
more and more studies apply other kinds of particle-membrane interactions. The
first alternative is equal-step-length random leap (ERL) (Xing et al., 2013):
a step encountering a membrane is canceled and another direction is chosen to
leap until the new step does not encounter any membranes. In this way, ERL
rejects some steps toward the membrane and effectively repels random walkers
away from the membrane. This repulsive effect leads to inhomogeneous particle
density around membranes (A, Fig. 2); therefore, ERL should not be applied to
simulations of MR contrasts requiring homogeneous particle density, such as
exchange, surface relaxation, and magnetization transfer. For the case of
impermeable, non-absorbing membrane, the bias of diffusivity transverse and
parallel to the membrane due to ERL is proportional to the step size (Figs.
3-4) and could be controlled by choosing a small step size in simulations.

The second alternative for the interaction with membrane is rejection sampling
(Ford and Hackney, 1997; Waudby and Christodoulou, 2011; Nguyen et al., 2018;
Palombo et al., 2018): a step encountering a membrane is canceled, and the
random walker stays still for the step. This simple approach has been shown to
be able to maintain a homogeneous particle density at all times (Szymczak and
Ladd, 2003) and is applicable to simulations of water exchange, as generalized
in B. However, it is rejecting some steps toward the membrane, which leads to
a small bias in diffusivity parallel to membranes (Fig. B.2), and thus it is
still not the best option for an accurate simulation pipeline.

To sum up, the criteria for choosing the particle-membrane interaction in an
MC simulator include (1) the maintenance of homogeneous particle density
around impermeable, permeable, and absorbing membranes, and (2) reliable
simulations of diffusion metrics without interaction-related bias. Rejection
sampling satisfies only the first criterion, and ERL fails both. In contrast,
elastic collision satisfies both criteria, and thus we implement it in RMS for
accurate simulations. Comparisons of other particle-membrane interactions can
be found in (Szymczak and Ladd, 2003; Jóhannesson and Halle, 1996), where
other approaches do not provide benefits in both calculation speed and
accuracy.

####  2.1.3 Impermeable, permeable and absorbing membranes

When a random walker encounters an impermeable membrane within a step, the
random walker is specularly reflected by the membrane, equivalently
experiencing an elastic collision (Szafer et al., 1995). The displacement
before and after the collision are summed up to the step size (Einstein et
al., 1905)

$\delta s=\sqrt{2dD_{0}\cdot\delta t}\,,$ (1)

where $d$ is dimensionality of space, $D_{0}$ is the intrinsic diffusivity,
and $\delta t$ is the time of each MC step.

Further, when a random walker encounters a permeable membrane, the walker has
a probability $1-P$ to be specularly reflected by the membrane, and a
probability $P$ to permeate through. In C, we derive the connection between
$P$ and the membrane permeability $\kappa$ in detail. Briefly, the permeation
probability $P$ can be determined in two ways:

  1. 1.

For the genuine membrane permeability $\kappa$, the permeation probability is
given by

$\displaystyle P$ $\displaystyle\simeq\frac{\kappa_{0}\delta s}{D_{0}}\cdot
C_{d}\,,$ (2a) $\displaystyle{\color[rgb]{0,0,0}C_{d}}$
$\displaystyle\equiv\begin{cases}1&d=1\,,\\\ \pi/4&d=2\,,\\\
2/3&d=3\,,\end{cases}$ (2b)

where $\kappa_{0}$ is the input permeability value (with
$\kappa_{0}\simeq\kappa$ if $P\ll 1$) (Powles et al., 1992; Szafer et al.,
1995; Fieremans and Lee, 2018). When $P$ is not small, the relation of genuine
permeability $\kappa$ and the input value $\kappa_{0}$ is given by Eq. (C.7),
leading to the probability of permeation from compartment 1 ($D_{1}$, $\delta
s_{1}$) to compartment 2 ($D_{2}$, $\delta s_{2}$):

$P_{1\to 2}=\frac{\frac{\kappa\delta s_{1}}{D_{1}}\cdot
C_{d}}{1+\frac{\kappa}{2}\left(\frac{\delta s_{1}}{D_{1}}+\frac{\delta
s_{2}}{D_{2}}\right)\cdot C_{d}}\,.$ (3)

  2. 2.

Alternatively, without assigning a nominal permeability, the permeation
probability can be determined by the ratio of tissue properties on both sides
of the membrane, such as the intrinsic diffusivity and spin concentration,
given by Eq. (C.15) (Baxter and Frank, 2013).

The theoretical background and limitations of the two approaches are discussed
in C.

Finally, for magnetization transfer (MT) effect, its MC simulation can be
effectively modeled as an “absorbing” membrane, or as a surface relaxation
effect. When a random walker encounters an absorbing membrane, the walker has
a probability

$P_{\text{MT}}\simeq\frac{\varrho\,\delta s}{D_{0}}\cdot C_{d}$

to lose all its magnetization (to be saturated) (Fieremans and Lee, 2018),
where $C_{d}$ is the same constant as in Eq. (2a), and $\varrho$ is the
surface relaxivity (with the units of velocity). A random walker’s
magnetization is effectively “saturated” when the exchange with the
macromolecule pool happens with the probability $P_{\text{MT}}$, equivalently
introducing the weighting

$\alpha_{\varrho}=\begin{cases}1&\text{not saturated;}\\\ 0&\text{after
saturation event}\end{cases}$ (4)

for the random walker’s contribution to the net signal.

The signal decay due to MT is caused by the exchange events between water
protons and saturated macromolecule protons during diffusion. The mean-field
estimate for the MT exchange rate $R$ (units of inverse time) between the
liquid pool and the macromolecular pool, with exchange happening at the
surfaces with the net surface-to-volume ratio $S/V$, is given by (Slijkerman
and Hofman, 1998)

$R=\varrho\cdot\frac{S}{V}\,.$

####  2.1.4 Other MR contrast mechanisms

The transverse magnetization experiences the $T_{2}$ NMR relaxation. Hence,
during all the time when the spin magnetization is in the transverse plane,
such as for the conventional spin-echo diffusion sequence, the weighting due
to the $T_{2}$ relaxation for each random walker’s contribution to the overall
signal is (Szafer et al., 1995)

$\alpha_{T}=\exp\left(-\sum_{i}\frac{t_{i}}{T^{(i)}}\right)\,,$ (5)

where $t_{i}$ is the total time of staying in the $i$-th compartment during
the echo time (TE), with the corresponding $T_{2}$ relaxation time
$T^{(i)}\equiv T_{2}^{(i)}$. Similarly, for simulations of a stimulated-echo
sequence, the weighting due to the $T_{1}$ relaxation during the mixing time
(when the magnetization is parallel to $B_{0}$) for each random walker’s
contribution is also given by Eq. (5) (Woessner, 1961), where $t_{i}$ is the
total time of staying in the $i$-th compartment during the mixing time, with
the corresponding $T_{1}$ relaxation time $T^{(i)}\equiv T_{1}^{(i)}$.

Other contrasts, such as susceptibility effect, blood oxygen level dependent
(BOLD) effect and intravoxel incoherent motion (IVIM), can be further added to
RMS by modifying each random walker’s phase factor $e^{-i\varphi(t)}$ due to
the flow velocity and the local Larmor frequency offset in a standard way.

####  2.1.5 Parallel computations and input-output

To achieve precise simulations of diffusion in complicated 3d microgeometry,
it is inevitable to employ a large number of random walkers. However, the
complexity of substrate and the combination of multiple MR contrast mechanisms
in recent studies substantially increase the computational load in
simulations, prompting the usage of parallel computations, either through
multiple nodes on a cluster, multiple threads on a CPU/GPU core, or their
combinations (Waudby and Christodoulou, 2011; Xu et al., 2018; Palombo et al.,
2018; Nguyen et al., 2018; Lee et al., 2020b).

The simulation in RMS is implemented in CUDA C++ and accelerated through the
parallel computation on GPU. Furthermore, by performing multiple simulations
on a GPU cluster, it is possible to further accelerate simulations through
multiple nodes, with each node equipped with a GPU core.

The input RMS files and the shell script to compile and run the CUDA C++
kernel are all generated in a MATLAB script. Moreover, both the input and
output files are text files, and all the codes are open-sourced and platform-
independent. These properties make the RMS congenial to even beginners in this
field to adapt RMS and program their own simulation pipelines.  In Table 1, we
summarize input parameters of RMS.

###  2.2 Monte Carlo simulations in realistic microstructure

Here we shortly summarize the diffusion simulations implemented in RMS:

  1. 1.

Random walkers’ initial positions ${\bf x}_{0}$ are randomly initialized to
achieve a homogeneous particle density. In RMS, we provide users the
flexibility to define a “dead” space, where no random walkers are initialized
and allowed to step in.

  2. 2.

Random walkers diffuse in a continuous space with voxelized microgeometry
(Section 2.1.1). For the $i$-th step, the random walker is at position ${\bf
x}_{i-1}$ before the random hop, and a step vector
${\color[rgb]{0,0,0}\delta}{\bf x}$ of constant length $\delta s$ ($<$ voxel
size) and random direction in 3d is generated.

  3. 3.

The random walker hops to a new position

${\bf x}_{i}={\bf x}_{i-1}+{\color[rgb]{0,0,0}\delta}{\bf x}$ (6)

if it does not encounter the edge of the voxel where it resides in before the
hopping. If the edge of the voxel is encountered, it will either permeate
though or be elastically reflected from the voxel edge based on the permeation
probability defined in Section 2.1.2, leading to a new position ${\bf x}_{i}$
accordingly. The permeation probability is set to 1 if the random walker
encounters the edge between two adjacent voxels belonging to the same
compartment.

  4. 4.

The above particle-membrane interaction will repeat at most three times during
each MC step due to the voxelized geometry and a small step size ($<$ voxel
size), as illustrated in Fig. 2.

  5. 5.

The diffusion metrics and signals are then calculated based on the diffusional
phase due to simulated diffusion trajectory ${\bf x}_{i}$ and signal weighting
($\alpha_{\varrho}$, $\alpha_{T}$ in Sections 2.1.3 and 2.1.4), with details
in the following Section 2.3.

The accuracy of RMS is guaranteed by the theoretical exploration of particle-
membrane interaction (A and B) and water exchange (C).

###  2.3 Diffusion metrics

Normalized diffusion signals of the pulse-gradient spin echo sequence are
calculated based on the diffusional phase $\varphi$ accumulated along the
diffusion trajectory ${\bf x}(\tau)$ (Szafer et al., 1995):

$S=\frac{\langle\alpha\cdot
e^{-i\varphi}\rangle_{N}}{\langle\alpha\rangle_{N}}\,,\quad\varphi=\int_{0}^{\text{TE}}{\bf
g}(\tau)\cdot{\bf x}(\tau)\,\text{d}\tau\,,$ (7)

where $\langle...\rangle_{N}$ is the average over $N$ random walkers,
$\alpha=\alpha_{\varrho}\cdot\alpha_{T}$ is the signal weighting of individual
walker in Eqs. (4) and (5), $\langle\alpha\rangle_{N}$ is the normalization
constant (the total non-diffusion-weighted signal), ensuring that the signal
$S\equiv 1$ in the absence of diffusion weighting $(b=0)$, TE is the echo
time, and ${\bf g}(\tau)$ is the diffusion-sensitizing gradient of the Larmor
frequency.

For the pulsed-gradient sequence, the cumulant expansion of Eq. (7) yields the
diffusivity and kurtosis in the narrow pulse limit, given by

$D=\frac{1}{2t}\cdot\frac{\left\langle\alpha\cdot\left(\Delta{\bf
x}\cdot\hat{\bf
g}\right)^{2}\right\rangle_{N}}{\langle\alpha\rangle_{N}}\,,\quad
K=\frac{\left\langle\alpha\cdot\left(\Delta{\bf x}\cdot\hat{\bf
g}\right)^{4}\right\rangle_{N}}{\left\langle\alpha\cdot\left(\Delta{\bf
x}\cdot\hat{\bf g}\right)^{2}\right\rangle_{N}^{2}}-3\,,$

where $\Delta{\bf x}={\bf x}(t)-{\bf x}(0)$ is the diffusion displacement,
$\hat{\bf g}$ is the gradient direction, and $t$ is the diffusion time. For an
ideal pulse-gradient sequence, $\text{TE}=t$ in the narrow pulse limit.

For the pulsed-gradient sequence with wide pulses (Stejskal and Tanner, 1965),
to obtain the diffusion and kurtosis tensors we simulate signals for multiple
diffusion weightings $b$ according to the gradient wave form ${\bf g}(\tau)$,
and fit the cumulant expansion in the powers of $b$ to the simulated signals
(Jensen et al., 2005):

$\ln S\simeq-b\sum_{ij}n_{ij}D_{ij}\\\
+\tfrac{1}{6}(b\bar{D})^{2}\sum_{ijkl}n_{ijkl}W_{ijkl}+{\cal O}(b^{3})\,,\quad
b=g^{2}\delta^{2}(\Delta-\delta/3)\,,$ (8)

where $n_{ij}=\hat{g}_{i}\hat{g}_{j}$,
$n_{ijkl}=\hat{g}_{i}\hat{g}_{j}\hat{g}_{k}\hat{g}_{l}$, $D_{ij}$ and
$W_{ijkl}$ are diffusion and kurtosis tensors,
$\bar{D}=\frac{1}{3}\sum_{i}D_{ii}$, and the diffusion time $t\approx\Delta$
is roughly the time interval $\Delta$ between two gradient pulses with the
precision of its definition given by the gradient pulse width $\delta$
(Novikov et al., 2019, Sec. 2). The above formula and the overall diffusion
attenuation $S=S[{\bf g}(\tau)]$ can be generalized to be a functional of the
multi-dimensional gradient trajectory ${\bf g}(\tau)$, with $bn_{ij}\to
B_{ij}$ being the $B$-matrix, and with the higher-order terms defiend
accordingly (Topgaard, 2017).

Voxel size | $l_{p}$ |  $<$ length scale in cell geometry (e.g., axon diameter, cell body size)
---|---|---
Intrinsic diffusivity | $D_{0}$ |  $\sim 2-3$ $\mu$m2/ms in vivo
Time for each step | $\delta t$ | small enough such that $P_{1\to 2}\ll 1$, $P_{\text{MT}}\ll 1$, and $\delta s<l_{p}$
Total time | $T_{\text{max}}$ |  $>t$ for narrow PG, $>(t+\delta)$ for wide PG
Number of steps | $N_{\text{step}}$ |  $=\left\lceil T_{\text{max}}/\delta t\right\rceil$, at least 1000 steps to lower discretization error
Number of particles | $N_{\text{par}}$ | as many as possible within the given calculation time
Step size | $\delta s$ |  $=\sqrt{2dD_{0}\cdot\delta t}$, $d=$dimension (typically $d=3$)
Optional |  |
Diffusion gradient | ${\bf g}(\tau)$ | limited by maximal gradient strength, slew rate, and b-value
Permeability | $\kappa$ |
Surface relaxivity | $\varrho$ |
MT exchange rate | $R$ | $=\varrho\cdot S/V$
Relaxation time |  $T_{1}$/$T_{2}$ |  $\sim 800/80$ ms in brain white matter at 3T
Table 1: Input parameters of RMS. The guideline of setting up simulation
parameters is based on (Fieremans and Lee, 2018). MT: magnetization transfer,
PG: pulsed-gradient, $t$: diffusion time, $\delta$: gradient pulse width,
$S/V$: surface-to-volume ratio.

##  3 RMS Applied to Intra-Axonal Microstructure

In this Section, we describe an RMS-compatible example of a realistic electron
microscopy (EM) tissue segmentation (Section 3.1), give an overview of the
related biophysical models (Section 3.2), describe the RMS settings for MC in
axonal geometry (Section 3.3), and outline our results for the diffusion along
(Section 3.4) and transvere (Section 3.5) to the axons.

All procedures performed in studies involving animals were in accordance with
the ethical standards of New York University School of Medicine. All mice were
treated in strict accordance with guidelines outlined in the National
Institutes of Health Guide for the Care and Use of Laboratory Animals, and the
experimental procedures were performed in accordance with the Institutional
Animal Care and Use Committee at the New York University School of Medicine.
This article does not contain any studies with human participants performed by
any of the authors.

###  3.1 Axon segmentation based on electron microscopy

A female 8-week-old C57BL/6 mouse was perfused trans-cardiacally, and the genu
of corpus callosum was fixed and analyzed with a scanning electron microscopy
(Zeiss Gemini 300). We selected a subset of the EM data ($36\times 48\times
20$ $\mu$m3 in volume, $6\times 6\times 100$ nm3 in voxel size), down-sampled
its voxel size to  $24\times 24\times 100$ nm3 slice-by-slice by using Lanczos
interpolation, and segmented the intra-axonal space (IAS) of 227 myelinated
axons using a simplified seeded-region-growing algorithm (Adams and Bischof,
1994).  The segmented IAS mask was further down-sampled into an isotropic
voxel size of ($100$ nm)3 for numerical simulations. More details can be found
in our previous work (Lee et al., 2019). For simulations in RMS, the voxelized
microgeometry based on the axon segmentation is shown in the bottom-right
corner of Fig. 1 (Lee et al., 2020b). In the following examples, we assume
that the tissue properties (i.e., diffusivity, relaxation time) are the same
in cytoplasm and organelles (e.g., mitochondria).

###  3.2 Tissue parameters and biophysical models for the axonal geometry

To quantify the axon geometry, we define axon’s inner diameter, caliber
variation, and axonal undulation as follows: The inner diameter $2r\equiv
2\sqrt{\Omega/\pi}$ of an axon cross-section is defined as that of an
equivalent circle with the same cross-sectional area $\Omega$ (West et al.,
2016); the caliber variation is defined as the coefficient of variation of
radius (ratio of standard deviation to mean), $\text{CV}(r)$ (Lee et al.,
2019); the axonal undulation $|{\bf w}|$ is defined as the shortest distance
between the axonal skeleton and the axon’s main axis (Lee et al., 2020a);
based on the analysis of an axonal skeleton, we can roughly estimate the
length scale of undulation amplitude $w_{0}$ and wavelength $\lambda$ using a
simplified 1-harmonic model (D).

In our previous studies, we showed that restrictions along axons are randomly
positioned with a finite correlation length $\ell_{c}$ of a few micrometers
(Lee et al., 2020b). This short-range disorder leads to the diffusivity along
axons

$D_{\parallel}(t)\simeq D_{\infty}+c\cdot t^{-1/2}\,,\quad t\gg
t_{c}={\ell_{c}^{2}\over 2D_{\infty}}$ (9)

approaching its $t\to\infty$ limit $D_{\infty}$ in a power-law fashion with
exponent $1/2$ (Novikov et al., 2014; Fieremans et al., 2016; Lee et al.,
2020b). Here, $c$ is the strength of the restrictions; the relation (9) holds
in the narrow pulse limit, and acquires further corrections for finite pulse
width $\delta$. The scaling (9) becomes valid when the diffusion length
exceeds the correlation length $\ell_{c}$ of the structural disorder (e.g.,
beads) along axons.

The bulk diffusivity along axons correlates with axon geometry via (Lee et
al., 2020b)

$\zeta\equiv\frac{D_{0}-D_{\infty}}{D_{\infty}}\propto\text{CV}^{2}(r)\,,$
(10)

indicating that the stronger the caliber variation $\text{CV}(r)$, the smaller
$D_{\infty}$.

Furthermore, the nonzero diffusivity transverse to axons is contributed by
axon caliber and undulation (Nilsson et al., 2012; Brabec et al., 2020; Lee et
al., 2020a). In the wide pulse limit, untangling the two effects in diffusion
transverse to axons is non-trivial. Instead, we translate the radial
diffusivity $D_{\perp}$ into the effective radius measured by MR,
$r_{\text{MR}}$ (Burcaw et al., 2015; Veraart et al., 2020), and identify
which contribution dominates (Lee et al., 2020a):

$\displaystyle r_{\text{MR}}$
$\displaystyle\equiv\left[\frac{48}{7}\delta\left(t-\frac{\delta}{3}\right)D_{0}D_{\perp}\right]^{1/4}$
(11) $\displaystyle\simeq\begin{cases}r_{\text{cal}}\equiv\left(\langle
r^{6}\rangle/\langle r^{2}\rangle\right)^{1/4}\,,&t_{D}\ll\delta\ll t_{u}\\\
r_{\text{und}}\equiv\left(\frac{6}{7\pi^{2}}\right)^{1/4}\cdot\sqrt{w_{0}\lambda}\,,&\delta\gg
t_{u}\end{cases}$ (12)

where $r_{\text{cal}}$ and $r_{\text{und}}$ are the apparent axon size due to
the axon caliber and undulation respectively, and

$t_{D}\equiv r^{2}/D_{0}\quad\mbox{and}\quad
t_{u}\equiv\lambda^{2}/(4\pi^{2}D_{0})$ (13)

are correlation times for the axon caliber and undulation respectively. In
realistic axonal shapes, $r\sim 0.5$ $\mu$m and $\lambda\sim$ 20 $\mu$m lead
to $t_{D}\sim 0.1$ ms and $t_{u}\sim 5$ ms for $D_{0}=2$ $\mu$m2/ms (Lee et
al., 2019, 2020a).

To factor out the effect of undulation, diffusion signals are spherically
(orientationally) averaged for each $b$-shell. Assuming that axonal segments
have cylindrical shapes, the spherically averaged intra-axonal signals are
given by (Jensen et al., 2016; Veraart et al., 2019, 2020)

$\bar{S}(b)\simeq\beta\,e^{-bD_{\perp}+{\cal O}(b^{2})}\,b^{-1/2}\,,$ (14)

where $\beta\equiv\sqrt{\pi/(4D_{a})}$, and $D_{a}$ is the diffusivity along
axon segments within the diffusion length $\sim\sqrt{D_{a}t}$. The effective
radius is then calculated based on $D_{\perp}$ and Eq. (11).

###  3.3 RMS settings for MC simulations in the intra-axonal space

To demonstrate the relation of diffusion metrics and microgeometry in intra-
axonal space in Eqs. (10) and (12), we performed diffusion simulations in 227
realistic axons segmented from mouse brain EM in Section 3.1: random walkers
diffused over $1\times 10^{6}$ steps with a duration $\delta t=2\times
10^{-4}$ ms and a length $\delta s=0.049$ $\mu$m (Eq. (1) with $D_{0}=2$
$\mu$m2/ms) for each of the following simulations of monopolar pulsed-gradient
sequences.

For simulations of narrow pulse sequence, $1\times 10^{7}$ random walkers per
axon were applied. The apparent diffusivity $D_{\parallel}(t)$ along axons was
calculated by using diffusion displacement along the axons’ main direction
(z-axis) at diffusion time $t\leq 100$ ms. Eq. (9) for the time range
$t=60-80$ ms was fit to the simulated $D_{\parallel}(t)$, and the fit
parameter $D_{\infty}$ was correlated with the caliber variation
$\text{CV}(r)$ in Eq. (10).

For simulations of wide pulse sequence, $2\times 10^{4}$ random walkers per
axon were applied. The duration between pulses $\Delta=1-100$ ms was equal to
the gradient pulse width $\delta$, such that diffusion time
$t\approx\Delta=\delta$. The diffusion signals were calculated based on the
accumulated diffusion phase (7) for 10 b-values $=0.2-2$ ms/$\mu$m2 along
directions (x- and y-axes) transverse to axons’ main direction (z-axis). The
simulated apparent diffusivity $D_{\perp}$ transverse to each axon was
translated to the effective radius by using Eq. (11) and compared with the
contribution of caliber variations and axonal undulations in Eq. (12).

For simulations of directionally averaged signal, $2\times 10^{4}$ random
walkers per axon were applied. The diffusion time and gradient pulse width
$(t,\delta)=(20,7.1)$, $(30,13)$ and $(50.9,35.1)$ ms were used to match the
experimental settings on animal 16.4T MR scanner (Bruker BioSpin), Siemens
Connectome 3T MR scanner, and clinical 3T MR scanner (Siemens Prisma) (Veraart
et al., 2020; Lee et al., 2020a). Diffusion signals were calculated based on
the accumulated diffusional phase for 18 b-values $=16-100$ ms/$\mu$m2 along
30 uniformly distributed directions for each b-shell. Eq. (14) was fit to the
spherically averaged signal $\bar{S}(b)$ from all axons (volume-weighted sum),
and the fit parameter $D_{\perp}$ was again translated to the effective radius
in Eq. (11) and compared with the contribution of caliber variations and
axonal undulations of all axons, i.e., $r_{\text{cal}}$ in Eq. (12) and
$\tilde{r}_{\text{und}}$ defined below.

The $r_{\text{und}}$ in Eq. (12) is defined for the individual axon and does
not take the volume differences between axons into account. To add up the
contribution of $D_{\perp}$ due to undulations of all axons, we define a
volume-weighted average of $r_{\text{und}}$:

$\tilde{r}_{\text{und}}=\left(\sum_{i}f_{i}\cdot
r_{\text{und},i}^{4}\right)^{1/4}\,,$ (15)

with the $i$-th axon’s volume fraction $f_{i}$, such that $\sum_{i}f_{i}=1$.

Fig. 3: Simulated apparent diffusivity parallel and transverse to 227
realistic axons segmented from a mouse brain corpus callosum. a-b In the
narrow pulse limit, the simulated time($t$)-dependent diffusivity
$D_{\parallel}(t)$ along axons scales as $1/\sqrt{t}$ in Eq. (9), and its bulk
diffusivity $D_{\infty}$ correlates with the caliber variation $\text{CV}(r)$
via Eq. (10). c In the wide pulse limit ($t=\delta$), the simulated
diffusivity $D_{\perp}(\delta)$ transverse to axons scales as $\delta^{-2}$
(dashed line) at very long times $\delta\gtrsim 50$ ms. d Based on the
analysis of realistic axonal skeleton (D), the estimated undulation wavelength
$\lambda$ positively correlates with the undulation amplitude $w_{0}$. e-f
Translating the simulated $D_{\perp}$ into an effective radius $r_{\text{MR}}$
measured by MR via Eq. (11), the value of $r_{\text{MR}}$ is compared with the
contributions of caliber variations and undulations, $r_{\text{cal}}$ and
$r_{\text{und}}$ in Eq. (12) respectively. The Pearson’s correlation
coefficient $R$ of $r_{\text{MR}}$ and each contribution is also shown as a
function of time $\delta$ (insets). The result in panel b is similar to Fig.
2d in (Lee et al., 2020b), whereas the short $T_{2}$ and slow diffusivity in
mitochondria are ignored here to focus on the effect of axonal shape. The
results in panels c and e-f are similar, but not identical, to the ones
obtained using ERL algorithm in Fig. 7 of (Lee et al., 2020a). Here, these
simulations were re-done using elastic collision with impermeable membranes,
and the estimation of undulation amplitude and wavelength is purely based on
the geometry of axonal skeleton, as opposed to using ERL in simulations and
fitting a 1-harmonic undulation model to simulated $D_{\perp}$ with fit
parameters $w_{0}$ and $\lambda$ in (Lee et al., 2020a). The effective theory
parameters and scaling relations are essentially the same for both MC
implementations (elastic collision vs ERL) because only diffusion simulations
within impermeable cells were performed.

###  3.4 Results: Diffusion metrics and caliber variations along axons

In 227 realistic axons segmented from a mouse brain, the simulation result of
narrow pulse sequence shows that the apparent diffusivity $D_{\parallel}(t)$
along axons scales as $1/\sqrt{t}$ in Eq. (9) (Fig. 3a), and the fitted bulk
diffusivity $D_{\infty}$ correlates with the caliber variation $\text{CV}(r)$
in Eq. (10) (Fig. 3b). This demonstrates that the diffusion along axons is
characterized by the short-range disorder in 1d, corresponding to randomly
positioned beads along axons (Novikov et al., 2014; Fieremans et al., 2016;
Lee et al., 2020b).

###  3.5 Results: Diffusion metrics transverse to axons

We compare the contribution of the axon caliber and axonal undulation to the
diffusivity $D_{\perp}$ transverse to axons (Lee et al., 2020a): At clinical
diffusion times ($t\sim\delta>10$ ms) of conventional wide-pulse sequences,
the apparent diffusivity transverse to axons is dominated by the contribution
of undulations (Fig. 3c-f). To release the requirement of very short diffusion
time $\sim 1$ ms for accurate axon diameter mapping, the modeling of
spherically averaged signals at multiple diffusion weightings are employed to
partially factor out the undulation effect (Fig. 4).

####  3.5.1 Diffusivity and axonal undulations

In realistic axons of a mouse brain, the simulation result of wide pulse
sequences shows that the apparent diffusivity $D_{\perp}$ transverse to axons
scales as $\delta^{-2}$ at very long time $\delta\gtrsim 50$ ms in most axons
(Fig. 3c). This scaling is consistent with both the undulation and caliber
contributions, since

$D_{\perp}\simeq-\tfrac{1}{b}\ln S\propto{r_{\rm
cal}^{4}\delta/D_{0}\over\delta^{2}(\Delta-\delta/3)}\sim\delta^{-2}\,,\quad\Delta=\delta\,.$

However, the onset of this scaling depends on the correlation time, Eq. (13).
If the caliber dominates, the $\delta^{-2}$ scaling will occur for $\delta\gg
t_{D}\sim 0.1$ ms (Neuman, 1974). However, the actual scaling happens at times
about 2 orders of magnitude longer. The onset for the undulations, $\delta\gg
t_{u}\sim 5$ ms, as estimated after Eq. (13), is far more consistent with our
RMS results.

To better understand the length scale of undulations, the undulation amplitude
$w_{0}$ and wavelength $\lambda$ are estimated based on the axonal skeleton
(D) (Fig. 3d). The scatter plot of $w_{0}$ and $\lambda$ shows positive
correlation; in other words, the longer undulation wavelength, the stronger
undulation amplitude. Furthermore, the slope
$\text{d}w_{0}/\text{d}\lambda\simeq 0.039$ and the simplified 1-harmonic
undulation model in Eq. (D.2) indicate an estimate of the intra-axon
undulation dispersion $\simeq 10^{\circ}$, which is about a factor of 2
smaller than the inter-axon fiber orientation dispersion $\simeq 22^{\circ}$,
leading to an overall dispersion angle
$\simeq\sqrt{(10^{\circ})^{2}+(22^{\circ})^{2}}\simeq 24^{\circ}$ (D) (Ronen
et al., 2014; Lee et al., 2019). Interestingly, the inter-axon dispersion
$\sim 22^{\circ}$ and overall dispersion $\sim 24^{\circ}$ are consistent with
the numerical calculation of time-dependent dispersion at long and short times
respectively in (Lee et al., 2019), where the same group of axons were
analyzed.

Furthermore, the simulated $D_{\perp}$ is translated to the effective radius
$r_{\text{MR}}$ measured by MR via Eq. (11) and compared with the contribution
of axon caliber and axonal undulation, i.e., $r_{\text{cal}}$ and
$r_{\text{und}}$ in Eq. (12) respectively (Fig. 3e-f). At very short time
$\delta\sim 1$ ms, the MR estimate $r_{\text{MR}}$ coincides with the
contribution of axon caliber $r_{\text{cal}}$; however, at long time
$\delta\sim 100$ ms, the MR estimate $r_{\text{MR}}$ is consistent with the
contribution of axonal undulations $r_{\text{und}}$. Finally, at clinical time
range $\delta\sim 10-50$ ms, the contribution of undulations $r_{\text{und}}$
shows much higher correlation with the MR estimate $r_{\text{MR}}$, compared
with the correlation of $r_{\text{MR}}$ and $r_{\text{cal}}$ (the insets of
Fig. 3e-f). These observations demonstrate that axonal undulations confound
the axonal diameter mapping at clinical time range and need to be factored out
for an accurate diameter estimation.

Fig. 4: Simulated spherically averaged signal of all 227 realistic axons
segmented from a mouse brain corpus callosum. The simulated spherically
averaged signal $\bar{S}(b)$ of all axons scales as $1/\sqrt{b}$ in Eq. (14),
and its negative intercept at $b\to\infty$ indicates an estimate of
diffusivity $D_{\perp}$ transverse to axons. For $(t,\delta)=(20,7.1)$ ms
(animal scanner), the negative intercept of $\bar{S}(b)$ corresponds to a
$D_{\perp}\simeq 1.6\times 10^{-4}$ $\mu$m2/ms and an MR-measured radius
$r_{\text{MR}}\simeq 0.72$ $\mu$m based on Eq. (11); however, for
$(t,\delta)=(30,13)$ and $(50.9,35.1)$ ms (human scanner), the positive
intercepts of $\bar{S}(b)$ leads to the non-physical, negative $D_{\perp}$.
The MR-measured radius $r_{\text{MR}}$ at $(20,7.1)$ ms is close to the
caliber contribution from histology, $r_{\text{cal}}=0.79$ $\mu$m in Eq. (12),
and much smaller than the undulation contribution from histology,
$\tilde{r}_{\text{und}}=2.15$ $\mu$m in Eq. (15), indicating that the
undulation contribution is factored out by orientational average. The result
of left panel is similar, but not identical, to the ones obtained using ERL
algorithm in Fig. 8g of (Lee et al., 2020a). Here, these simulations were re-
done using elastic collision with impermeable membranes. The effective theory
parameters and scaling relations are essentially the same for both MC
implementations because only simulations within impermeable axons were
performed.

####  3.5.2 Spherically averaged signals and axonal diameter mapping

The spherically averaged signal $\bar{S}(b)$ of all axons from a mouse brain
scales as $1/\sqrt{b}$ in Eq. (14) (Fig. 4), and its negative intercept at
$b\to\infty$ corresponds to an estimate of radial diffusivity $D_{\perp}$,
whose contribution of undulations is partially factored out: On the one hand,
at short time $(t,\delta)=(20,7.1)$ ms (animal scanner), the $\bar{S}(b)$ has
a negative intercept, indicating a $D_{\perp}\simeq 1.6\times 10^{-4}$
$\mu$m2/ms and an MR-measured radius $r_{\text{MR}}\simeq 0.72$ $\mu$m in Eq.
(11), which is close to the histology-based caliber contribution
$r_{\text{cal}}=0.79$ $\mu$m in Eq. (12) and smaller than the histology-based
undulation contribution $r_{\text{und}}=2.15$ $\mu$m in Eq. (15). The fitted
diffusivity along axon segments, $D_{a}\simeq 0.19\cdot D_{0}$, is
unexpectedly low since the variation of local diffusivities along axons is not
considered in Eq. (14). On the other hand, at longer diffusion times
$(t,\delta)=(30,13)$ ms and $(50.9,35.1)$ ms (human scanners), the
$\bar{S}(b)$ has a positive intercept, leading to the non-physical, negative
$D_{\perp}$.

The simulation result demonstrates that the requirement of short diffusion
time $\sim 1$ ms for an accurate axon size estimation can be partly released
by spherically averaging signals.

##  4 Outlook

Performing MC simulations in realistic cell geometries using the proposed RMS
helps to test the sensitivity of diffusion MRI to tissue features and validate
the biophysical models. RMS is an open-source platform for the Monte Carlo
simulations of diffusion in realistic tissue microstructure.  In addition to
the examples of diffusion within intra-cellular space, it is also possible to
perform simulations of diffusion in the extra-cellular space, as well as of
the exchange between intra- and extra-cellular spaces. The tissue preparation
preserving extra-cellular space in histology is non-trivial, prompting the
development of pipelines to generate semi-realistic tissue microstructure by
packing multiple artificial “cells”, such as MEDUSA (Ginsburger et al., 2019)
and ConFiG (Callaghan et al., 2020), whose generated microgeometry potentially
could be transformed to 3d voxelized data, compatible with RMS.

Beyond the proposed pipeline, here we provide an outlook for the next
generation simulation tools. So far, most of the implementations of diffusion
simulations focus on the MR sequence with radiofrequency pulses of 90∘ and
180∘, such as the spin-echo (Stejskal and Tanner, 1965) and stimulated-echo
(Tanner, 1970) sequences. To simulate other MR sequences, such as the steady-
state free precession sequence (McNab and Miller, 2008), it is required to
combine the diffusion simulation with the Bloch simulator in MR system.
Furthermore, the generation of accurate cell segmentation for simulations in
microscopy-based geometry is time-consuming and labor-intensive; performing
diffusion simulations directly in a mesoscopic diffusivity map, obtained via a
transformation of the microscopy intensity, could largely simplify the
validation pipeline — however, such segmentation-less simulation must itself
be thoroughly validated.

As tissue microstructure imaging with MRI is a vastly expanding research area
of quantitative MRI, the synergies between microstructure model validation,
hardware and acquisition techniques have a potential to transform MRI into a
truly quantitative non-invasive microscopy/histology technique, as discussed
in (Novikov, 2020). Hence, the development of next-generation simulation tools
will benefit not only the field of microstructural MR imaging, but the whole
MRI community.

##  5 Conclusions

Numerical simulations in realistic 3d microgeometry based on microscopy data
serve as a critical validation step for biophysical models, in order to obtain
quantitative biomarkers, e.g., axonal diameter, the degree of caliber
variations and axonal undulations, for potential clinical applications. With
the help of the proposed RMS pipeline, it is possible to achieve fast and
accurate simulations of diffusion in realistic tissue microstructure, as well
as the interplay with other MR contrasts. RMS enables ab initio simulations
from realistic microscopy data, facilitating model validation and experiment
optimization for microstructure MRI in both clinical and preclinical settings.

## Acknowledgements

We would like to thank Sune Jespersen for fruitful discussions about theory,
and the BigPurple High Performance Computing Center of New York University
Langone Health for numerical computations on the cluster. Research was
supported by the National Institute of Neurological Disorders and Stroke of
the NIH under awards R01 NS088040 and R21 NS081230, by the National Institute
of Biomedical Imaging and Bioengineering (NIBIB) of the NIH under award number
U01 EB026996, and by the Irma T. Hirschl fund, and was performed at the Center
of Advanced Imaging Innovation and Research (CAI2R, www.cai2r.net), a
Biomedical Technology Resource Center supported by NIBIB with the award P41
EB017183.

Fig. 1 is adapted from (Chin et al., 2002; Xu et al., 2018; Nguyen et al.,
2018; Palombo et al., 2019; Lee et al., 2020b) with permission from Wiley,
Elsevier, and Springer Nature.

Fig. 1 is adapted from (Fieremans and Lee, 2018) with permission from
Elsevier.

## Declarations of interest

None.

## Data and code availability statement

The SEM data and axon segmentation can be downloaded on our web page
(www.cai2r.net/resources/software).

The source codes of Monte Carlo simulations can be downloaded on our Github
page  (github.com/NYU-DiffusionMRI/monte-carlo-simulation-3D-RMS). The first
release of RMS supports simulations of multiple MR contrasts (diffusion,
$T_{2}$ relaxation, water exchange) and the pulsed-gradient spin-echo
sequence. Simulations of other MR contrasts (MT and $T_{1}$ relaxation),
sequences (stimulated-echo), arbitrary diffusion gradient waveforms, as well
as flow and susceptibility-induced effects will be updated in the future.

##  Appendix A Unbiased simulation of particle-membrane interaction — equal-
step-length random leap

To largely reduce the computational load of MC simulations, Xing et al. (2013)
proposed to model the interaction of random walkers and membranes by using the
equal-step-length random leap (ERL): a step crossing the membrane is canceled,
and another direction is then chosen to leap until the new step does not cross
any membranes.

However, random walkers are effectively repulsed away from the membrane since
the step encountering the membrane is canceled by the algorithm. This
repulsion effect close to the membrane in a thickness of the step size $\delta
s$ introduces the bias of diffusivity transverse and parallel to the membrane,
as discussed below.

###  A.1 Bias in the pore size estimation caused by equal-step-length random
leap

To understand the origin of the bias caused by ERL, we firstly discuss the
particle density between two parallel planes, composed of two impermeable
membranes with a spacing $a$. Due to the time-reversal symmetry of diffusion,
the particle density $\rho(h)$ with a distance $h$ to the membrane is
proportional to the fraction $p_{\Omega}(h)$ of the hopping orientation, along
which particles do not encounter the membrane (blue solid angle in Fig. 1):

$\displaystyle p_{\Omega}(h\leq\delta s)$
$\displaystyle=\begin{cases}\frac{1}{2}&d=1\,,\\\
1-\frac{1}{\pi}\cos^{-1}\left(\frac{h}{\delta s}\right)&d=2\,,\\\
\frac{1}{2}\left(1+\frac{h}{\delta s}\right)&d=3\,,\end{cases}$ (A.1)
$\displaystyle p_{\Omega}(h>\delta s)$ $\displaystyle=1\,.$

The fact that $\rho(h)\propto p_{\Omega}(h)$ in ERL is demonstrated by
simulations of diffusion between two parallel planes in 1d, 2d, and 3d in Fig.
2, where the density $\tilde{\rho}(h)$ is normalized by the number of random
walkers in simulations, and the step size $\delta s$ in Eq. (1) is tuned by
varying the time step $\delta t=0.1-2.5$ $\mu$s with the intrinsic diffusivity
$D_{0}=2$ $\mu$m2/ms.

The knowledge of the particle density further enables us to estimate the bias
in the pore size estimation due to ERL. Considering the diffusion between two
parallel planes at $x=\pm a/2$, the diffusion signal measured by using narrow-
pulse monopolar sequence is given by (Callaghan, 1993)

$S(q,t)=\int\rho_{0}(x)\,{\cal{G}}(x,x^{\prime};t)\,e^{-iq\cdot(x-x^{\prime})}\text{d}x\,\text{d}x^{\prime}\,,$
(A.2)

where $q$ is the diffusion wave vector normal to the plane surface,
$\rho_{0}(x)$ is the particle density at $t=0$, and $\cal G$ is the diffusion
propagator, i.e., the probability of a spin at $x$ diffusing to $x^{\prime}$
during the time $t$.

In MC simulations, on the one hand, we initialize a homogeneous particle
density at $t=0$, i.e., $\rho_{0}(x)=$ const if $x\in[-a/2,a/2]$. On the other
hand, at long times ($t\gg t_{D}=a^{2}/D_{0}$), the random walker loses its
memory of the initial position and has an equal probability of being anywhere
between two planes (Callaghan, 1993):

${\cal G}(x,x^{\prime};t\gg t_{D})\simeq\rho(x)\propto
p_{\Omega}\left(h(x)\right)\,,$

where $p_{\Omega}(h)$ is given by Eq. (A.1) with the distance $h(x)=a/2-|x|$.
Substituting into Eq. (A.2), we obtain

$S(q,t\gg t_{D})\propto\rho_{0}(q)\cdot p_{\Omega}(q)\,,$ (A.3)

where $\rho_{0}(q)$ and $p_{\Omega}(q)$ are the corresponding Fourier
transform quantities:

$\rho_{0}(q)=\text{sinc}\left(\frac{qa}{2}\right)\,,$ (A.4)

and

$p_{\Omega}(q)=\frac{1}{2}\left[\text{sinc}\left(\frac{qa}{2}\right)+\text{sinc}\left(q\cdot\left(\frac{a}{2}-\delta
s\right)\right)\right]\,,$ (A.5)

for 1-dimensional parallel planes, and

$p_{\Omega}(q)=\frac{1}{2}\left[\text{sinc}\left(\frac{qa}{2}\right)\cdot(J_{0}(q\delta
s)+1)\right.\\\ \left.-\frac{\cos(qa/2)}{qa/2}\cdot H_{0}(q\delta
s)\right]\,,$ (A.6)

for 2-dimensional parallel planes with $H_{0}(\cdot)$ the Struve function and
$J_{0}(\cdot)$ the Bessel function of the first kind, and

$p_{\Omega}(q)=\frac{1}{2}\left[\text{sinc}\left(\frac{qa}{2}\right)+\text{sinc}\left(\frac{q\delta
s}{2}\right)\cdot\frac{\sin\left(q\cdot(a-\delta s)/2\right)}{qa/2}\right]\,,$
(A.7)

for 3-dimensional parallel planes. Substituting Eqs. (A.4), (A.5), (A.6) and
(A.7) into Eq. (A.3) and applying Taylor series for $q\delta s\ll qa\ll 1$, we
have

$-\ln\left(\frac{S}{S_{0}}\right)_{t\gg
t_{D}}\simeq\begin{cases}\frac{1}{12}q^{2}\cdot\left(a^{2}-a\delta
s\right)&d=1\,,\\\ \frac{1}{12}q^{2}\cdot\left(a^{2}-\frac{2}{\pi}a\delta
s\right)&d=2\,,\\\ \frac{1}{12}q^{2}\cdot\left(a^{2}-\frac{1}{2}a\delta
s\right)&d=3\,,\end{cases}$

where $S_{0}\equiv S|_{q\to 0}$.

Using the definition of the apparent diffusivity,
$D\equiv-\frac{1}{b}\ln\frac{S}{S_{0}}$ with $b=q^{2}t$ in narrow pulse, the
diffusivity transverse to parallel planes in simulations of ERL is given by
$D(t)=(a^{\prime})^{2}/(12t)$, where

$\frac{a-a^{\prime}}{2}\equiv\delta r\simeq\begin{cases}\frac{1}{4}\delta
s&d=1\,,\\\ \frac{1}{2\pi}\delta s&d=2\,,\\\ \frac{1}{8}\delta
s&d=3\,.\end{cases}$ (A.8)

Comparing with the unbiased solution $D(t)=a^{2}/(12t)$ (Callaghan, 1993), the
bias of the pore size estimation due to ERL can be considered as an effective
pore shrinkage, and the $\delta r$ serves as the shrinkage distance from the
membrane. To demonstrate Eq. (A.8), we performed MC simulations of diffusion
between two parallel planes in 1d, 2d, and 3d (implemented by using ERL), as
shown in Fig. 3, where the step size $\delta s$ is tuned by varying the time
step $\delta t$, as in Fig. 2.

Fig. 1: The solid angle of the hopping orientation, along which particles do
not encounter the membrane (blue). The step size $\delta s$ and the distance
$h$ to the membrane ($h\leq\delta s$) are shown in a 1d, b 2d, and c 3d. The
figure is adapted from (Fieremans and Lee, 2018) with the permission from
Elsevier.  Fig. 2: Inhomogeneous particle density close to the membrane due to
equal-step-length random leap. For elastic collisions, random walkers can not
only leap away from the membrane but also leap toward and be reflected by the
membrane. However, for the equal-step-length random leap, only the steps
jumping away from the boundary are allowed; this repulsive effect results in
lower particle density close to the membrane in a thickness of the step size
$\delta s$ (indicated by vertical lines). The theoretical prediction of this
effect (left column) is consistent with the simulation result (right column)
between two parallel planes in a distance $a$ in 1d, 2d and 3d. Fig. 3: The
bias of diffusivity transverse to membranes due to equal-step-length random
leap. The inhomogeneous particle density due to equal-step-length random leap
causes the bias in diffusivity $D(t)$ transverse to membranes and subsequently
the pore shape estimation, for example, the distance between two parallel
planes. The ground truth of the planes distance is $a=1$ $\mu$m, and its
biased estimation based on $D(t)$ is $a^{\prime}$. This bias $a-a^{\prime}$ is
roughly proportional to the step size $\delta s$.

###  A.2 Bias in the diffusivity parallel to the membrane caused by equal-
step-length random leap

Considering the diffusion in 3d parallel to a membrane and its simulation
implemented with ERL, the random walker close to the membrane has a higher
probability to leap parallel to than perpendicular to the membrane due to the
repulsive effect of ERL in A.1. As a result, the diffusion displacement as
well as the diffusivity $D_{\parallel}(t)$ parallel to the membrane are
slightly overestimated.

Given that a random walker close to a membrane with a distance $h\leq\delta s$
leaps into a direction of a polar angle $\theta$ and an azimuthal angle $\phi$
(spherical coordinate with zenith direction along $x$-axis), its diffusion
displacements parallel to the membrane (along, e.g., $z$-axis) is $\delta
z=\delta s\sin\theta\cos\phi$. Averaging over all possible directions allowed
by the ERL, the second order cumulant of $\delta z$ is given by

$\displaystyle\langle\delta z^{2}\rangle_{\Omega}$
$\displaystyle=\int_{\cos^{-1}(\frac{h}{\delta s})}^{0}\int_{0}^{2\pi}(\delta
s\sin\theta\cos\phi)^{2}\,\frac{\sin\theta\,\text{d}\theta}{1+\frac{h}{\delta
s}}\,\frac{\text{d}\phi}{2\pi}$ $\displaystyle=\tfrac{1}{6}(2\delta
s^{2}+h\delta s-h^{2})\,,\quad h\leq\delta s\,.$

Further averaging over the thickness $\delta s$ surrounding the membrane with
the consideration of the inhomogeneous density $\rho(h)\propto p_{\Omega}(h)$
in 3d in Eq. (A.1), we obtain the second order cumulant of $\delta z$, given
by

$\langle\delta z^{2}\rangle_{h\leq\delta s}=\frac{\int_{0}^{\delta
s}\langle\delta z^{2}\rangle_{\Omega}\cdot
p_{\Omega}\,\text{d}h}{\int_{0}^{\delta
s}p_{\Omega}\,\text{d}h}=\frac{13}{36}\delta s^{2}\,,$

leading to a distinct diffusivity parallel to the membrane for walkers close
to the membrane:

$D_{0}^{\prime}=\frac{\langle\delta z^{2}\rangle_{h\leq\delta s}}{2\delta
t}=\frac{13}{12}D_{0}\,,$

where the step size in Eq. (1) in 3d is applied.

Random walkers surrounding the membrane ($h\leq\delta s$) have a diffusivity
$D_{0}^{\prime}$ parallel to the membrane and a (density-weighted) volume
fraction

$\displaystyle f^{\prime}$ $\displaystyle=\int_{0}^{\delta
s}\rho(h)\cdot\frac{S}{V}\,\text{d}h$ $\displaystyle\simeq\int_{0}^{\delta
s}p_{\Omega}(h)\cdot\frac{S}{V}\,\text{d}h=\frac{3}{4}\cdot\frac{S}{V}\delta
s\,,$

where $\frac{S}{V}$ is the surface-to-volume ratio. Similarly, random walkers
away from the membrane have a diffusivity $D_{0}$ parallel to the membrane and
a volume fraction $f=1-f^{\prime}$. Their volume-weighted sum yields the
overall diffusivity parallel to the membrane, given by

$\displaystyle D_{\parallel}$ $\displaystyle=f\cdot D_{0}+f^{\prime}\cdot
D_{0}^{\prime}$ $\displaystyle\simeq
D_{0}\cdot\left(1+\frac{1}{16}\frac{S}{V}\delta s\right)>D_{0}\,,$ (A.9)

where the correction term is negligible when a small step size is applied.

The above Eq. (A.9) was demonstrated by performing MC simulations of diffusion
between two parallel planes in 3d (implemented by using ERL), as shown in Fig.
4, where the step size $\delta s$ is tuned by varying the time step $\delta
t$, as in Fig. 2.

Fig. 4: The bias of diffusivity parallel to membranes due to equal-step-length
random leap. The inhomogeneous particle density due to equal-step-length
random leap causes the bias in diffusivity $D_{\parallel}$ parallel to the
membrane (e.g., two parallel planes in 3d). This bias is roughly proportional
to the step size $\delta s$ and closely related to the surface-to-volume ratio
$S/V$ of the geometry.

##  Appendix B Unbiased simulation of particle-membrane interaction —
rejection sampling

Another alternative of elastic collision is rejection sampling (Ford and
Hackney, 1997; Waudby and Christodoulou, 2011; Nguyen et al., 2018; Palombo et
al., 2018): a step crossing a membrane is canceled, and the random walker does
not move for the step. This simple approach properly maintains a homogeneous
particle density in simulations of diffusion in the substrate composed of
impermeable membranes (Szymczak and Ladd, 2003) and can be generalized to the
case of permeable membranes as follows.

###  B.1 Generalization of rejection sampling to the case of permeable
membranes

Here, we use the theoretical framework similar to that in Appendix 1 and
Figure 7 of (Fieremans et al., 2010): Given that, in 1d, a permeable membrane
is positioned at $x_{M}=0$, and the intrinsic diffusivity at both sides is
$D_{0}$, the probability of a random walker showing at position $x_{P}$ on the
right side of the membrane at time $t+\delta t$, such that $x_{P}-x_{M}<\delta
s$, is related to the probabilities of previous steps at time $t$, i.e.,

${\cal P}(x_{P},t+\delta t)=\frac{1}{2}{\cal P}(x_{P}+\delta
s,t)+\frac{P_{T}}{2}{\cal P}(x_{P}-\delta s,t)+\frac{1-P_{T}}{2}{\cal
P}(x_{P},t)\,,$ (B.1)

where $\delta t$ is the time step, $\delta s$ is the step size in Eq. (1) in
1d, $P_{T}$ is the permeation probability, and ${\cal P}(\cdot)$ is the
probability density function (PDF) of the particle population at a given
position and time. Notably, the third right-hand-side term in Eq. (B.1)
corresponds to the canceled step without update due to the rejection sampling.
This minor difference from the Eq. [40] in (Fieremans et al., 2010) leads to
slightly different result in the permeation probability.

The PDF obeys the diffusion equation,

$D_{0}\partial_{x}^{2}{\cal P}=\partial_{t}{\cal P}\,,$

and the boundary condition at the permeable membrane of permeability $\kappa$
is given by

$\kappa\left({\cal P}^{-}-{\cal P}^{+}\right)=-D_{0}\partial_{x}{\cal
P}^{-}\,,\quad D_{0}\partial_{x}{\cal P}^{-}=D_{0}\partial_{x}{\cal P}^{+}\,.$

Substituting into the Taylor expansion of Eq. (B.1) at $x_{M}$ and $t$ and
ignoring higher order terms, we obtain the permeation probability specifically
for the implementation of rejection sampling:

$\frac{P_{T}}{1-P_{T}}=\frac{\kappa\delta s}{D_{0}}\,,$

similar to the functional form in Eq. (C.1), except that $P_{T}$ is related to
the step size, rather than the particle’s distance to the membrane. Likewise,
the above expression of $P_{T}$ can be generalized in 2d and 3d by using Eqs.
(2a) and (C.7) in C.

It has been shown that both rejection sampling and elastic collision maintain
homogeneous particle density around membranes (Szymczak and Ladd, 2003).
Furthermore, to demonstrate the permeation probability in Eq. (3) for both
interactions, we performed MC simulations of diffusion between two parallel
planes separated by a distance $a=1$ $\mu$m in 1$d$, 2$d$, and 3$d$
(implemented by using rejection sampling or elastic collision). We created a
localized density source of Dirac delta function at time $t=0$, i.e.,
$\rho_{0}(x)=\delta_{D}(x)$, half-way between two permeable membranes of
permeability $=0.1$ $\mu$m/ms at $x_{M}=\pm a/2$. Then we calculated the
permeability $\kappa$ based on the particle density $\rho(x)$ around membranes
(Powles et al., 1992; Fieremans et al., 2008):

$-D_{0}^{+}\cdot\partial_{x}\,\rho(x_{M}^{+})=-D_{0}^{-}\cdot\partial_{x}\,\rho(x_{M}^{-})=\kappa\cdot\left[\rho(x_{M}^{-})-\rho(x_{M}^{+})\right]\,,$

where $D_{0}^{+}$ and $D_{0}^{-}$ were diffusivities on the right and left
side of the membrane, and $x_{M}^{+}$ and $x_{M}^{-}$ were positions at the
right and left side of the membrane. Here we applied $2\times 10^{7}$ random
walkers diffusing over 1000 steps with a step duration $\delta t=4\times
10^{-4}$ ms and a step size given by Eq. (1), where
$D_{0}=D_{0}^{+}=D_{0}^{-}=0.5$ $\mu$m2/ms.

For both rejection sampling and elastic collision, the permeability calculated
based on the particle density is consistent with the theoretical value input
in Eq. (3) (Fig. B.1). Ideally, the calculated permeability is independent of
diffusion time, and the spurious permeability time-dependence at short time
is, to the best of our understanding, caused by the discretization error due
to the small number of steps, when too few particles reach the membrane and
contribute to the flux.

Fig. B.1: Permeability calculated based on the particle density around the
membrane, which interacts with random walkers through rejection sampling (RS)
or elastic collision (EC). When a random walker encounters a membrane of
permeability $\kappa$, the permeation probability is determined by Eq. (3),
whose applicability for RS and EC is demonstrated by permeability calculated
based on the simulated particle density around two parallel planes (membranes)
separated by a distance $a=1$ $\mu$m in 1$d$, 2$d$ and 3$d$. To create a
density impulse at time $t=0$, random walkers are initialized at $x=0$ with
membranes at $x_{M}=\pm a/2$. The calculated permeability is consistent with
the theoretical value input in Eq. (3). To show the trend of data without
distorting the tendency, the curve of $\kappa$ over $t$ is smoothed by
applying a Savitzky-Golay filter of 11-point cubic polynomial.  Fig. B.2: The
bias of diffusivity parallel to membranes due to rejection sampling. The
algorithm of rejecting steps toward a membrane due to rejection sampling
causes the bias in diffusivity $D_{\parallel}$ parallel to the membrane (e.g.,
two parallel planes in 3d). This bias is roughly proportional to the step size
$\delta s$ and closely related to the surface-to-volume ratio $S/V$ of the
geometry.

###  B.2 Bias in the diffusivity parallel to the membrane caused by rejection
sampling

In this section, we will follow the framework in A.2 to evaluate the bias in
diffusivity parallel to the membrane due to rejection sampling. Considering
that the diffusion around a membrane in 3d obeys the rejection sampling, the
random walker cancels some steps (without updates) toward but not completely
transverse to the membrane, leading to smaller diffusion displacement cumulant
as well as diffusivity $D_{\parallel}(t)$ parallel to the membrane.

Given that a random walker close to a membrane with a distance $h\leq\delta s$
leaps into a direction of a polar angle $\theta$ and an azimuthal angle $\phi$
(spherical coordinate with zenith direction along $x$-axis), its diffusion
displacements parallel to the membrane (along, e.g., $z$-axis) is $\delta
z=\delta s\sin\theta\cos\phi$. Averaging over all possible directions allowed
by the rejection sampling, the second order cumulant of $\delta z$ is given by

$\displaystyle\langle\delta z^{2}\rangle_{\Omega}$
$\displaystyle=\int_{\cos^{-1}(\frac{h}{\delta s})}^{0}\int_{0}^{2\pi}(\delta
s\sin\theta\cos\phi)^{2}\,\frac{\sin\theta\,\text{d}\theta\,\text{d}\phi}{4\pi}\,$
$\displaystyle=\frac{\delta s^{2}}{12}\left[2+3\frac{h}{\delta
s}-\left(\frac{h}{\delta s}\right)^{3}\right]\,,\quad h\leq\delta s\,.$

Further averaging over the thickness $\delta s$ surrounding the membrane with
a homogeneous particle density maintained by rejection sampling, we obtain the
second order cumulant of $\delta z$, given by

$\langle\delta z^{2}\rangle_{h\leq\delta s}=\int_{0}^{\delta s}\langle\delta
z^{2}\rangle_{\Omega}\cdot p_{\Omega}\,\frac{\text{d}h}{\delta
s}=\frac{13}{48}\delta s^{2}\,,$

leading to a distinct diffusivity parallel to the membrane for walkers close
to the membrane:

$D_{0}^{\prime}=\frac{\langle\delta z^{2}\rangle_{h\leq\delta s}}{2\delta
t}=\frac{13}{16}D_{0}\,,$

where the step size in Eq. (1) in 3d is applied.

Random walkers surrounding the membrane ($h\leq\delta s$) have a diffusivity
$D_{0}^{\prime}$ parallel to the membrane and a volume fraction
$f^{\prime}=\frac{S}{V}\delta s$. Similarly, random walkers away from the
membrane have a diffusivity $D_{0}$ parallel to the membrane and a volume
fraction $f=1-f^{\prime}$. Their volume-weighted sum yields the overall
diffusivity parallel to the membrane, given by

$\displaystyle D_{\parallel}$ $\displaystyle=f\cdot D_{0}+f^{\prime}\cdot
D_{0}^{\prime}$ $\displaystyle\simeq
D_{0}\cdot\left(1-\frac{3}{16}\frac{S}{V}\delta s\right)<D_{0}\,,$ (B.2)

where the correction term is negligible when a small step size is applied.

The above Eq. (B.2) was demonstrated by performing MC simulations of diffusion
between two parallel planes in 3d (implemented by using rejection sampling),
as shown in Fig. B.2, where the step size $\delta s$ is tuned by varying the
time step $\delta t$.

##  Appendix C Unbiased simulation of membrane permeability

Here we introduce the theoretical background of unbiased simulations of
permeable membranes and provide a first order correction of permeation
probability by considering the particle density flux around membranes. The
theoretical results extend the applicability of related simulation models and
offer a guide to choose simulation parameters. In this section, the particle
density at position $x$ is denoted by $\rho(x)$.

###  C.1 The physics of the permeability correction: Equal molecular
concentration

The discussion in this section follows the logic of the Appendix B in our
previous work (Lee et al., 2020c). Given that a random walker encounters a
membrane of permeability $\kappa$, the permeation probability $P$ through the
membrane is related with the distance $h$ between the random walker and the
encountered membrane when $h\leq\delta s$, with $\delta s$ the step size in
Eq. (1), $D_{0}$ the intrinsic diffusivity and $\delta t$ the time step in $d$
dimensional space, as derived in Appendix A of (Fieremans et al., 2010), Eq.
(43):

$\frac{P}{1-P}=\frac{2\kappa h}{D_{0}}\,,$ (C.1)

which is a well-regularized functional form of $P$ even for a highly permeable
membrane: As expected, the limit $\kappa\to\infty$ corresponds to the
probability $P\to 1$.

In actual implementations, to reduce the computational load due to the
calculation of the distance $h$ from a random walker to the encountered
membrane, the distance $h$ could be approximated by the step size $\delta s$,
by averaging the particle density flux $j(h)$ around the membrane over
$h\leq\delta s$. To do so, a low probability is assumed ($P\ll 1$), such that
the denominator in the left-hand-side term in Eq. (C.1) is about $1$, i.e.,
$P(h)\simeq 2\kappa_{0}h/D_{0}$, where the permeability is redefined
($\kappa\to\kappa_{0}$) since the approximate relation is not exact. Averaging
the particle density flux $j(h)$ is equivalent to averaging the permeation
probability $P(h)$ because $j\propto P$, leading to

$P(h)\to P(\delta s)\simeq\frac{\int_{0}^{\delta
s}P(h)\,\tilde{p}_{\Omega}(h\leq\delta s)\,\text{d}h}{\int_{0}^{\delta
s}\tilde{p}_{\Omega}(h\leq\delta s)\,\text{d}h}=\frac{\kappa_{0}\delta
s}{D_{0}}\cdot C_{d}\ll 1\,,$

where $\tilde{p}_{\Omega}=1-p_{\Omega}$ is the fraction of directions
encountering the membrane, corresponding to the red solid angle in Fig. 1, and
$p_{\Omega}$ is defined in Eq. (A.1).

In this case, the permeation probability is given by Eq. (2a), where its
presumption, i.e., $P\ll 1$, yields a condition to be satisfied:

$\kappa_{0}\ll\sqrt{\frac{D_{0}}{2d\delta t}}\cdot\frac{1}{C_{d}}\,,$

indicating that using a short time step $\delta t$ is required in simulations
of exchange through a highly permeable membrane (large $\kappa_{0}$); in this
case, the genuine permeability $\kappa$ is close to the input value
$\kappa_{0}$.

For larger $\kappa$, the input $\kappa_{0}$ would be significantly different
from the genuine value $\kappa$ in simulations while extending the
approximation of $h$ by $\delta s$. We will show that averaging over
$h\leq\delta s$ simply renormalizes the input $\kappa_{0}$ entering Eq. (2a),
leading to a genuine $\kappa$ in Eq. (C.7) given below. The idea behind is
that averaging over $\delta s$ affects not just the permeation probability but
also the particle flux density. We here demand the Fick’s first law satisfying
the permeation probability in Eq. (2a) and derive a correction factor
renormalizing the permeability $\kappa_{0}\to\kappa$.

The particle density flux from left to right ($1\to 2$) is given by

$j_{1\to 2}\simeq\frac{(2dC_{d})^{-1}\cdot\langle\rho_{1}\rangle\cdot(S\delta
s_{1})\cdot P_{1\to 2}}{S\cdot\delta t}\,,$ (C.2)

where $S$ is the surface area, $P_{1\to 2}$ is the permeation probability from
left to right given by Eq. (2a), $\langle\rho_{1}\rangle$ is the particle
density averaged over the layer (of thickness $\delta s_{1}$) on the left side
of the membrane, and $(2dC_{d})^{-1}$ is the mean fraction of the hopping
orientation (averaged over $\delta s_{1}$ as well) along which particles
encounter the membrane (Fieremans and Lee, 2018).

Substituting Eq. (2a) into Eq. (C.2) and using $(\delta s_{1})^{2}/\delta
t=2dD_{1}$, we obtain $j_{1\to 2}\simeq\kappa_{0}\cdot\langle\rho_{1}\rangle$.
Similarly, the particle density flux from right to left side is $j_{2\to
1}\simeq-\kappa_{0}\cdot\langle\rho_{2}\rangle$. Then the net particle flux
density is given by

$j=j_{1\to 2}+j_{2\to
1}\simeq\kappa_{0}\cdot\left(\langle\rho_{1}\rangle-\langle\rho_{2}\rangle\right)\,.$
(C.3)

Given that the particle density right at left and right surface of the
membrane is $\rho_{0,1}$ and $\rho_{0,2}$ (without spatial averaging), the net
particle density flux is (by definition of the genuine $\kappa$)

$\displaystyle j$ $\displaystyle=\kappa\cdot(\rho_{0,1}-\rho_{0,2})$ (C.4)
$\displaystyle=-D_{1}\cdot\partial_{x}\,\rho_{0,1}=-D_{2}\cdot\partial_{x}\,\rho_{0,2}\,,$
(C.5)

where $\kappa$ is the genuine permeability that we would like to achieve with
simulations, different from the input value $\kappa_{0}$, and
$\partial_{x}\,\rho_{0,1}$ and $\partial_{x}\,\rho_{0,2}$ are density
gradients right at left and right surface of the membrane.

Here we average the density $\langle\rho_{1}\rangle$ (and
$\langle\rho_{2}\rangle$) over the layer on the left (and right) side of the
membrane, of thickness $\delta s_{1}$ (and $\delta s_{2}$), and equate the
flux density in Eq. (C.3) to that in Eq. (C.4) to obtain the genuine
permeability $\kappa$.

Approximating the particle density ($\rho_{1}$, $\rho_{2}$) variation close to
the membrane with a linear function of the distance from the membrane, we have

$\rho_{1}(x)\simeq\rho_{0,1}+\partial_{x}\,\rho_{0,1}\cdot
x\,,\quad\rho_{2}(x)\simeq\rho_{0,2}+\partial_{x}\,\rho_{0,2}\cdot x\,.$

Considering the fraction $\tilde{p}_{\Omega}(x)$ of the hopping orientation
along which particles encounter the membrane, as shown in Fig. A1 in
(Fieremans and Lee, 2018), the particle density averaged within the thickness
of step size is given by

$\displaystyle\langle\rho_{1}\rangle$ $\displaystyle=\frac{\int_{-\delta
s_{1}}^{0}\rho_{1}(x)\tilde{p}_{\Omega}(x)\,\text{d}x}{\int_{-\delta
s_{1}}^{0}\tilde{p}_{\Omega}(x)\,\text{d}x}=\rho_{0,1}-\partial_{x}\,\rho_{0,1}\cdot\delta
s_{1}\cdot\frac{C_{d}}{2}\,,$ (C.6a) $\displaystyle\langle\rho_{2}\rangle$
$\displaystyle=\frac{\int_{0}^{\delta
s_{2}}\rho_{2}(x)\tilde{p}_{\Omega}(x)\,\text{d}x}{\int_{0}^{\delta
s_{2}}\tilde{p}_{\Omega}(x)\,\text{d}x}=\rho_{0,2}+\partial_{x}\,\rho_{0,2}\cdot\delta
s_{2}\cdot\frac{C_{d}}{2}\,,$ (C.6b)

where

$\tilde{p}_{\Omega}(x)=\begin{cases}1-p_{\Omega}(|x|\leq\delta
s_{1})&x<0\,,\\\ 1-p_{\Omega}(|x|\leq\delta s_{2})&x\geq 0\,,\end{cases}$

and $p_{\Omega}$ is defined in Eq. (A.1).

Substituting Eqs. (C.3) and (C.5)–(C.6) into Eq. (C.4) yields

$\kappa=\frac{\kappa_{0}}{1-\alpha}>\kappa_{0}\,,$ (C.7)

where

$\displaystyle\alpha$ $\displaystyle=\frac{1}{2}\kappa_{0}\left(\frac{\delta
s_{1}}{D_{1}}+\frac{\delta s_{2}}{D_{2}}\right)\cdot C_{d}$
$\displaystyle=\frac{P_{1\to 2}+P_{2\to 1}}{2}\,.$ (C.8)

Interestingly, the correction factor $\alpha$ is the permeation probability
averaged for both directions, i.e., $\alpha\in[0,1]$. Therefore, the genuine
permeability $\kappa$ in Monte Carlo simulations of any dimension $d$ is
always larger than the input value $\kappa_{0}$, as in Eq. (C.7), where the
correction factor $\alpha$ is essential especially when simulating the
diffusion across a highly permeable membrane. To minimize $\alpha$ and reduce
the bias, a smaller time-step and larger intrinsic diffusivity should be used.

Practically, to simulate a membrane of permeability $\kappa$, we have to tune
the input permeability $\kappa_{0}$ for the permeation permeability in Eq.
(2a) based on

$\kappa_{0}=\frac{\kappa}{1+\kappa\cdot(\alpha/\kappa_{0})}\,,$

where the right-hand side is independent of $\kappa_{0}$ due to Eq. (C.8).
Substituting the above relation into Eq. (2a) yields the corrected permeation
probability in Eq. (3). The above correction ensures the genuine permeability
$\kappa$ in simulations.

Furthermore, the corrected permeation probability in Eq. (3) should still be
$\ll 1$, leading to the following constraint, as a guidance of choosing
simulation parameters:

$\kappa\ll\sqrt{\frac{2}{d\delta
t}}\cdot\frac{1}{C_{d}}\cdot\frac{\sqrt{D_{1}D_{2}}}{\left|\sqrt{D_{1}}-\sqrt{D_{2}}\right|}\,.$

In other words, Eq. (3) works particularly well for a small time-step $\delta
t$, large intrinsic diffusivities $\sqrt{D_{1}D_{2}}$, and similar intrinsic
diffusivities between compartments ($D_{1}\simeq D_{2}$).

Fig. C.1: Non-zero density flux close to the membrane leads to the bias in
membrane permeability. The particle density across a permeable membrane has an
offset due to the finite permeability. The input permeability $\kappa_{0}$ is
related to the difference of particle density ($\langle\rho_{1}\rangle$,
$\langle\rho_{2}\rangle$) averaged within the thickness of the hopping step
size ($\delta s_{1}$, $\delta s_{2}$); however, the genuine permeability
$\kappa$, to be achieved in simulations, is related to the difference of the
particle density ($\rho_{0,1}$, $\rho_{0,2}$) at the membrane. As a result,
the genuine permeability $\kappa$ is always larger than the input value
$\kappa_{0}$. To clarify, only density flux to the right is plotted in this
figure. In actual simulations, the balance of the density flux in both
directions yields a homogeneous particle density in each compartment.

###  C.2 General case: Different spin concentration at both sides of the
membrane

In the previous section, the medium is assumed to have the same spin
concentration in all compartments. However, a lower spin concentration is
expected for some tissue microstructure, such as myelin water. To generalize
for different spin concentrations in each compartment, the permeation
probability in Eq. (2a) is re-written as

$\displaystyle P_{1\to 2}$ $\displaystyle\simeq\frac{\kappa_{0}\delta
s_{1}}{D_{1}}\cdot C_{d}\cdot\left(\frac{c_{2}}{c_{1}}\right)^{\lambda}\,,$
(C.9a) $\displaystyle P_{2\to 1}$ $\displaystyle\simeq\frac{\kappa_{0}\delta
s_{2}}{D_{2}}\cdot C_{d}\cdot\left(\frac{c_{2}}{c_{1}}\right)^{\lambda-1}\,,$
(C.9b)

where $c_{1}$ and $c_{2}$ are spin concentrations over the left and right
sides of the membrane, and $\lambda\in[0,1]$ is an exponent determined later.
It is worthwhile to notice that the probability ratio $P_{1\to 2}/P_{2\to
1}=c_{2}\sqrt{D_{2}}/c_{1}\sqrt{D_{1}}$ is maintained to ensure the particle
density equilibrium for all diffusion times.

Similar to the derivation in previous section, substituting Eq. (C.9) into Eq.
(C.2) and calculating $j_{1\to 2}$ and $j_{2\to 1}$, we obtain

$\displaystyle j$ $\displaystyle=j_{1\to 2}+j_{2\to 1}$
$\displaystyle\simeq\kappa_{0}\cdot\left[\langle\rho_{1}\rangle\cdot\left(\frac{c_{2}}{c_{1}}\right)^{\lambda}-\langle\rho_{2}\rangle\cdot\left(\frac{c_{2}}{c_{1}}\right)^{\lambda-1}\right]\,.$
(C.10)

Considering the ratio $c_{2}/c_{1}$ of spin concentrations, the net particle
flux density is given by

$\displaystyle j$
$\displaystyle=\kappa\cdot\left[\rho_{0,1}\cdot\left(\frac{c_{2}}{c_{1}}\right)^{\lambda}-\rho_{0,2}\cdot\left(\frac{c_{2}}{c_{1}}\right)^{\lambda-1}\right]$
(C.11)
$\displaystyle=-D_{1}\partial_{x}\,\rho_{0,1}=-D_{2}\partial_{x}\,\rho_{0,2}\,,$
(C.12)

where the unbiased permeability $\kappa$ is re-defined accordingly.

Substituting Eqs. (C.6), (C.10) and (C.12) into Eq. (C.11) yields

$\kappa=\frac{\kappa_{0}}{1-\alpha_{\lambda}}\,,$ (C.13)

where

$\displaystyle\alpha_{\lambda}$
$\displaystyle=\frac{1}{2}\kappa_{0}\left[\frac{\delta
s_{1}}{D_{1}}\cdot\left(\frac{c_{2}}{c_{1}}\right)^{\lambda}+\frac{\delta
s_{2}}{D_{2}}\cdot\left(\frac{c_{2}}{c_{1}}\right)^{\lambda-1}\right]\cdot
C_{d}$ $\displaystyle=\frac{P_{1\to 2}+P_{2\to 1}}{2}\,.$ (C.14)

The choice of $\lambda$ is essential not only for the generalization of
permeability definition, as in Eqs. (C.9)–(C.11), but also for the
permeability bias in Eq. (C.13). On the one hand, to maintain the same
permeability definition for all membranes in the medium, we can fix $\lambda$
as a constant, e.g., 0, 1/2 or 1. On the other hand, to minimize the
correction factor $\alpha_{\lambda}$, we can chose the $\lambda$ based on

$\displaystyle\lambda$ $\displaystyle=\underset{\lambda\in[0,1]}{\text{arg
min}}(\alpha_{\lambda})$
$\displaystyle=\left(1+\frac{c_{2}\sqrt{D_{2}}}{c_{1}\sqrt{D_{1}}}\right)^{-1}\,,$

which is well-regularized for even extreme cases, e.g., $c_{1}=0$ or
$c_{2}=0$.

###  C.3 Alternative approach for simulations of permeable membrane

Instead of assigning a nominal permeability $\kappa$ for a permeable membrane,
Baxter and Frank (2013) defined the permeation probability based on the spin
concentration ($c_{1}$, $c_{2}$) and intrinsic diffusivity ($D_{1}$, $D_{2}$)
over the left and right side of the membrane:

$\displaystyle P_{1\to 2}$
$\displaystyle=\frac{c_{2}\sqrt{D_{2}}}{c_{1}\sqrt{D_{1}}}\,,$ (C.15a)
$\displaystyle P_{2\to 1}$ $\displaystyle=1\,,$ (C.15b)

where the left side compartment 1 is a “high-flux medium”, compared with the
right side compartment 2, i.e., $c_{1}\sqrt{D_{1}}\geq c_{2}\sqrt{D_{2}}$.
This approach has been applied in simulations of, e.g., the exchange between
intra-/extra-axonal water an myelin water (Harkins and Does, 2016) and the
water exchange between intra-axonal cytoplasm and mitochondria (Lee et al.,
2020b). It seems that this method introduces neither adjustable parameters for
membrane permeability nor particle density transition over the membrane;
however, this is true only for infinitely small time-step $\delta t$. For
finite $\delta t$, a $\delta t$-dependent permeability may emerge in
simulations.

The derivation of this extra permeability is similar to those in previous
sections. Substituting Eqs. (C.6) and (C.15) into Eq. (C.2), the particle
density flux for both directions is given by

$\displaystyle j_{1\to 2}$
$\displaystyle\simeq\kappa_{0}^{\prime}\cdot\left(\rho_{0,1}-\partial_{x}\,\rho_{0,1}\cdot\delta
s_{1}\cdot\frac{C_{d}}{2}\right)\cdot\frac{c_{2}}{c_{1}}\,,$ $\displaystyle
j_{2\to 1}$
$\displaystyle\simeq-\kappa_{0}^{\prime}\cdot\left(\rho_{0,2}+\partial_{x}\,\rho_{0,2}\cdot\delta
s_{2}\cdot\frac{C_{d}}{2}\right)\,,$

where

$\kappa_{0}^{\prime}=\sqrt{\frac{D_{2}}{2d\cdot\delta
t}}\cdot\frac{1}{C_{d}}.$

Therefore, the net particle density flux is

$\displaystyle j$ $\displaystyle=j_{1\to 2}+j_{2\to 1}$
$\displaystyle\simeq\kappa_{0}^{\prime}\cdot\left[\left(\rho_{0,1}\cdot\frac{c_{2}}{c_{1}}-\rho_{0,2}\right)-\frac{C_{d}}{2}\left(\partial_{x}\,\rho_{0,1}\cdot\delta
s_{1}\cdot\frac{c_{2}}{c_{1}}+\partial_{x}\,\rho_{0,2}\cdot\delta
s_{2}\right)\right]\,,$ (C.16)

which indicates the equilibrium condition at $t\to\infty$ limit:

$j\to 0\,,\quad\partial_{x}\,\rho_{0,1}\to
0\,,\quad\partial_{x}\,\rho_{0,2}\to
0\,,\quad\frac{\rho_{0,1}}{\rho_{0,2}}\to\frac{c_{1}}{c_{2}}\,.$

Considering the ratio of spin concentrations, we have the net particle density
flux

$\displaystyle j$
$\displaystyle=\kappa^{\prime}\cdot\left(\rho_{0,1}\cdot\frac{c_{2}}{c_{1}}-\rho_{0,2}\right)$
(C.17)
$\displaystyle=-D_{1}\cdot\partial_{x}\,\rho_{0,1}=-D_{2}\cdot\partial_{2}\,\rho_{0,2}\,,$
(C.18)

where $\kappa^{\prime}$ is the effective permeability due to the finite time-
step. Substituting Eqs. (C.16) and (C.18) into Eq. (C.17) yields

$\kappa^{\prime}=\frac{\kappa_{0}^{\prime}}{1-\beta}\propto\frac{1}{\sqrt{\delta
t}}\,,$ (C.19)

where

$\displaystyle\beta$
$\displaystyle=\frac{1}{2}\left(\frac{c_{2}\sqrt{D_{2}}}{c_{1}\sqrt{D_{1}}}+1\right)$
$\displaystyle=\frac{P_{1\to 2}+P_{2\to 1}}{2}\,.$

Interestingly, similar to $\alpha$ in Eq. (C.8) and $\alpha_{\lambda}$ in Eq.
(C.14), the correction factor $\beta\in[0.5,1]$ is also the permeation
probability averaged for both directions.

For an infinitely small time-step ($\delta t\to 0$), the effective
permeability is infinitely large ($\kappa^{\prime}\to\infty$), as predicted by
Eq. (C.19). In this case, the finite density flux $j$ in Eq. (C.17) indicates
no particle density transition over the membrane, i.e.
$\rho_{0,1}/\rho_{0,2}\to c_{1}/c_{2}$. Similarly, when
$c_{1}\sqrt{D_{1}}/c_{2}\sqrt{D_{2}}=1$ in Eq. (C.15), the $P_{1\to 2}=1$
leads to an infinitely large $\kappa^{\prime}$, based on Eq. (C.19), and
$\rho_{1}/\rho_{2}\to c_{1}/c_{2}$ due to the finite density flux $j$ in Eq.
(C.17).

In contrast, when $P_{1\to 2}<1$, a finite time-step $\delta t$ results in an
extra effective permeability $\kappa^{\prime}$, hindering the permeation
through membranes. To reduce this unwanted effect, the applied time step needs
to be small. For example, considering a multi-compartmental system in $1d$,
the size and intrinsic diffusivity in the $i$-th compartment are $l_{i}$ and
$D_{i}$. Then we can ignore the hindrance through membranes caused by
$\kappa^{\prime}$, if the time-step is sufficiently small, such that

$\kappa^{\prime}\gg\text{max}\left(\frac{D_{i}}{l_{i}}\right)\,,$

where $D_{i}/l_{i}$ is the intrinsic permeability of the $i$-th compartment
(Novikov et al., 2011). However, in $2d$ and $3d$, the compartment length
scale $l_{i}$ could be ill-defined, complicating the choice of time step
$\delta t$.

##  Appendix D Estimation of the undulation amplitude and wavelength from a 3d
axonal skeleton

In this Appendix, we will introduce how to estimate the length scale of
undulation amplitude and wavelength based on a 3d axonal skeleton.

Considering an axonal skeleton aligned to its main axis (z-axis), the skeleton
can be quantified as ${\bf l}={\bf w}+z\hat{\bf z}$, with $\bf w$ the
deviation of the skeleton from the main axis. We can decompose the axonal
skeleton with multiple harmonics (Lee et al., 2020a):

${\bf w}=\sum_{n}w_{xn}\cos(k_{n}z+\phi_{xn})\,\hat{\bf
x}+w_{yn}\cos(k_{n}z+\phi_{yn})\,\hat{\bf y}\,,$

where, for the $n$-th harmonic ($n\in\mathbb{N}$), $w_{xn}$ and $w_{yn}$ are
the undulation amplitudes along x- and y-axes, and $\phi_{xn}$ and $\phi_{yn}$
are the phases. Here we focus on the case of $k_{n}=n\cdot 2\pi/L_{z}$ (with
$L_{z}$ the axonal length along $z$) such that all harmonics are orthogonal.
This allows us to define a Euclidean measure of the undulation amplitude,
given by

$w_{0}\equiv\sqrt{\sum_{n}w_{0n}^{2}}\,,\quad
w_{0n}\equiv\sqrt{w_{xn}^{2}+w_{yn}^{2}}\,.$

The values of $w_{xn}$ and $w_{yn}$ can be estimated by applying discrete
cosine transform to $\bf w$, yielding a length scale of the undulation
amplitude $\sim w_{0}$.

Furthermore, the projection factor transverse to the main direction has been
shown to be (Lee et al., 2020a)

$\langle\sin^{2}\delta\theta\rangle_{z}\simeq\frac{2}{\xi^{2}}\sum_{n}\left(\frac{\pi
w_{0n}}{\lambda_{n}}\right)^{2}\,,\quad\xi\simeq 1+\sum_{n}\left(\frac{\pi
w_{0n}}{\lambda_{n}}\right)^{2}\gtrsim 1\,,$ (D.1)

where $\delta\theta=\delta\theta(z)$ is the angle between the individual
axon’s skeleton segment at $z$ and its main axis, $\langle...\rangle_{z}$ is
the average along each axon’s main axis, and $\lambda_{n}=2\pi/k_{n}$ is the
corresponding wavelength. To have a rough estimate of the undulation
wavelength, we impose a simplified 1-harmonic model ($n\leq 1$) to Eq. (D.1),
leading to

$\lambda\simeq\pi
w_{0}\sqrt{\frac{2}{\langle\sin^{2}\delta\theta\rangle_{z}}}\,,$ (D.2)

where we approximate $\xi\simeq 1$ due to the small $w_{0}/\lambda\sim 0.02$
in realistic axons of a mouse brain (Lee et al., 2020a).

Note that Eq. (D.2) should not generally be interpreted as a proportionality
relation between the undulation wavelength and its amplitude, since each axon
can have its own dispersion $\delta\theta(z)$. However, empirically, this
proportionality seems to be valid based on the high correlation between these
two geometric quantities in Fig. 3d. In other words, the slope
$\text{d}w_{0}/\text{d}\lambda\simeq\langle\sin^{2}\delta\theta\rangle_{z}^{1/2}/\left(\pi\sqrt{2}\right)$
seems to be sufficiently axon-independent. From the estimated slope in that
plot, we find the intra-axon undulation dispersion
$\sqrt{\langle\delta\theta^{2}\rangle}\simeq\langle\sin^{2}\delta\theta\rangle_{z}^{1/2}\approx
10^{\circ}$, which is about a factor of 2 smaller than the inter-axon fiber
orientation dispersion $\sqrt{\langle\bar{\theta}^{2}\rangle}\simeq
22^{\circ}$ estimated from histology (Ronen et al., 2014; Lee et al., 2019).
By using the Rodrigues’ rotation formula and the small angle approximation, it
is straightforward to combine the two contributions and estimate the overall
dispersion angle
$\sqrt{\langle\theta^{2}\rangle}\simeq\sqrt{\langle\bar{\theta}^{2}\rangle+\langle\delta\theta^{2}\rangle}\approx
24^{\circ}$. This means that the “Standard Model” estimates of dispersion
(Novikov et al., 2018; Dhital et al., 2019) is dominated by the inter-axonal
contribution, in agreement with the Standard Model assumptions.

## References

  * Adams and Bischof (1994) Adams, R., Bischof, L., 1994\.  Seeded region growing.  IEEE Transactions on pattern analysis and machine intelligence 16, 641–647.
  * Alexander et al. (2018) Alexander, D.C., Dyrby, T.B., Nilsson, M., Zhang, H., 2018\.  Imaging brain microstructure with diffusion MRI: practicality and applications.  NMR in Biomedicine doi:10.1002/NBM.3841.
  * Baxter and Frank (2013) Baxter, G.T., Frank, L.R., 2013\.  A computational model for diffusion weighted imaging of myelinated white matter.  NeuroImage 75, 204–212.
  * Brabec et al. (2020) Brabec, J., Lasič, S., Nilsson, M., 2020.  Time-dependent diffusion in undulating thin fibers: Impact on axon diameter estimation.  NMR in Biomedicine 33, e4187.
  * Burcaw et al. (2015) Burcaw, L.M., Fieremans, E., Novikov, D.S., 2015.  Mesoscopic structure of neuronal tracts from time-dependent diffusion.  NeuroImage 114, 18–37.
  * Callaghan (1993) Callaghan, P.T., 1993.  Principles of nuclear magnetic resonance microscopy.  Oxford University Press on Demand.
  * Callaghan et al. (2020) Callaghan, R., Alexander, D.C., Palombo, M., Zhang, H., 2020\.  ConFiG: Contextual fibre growth to generate realistic axonal packing for diffusion MRI simulation.  NeuroImage 220, 117107\.
  * Chin et al. (2002) Chin, C.L., Wehrli, F.W., Hwang, S.N., Takahashi, M., Hackney, D.B., 2002.  Biexponential diffusion attenuation in the rat spinal cord: computer simulations based on anatomic images of axonal architecture.  Magnetic Resonance in Medicine 47, 455–460.
  * Dhital et al. (2019) Dhital, B., Reisert, M., Kellner, E., Kiselev, V.G., 2019\.  Intra-axonal diffusivity in brain white matter.  NeuroImage 189, 543–550.
  * Einstein et al. (1905) Einstein, A., et al., 1905.  On the motion of small particles suspended in liquids at rest required by the molecular-kinetic theory of heat.  Annalen der physik 17, 208\.
  * Fieremans et al. (2016) Fieremans, E., Burcaw, L.M., Lee, H.H., Lemberskiy, G., Veraart, J., Novikov, D.S., 2016\.  In vivo observation and biophysical interpretation of time-dependent diffusion in human white matter.  NeuroImage 129, 414–427.
  * Fieremans et al. (2008) Fieremans, E., De Deene, Y., Delputte, S., Özdemir, M.S., D’Asseler, Y., Vlassenbroeck, J., Deblaere, K., Achten, E., Lemahieu, I., 2008.  Simulation and experimental verification of the diffusion in an anisotropic fiber phantom.  Journal of magnetic resonance 190, 189–199.
  * Fieremans and Lee (2018) Fieremans, E., Lee, H.H., 2018\.  Physical and numerical phantoms for the validation of brain microstructural MRI: A cookbook.  NeuroImage 182, 39–61.
  * Fieremans et al. (2010) Fieremans, E., Novikov, D.S., Jensen, J.H., Helpern, J.A., 2010\.  Monte Carlo study of a two-compartment exchange model of diffusion.  NMR in Biomedicine 23, 711–724.
  * Ford and Hackney (1997) Ford, J.C., Hackney, D.B., 1997\.  Numerical model for calculation of apparent diffusion coefficients (adc) in permeable cylinders—comparison with measured adc in spinal cord white matter.  Magnetic Resonance in Medicine 37, 387–394.
  * Ginsburger et al. (2019) Ginsburger, K., Matuschke, F., Poupon, F., Mangin, J.F., Axer, M., Poupon, C., 2019\.  MEDUSA: A GPU-based tool to create realistic phantoms of the brain microstructure using tiny spheres.  NeuroImage 193, 10–24.
  * Grebenkov (2007) Grebenkov, D.S., 2007.  NMR survey of reflected brownian motion.  Reviews of Modern Physics 79, 1077\.
  * Hall and Alexander (2009) Hall, M.G., Alexander, D.C., 2009\.  Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI.  IEEE transactions on medical imaging 28, 1354–1364.
  * Harkins and Does (2016) Harkins, K.D., Does, M.D., 2016\.  Simulations on the influence of myelin water in diffusion-weighted imaging.  Physics in Medicine & Biology 61, 4729.
  * Jelescu and Budde (2017) Jelescu, I.O., Budde, M.D., 2017\.  Design and Validation of Diffusion MRI Models of White Matter.  Frontiers in Physics 5.  doi:10.3389/fphy.2017.00061.
  * Jensen et al. (2016) Jensen, J.H., Glenn, G.R., Helpern, J.A., 2016.  Fiber ball imaging.  NeuroImage 124, 824–833.
  * Jensen et al. (2005) Jensen, J.H., Helpern, J.A., Ramani, A., Lu, H., Kaczynski, K., 2005.  Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging.  Magnetic Resonance in Medicine 53, 1432–1440.
  * Jóhannesson and Halle (1996) Jóhannesson, H., Halle, B., 1996\.  Solvent diffusion in ordered macrofluids: a stochastic simulation study of the obstruction effect.  The Journal of chemical physics 104, 6807–6817.
  * Jones (2010) Jones, D.K., 2010.  Diffusion MRI.  Oxford University Press.
  * Kiselev (2017) Kiselev, V.G., 2017.  Fundamentals of diffusion MRI physics.  NMR in Biomedicine 30, e3602.
  * Lee et al. (2020a) Lee, H.H., Jespersen, S.N., Fieremans, E., Novikov, D.S., 2020a.  The impact of realistic axonal shape on axon diameter estimation using diffusion MRI.  NeuroImage , 117228.
  * Lee et al. (2020b) Lee, H.H., Papaioannou, A., Kim, S.L., Novikov, D.S., Fieremans, E., 2020b.  A time-dependent diffusion MRI signature of axon caliber variations and beading.  Communications biology 3, 1–13.
  * Lee et al. (2020c) Lee, H.H., Papaioannou, A., Novikov, D.S., Fieremans, E., 2020c.  In vivo observation and biophysical interpretation of time-dependent diffusion in human cortical gray matter.  NeuroImage , 117054.
  * Lee et al. (2020d) Lee, H.H., Tian, Q., Ngamsombat, C., Berger, D.R., Lichtman, J.W., Huang, S.Y., Novikov, D.S., Fieremans, E., 2020d.  Random walk simulations of diffusion in human brain white matter from 3d EM validate diffusion time-dependence transverse and parallel to axons.  28th Annual Meeting of the International Society for Magnetic Resonance in Medicine, Paris, France 28.
  * Lee et al. (2019) Lee, H.H., Yaros, K., Veraart, J., Pathan, J.L., Liang, F.X., Kim, S.G., Novikov, D.S., Fieremans, E., 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 Structure and Function 224, 1469–1488.
  * McNab and Miller (2008) McNab, J.A., Miller, K.L., 2008\.  Sensitivity of diffusion weighted steady state free precession to anisotropic diffusion.  Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 60, 405–413.
  * Neuman (1974) Neuman, C., 1974.  Spin echo of spins diffusing in a bounded medium.  The Journal of Chemical Physics 60, 4508–4511.
  * Nguyen et al. (2018) Nguyen, K.V., Hernández-Garzón, E., Valette, J., 2018.  Efficient GPU-based Monte-Carlo simulation of diffusion in real astrocytes reconstructed from confocal microscopy.  Journal of Magnetic Resonance 296, 188–199.
  * Nilsson et al. (2012) Nilsson, M., Lätt, J., Ståhlberg, F., van Westen, D., Hagslätt, H., 2012.  The importance of axonal undulation in diffusion MR measurements: a Monte Carlo simulation study.  NMR in Biomedicine 25, 795–805.
  * Novikov (2020) Novikov, D.S., 2020.  The present and the future of microstructure mri: From a paradigm shift to “normal science”.  Journal of Neuroscience Methods , 108947\.
  * Novikov et al. (2011) Novikov, D.S., Fieremans, E., Jensen, J.H., Helpern, J.A., 2011\.  Random walk with barriers.  Nature Physics 7, 508–514.
  * Novikov et al. (2019) Novikov, D.S., Fieremans, E., Jespersen, S.N., Kiselev, V.G., 2019\.  Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation.  NMR in Biomedicine 32, e3998.
  * Novikov et al. (2014) Novikov, D.S., Jensen, J.H., Helpern, J.A., Fieremans, E., 2014\.  Revealing mesoscopic structural universality with diffusion.  Proceedings of the National Academy of Sciences of the United States of America 111, 5088–5093.
  * Novikov et al. (2018) Novikov, D.S., Kiselev, V.G., Jespersen, S.N., 2018.  On modeling.  Magnetic Resonance in Medicine 79, 3172–3193.  doi:10.1002/mrm.27101.
  * Novikov et al. (2018) Novikov, D.S., Veraart, J., Jelescu, I.O., Fieremans, E., 2018\.  Rotationally-invariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion MRI.  NeuroImage 174, 518–538.
  * Palombo et al. (2019) Palombo, M., Alexander, D.C., Zhang, H., 2019.  A generative model of realistic brain cells with application to numerical simulation of the diffusion-weighted MR signal.  NeuroImage 188, 391–402.
  * Palombo et al. (2018) Palombo, M., Ligneul, C., Hernandez-Garzon, E., Valette, J., 2018\.  Can we detect the effect of spines and leaflets on the diffusion of brain intracellular metabolites?  NeuroImage 182, 283–293.
  * Powles et al. (1992) Powles, J.G., Mallett, M., Rickayzen, G., Evans, W., 1992\.  Exact analytic solutions for diffusion impeded by an infinite array of partially permeable barriers.  Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 436, 391–403.
  * Ronen et al. (2014) Ronen, I., Budde, M., Ercan, E., Annese, J., Techawiboonwong, A., Webb, A.G., 2014\.  Microstructural organization of axons in the human corpus callosum quantified by diffusion-weighted magnetic resonance spectroscopy of n-acetylaspartate and post-mortem histology.  Brain Structure & Function 219, 1773–1785.
  * Slijkerman and Hofman (1998) Slijkerman, W., Hofman, J., 1998\.  Determination of surface relaxivity from NMR diffusion measurements.  Magnetic resonance imaging 16, 541–544.
  * Stejskal and Tanner (1965) Stejskal, E.O., Tanner, J.E., 1965\.  Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient.  The journal of chemical physics 42, 288–292.
  * Szafer et al. (1995) Szafer, A., Zhong, J., Gore, J.C., 1995.  Theoretical model for water diffusion in tissues.  Magnetic Resonance in Medicine 33, 697–712.
  * Szymczak and Ladd (2003) Szymczak, P., Ladd, A., 2003\.  Boundary conditions for stochastic solutions of the convection-diffusion equation.  Physical Review E 68, 036704\.
  * Tanner (1970) Tanner, J.E., 1970.  Use of the stimulated echo in NMR diffusion studies.  The Journal of Chemical Physics 52, 2523–2526.
  * Topgaard (2017) Topgaard, D., 2017.  Multidimensional diffusion MRI.  Journal of Magnetic Resonance 275, 98–113.  doi:10.1016/j.jmr.2016.12.007.
  * Veraart et al. (2019) Veraart, J., Fieremans, E., Novikov, D.S., 2019.  On the scaling behavior of water diffusion in human brain white matter.  NeuroImage 185, 379–387.
  * Veraart et al. (2020) Veraart, J., Nunes, D., Rudrapatna, U., Fieremans, E., Jones, D.K., Novikov, D.S., Shemesh, N., 2020.  Noninvasive quantification of axon radii using diffusion MRI.  Elife 9, e49855.
  * Waudby and Christodoulou (2011) Waudby, C.A., Christodoulou, J., 2011\.  GPU accelerated Monte Carlo simulation of pulsed-field gradient NMR experiments.  Journal of Magnetic Resonance 211, 67–73.
  * West et al. (2016) West, K.L., Kelm, N.D., Carson, R.P., Does, M.D., 2016\.  A revised model for estimating g-ratio from MRI.  NeuroImage 125, 1155–1158.
  * Woessner (1961) Woessner, D.E., 1961.  Effects of diffusion in nuclear magnetic resonance spin-echo experiments.  The Journal of Chemical Physics 34, 2057–2061.
  * Xing et al. (2013) Xing, H., Lin, F., Wu, Q., Gong, Q., 2013\.  Investigation of different boundary treatment methods in Monte‐Carlo simulations of diffusion NMR.  Magnetic Resonance in Medicine 70, 1167–1172.
  * Xu et al. (2018) Xu, T., Foxley, S., Kleinnijenhuis, M., Chen, W.C., Miller, K.L., 2018.  The effect of realistic geometries on the susceptibility-weighted MR signal in white matter.  Magnetic Resonance in Medicine 79, 489–500.
  * Yeh et al. (2013) Yeh, C.H., Schmitt, B., Le Bihan, D., Li-Schlittgen, J.R., Lin, C.P., Poupon, C., 2013\.  Diffusion microscopist simulator: a general Monte Carlo simulation system for diffusion magnetic resonance imaging.  PloS One 8, e76626.