Monte Carlo Simulator for Diffusion-weighted Imaging Sequences.

We developed a Monte Carlo simulator for diffusion-weighted imaging sequences which displays the motion of water molecules and computes the dynamic phase dispersion due to the applied motion probing gradients. This simulator can be used to validate the analytical equations of diffusion models and understand their limitations due to their approximations. Here, we introduce the software and some specific use cases. The software can be downloaded from the following website: https://www.nirs.qst.go.jp/amr_diag.

Introduction

In vivo exploration of the microstructures of the central nervous system by using diffusion-weighted images has recently become both more popular and more complicated by the development of the advanced imaging techniques. For examples, multi-diffusion-time encoding[1] in addition to multiple b-value encoding,[2] and double diffusion encoding (DDE).[3] These new techniques enable us to find out more about the detailed structure of in vivo tissues than before, but it can be difficult to intuitively understand which diffusional phenomenon is being visualized and how it reflects the microstructure of the tissue. Therefore, we developed a Monte Carlo simulation software to help validate and understand these imaging techniques. The source code works on Matlab R2017b equipped with a parallel computing toolbox (Mathworks, Natick, MA, USA). It is available at our website (https://www.nirs.qst.go.jp/amr_diag). In addition, the software has been combined with a 2D axon packing software[4] and included to the qMRlab project[5,6] (https://qmrlab.org/).

Table 1 summarizes the parameters that can be arbitrarily adjusted for the simulation. We describe how the simulation is performed.

Table 1

Arbitrarily adjustable parameters for Monte Carlo simulation

Name	Brief instruction	Remarks
Mode	2D or 3D. The field to be simulated	Select one
Echo time [ms]	The total time span of the simulation
Time-step length [ms]	The unit time span for simulating molecular displacement and MPG
Obstacles
Number	The number of obstacles
Position and size [μm]	The center and radius of each circle (2D mode) or the center, radius, and height of each cylinder (3D mode)	*
Penetration	The probability that a molecule will penetrate through obstacles without being reflected	Zero to one*
Molecules
Number	The number of molecules
Positions [μm]	The initial positions of the molecules. Usually defined at random so that they are distributed within a certain target area	*
Directions	The initial directions of the molecules to be displaced. Usually defined at random	*
Flight [μm]	Mean and standard deviation of the distance that the molecules moves within each time-step
MPG
Strength [mT/m]	Strength of the gradient in each time-step
Direction	Matrix that defines the direction of the gradient (vector) in each time-step

MPG: motion probing gradient. Remarks: limitation for parameters are indicated, blank means that there is no limitation.

the parameter can be defined for each obstacle or molecule independently, but for convenience, assigning random number according to appropriate rules (e.g. define only mean and standard deviation) is supposed.

Dimensionality

Either 2D mode or 3D mode can be selected. 2D mode simulates a plane field that has arbitrarily arranged circles of arbitrary sizes, where the circles work as obstacles that represent cell membranes. The circles can be overlapped if needed, but in general use, an axon-packing algorithm[4] automatically and randomly packs the circles so that their radii have an arbitral mean and variance, with arbitral shortest distance between circles. In 3D mode, each circle becomes a cylinder that can be assumed to represent parallel neural fibers.

Brownian Motion

The proposed software simulates the water diffusion using the simple model. First, the water molecules are assumed as random walkers that do not interact with each other. Second, the water molecules interact with cell membranes (i.e. circles and cylinders defined in the previous section) by either total internal reflection or transmission, which occurs randomly according to the arbitrary defined proportion. The speed of random walking at intra- and extracellular spaces are the same.

The entire simulation is performed by repeating a short simulation corresponding to an arbitrarily defined unit of time (a time-step). Fig. 1a is the flowchart that explains how each water molecule movements are calculated within each time-step. In each time-step, the water molecules start diffusing from their existing positions in directions determined randomly for each molecule at the beginning of each time-step. The molecules move linearly, but when they encounter the borders of circles or cylinders, they either are reflected or they penetrate the borders with an arbitrarily determined probability (Fig. 1b). The molecules continue to move and to repeat these reflections and penetrations until the given distances for the time-step (that is determined automatically according to the mean and variance given be the user) expire. The time-step is repeated until the total time reaches the echo time.

Fig. 1

Outline of the simulation algorithm.

Simulation of the Motion Probing Gradient

The direction and magnetic gradient strength that form the motion probing gradient (MPG) can be defined at each time-step. In each time-step, the relative phase of each molecule shifts depending on the strength of the additional magnetization (i.e., MPG) it has experienced during that time-step. This phase shift is calculated separately for each linear movement within the time-step (i.e., for movements separated by reflections), so that the effect of the MPG can be accurately simulated even when the molecule has been reflected many times during the time-step.

Outputs of the Simulation

The positions of the water molecules in each time-step are visualized as animated dots moving in the fields of circles or cylinders. In addition, graphs of mean squared displacement, strength and direction of the MPG, and phase dispersion level (S) reflecting the signal intensity at each time-step are plotted. Here, S is calculated as, where n is the number of molecules and θ is the relative phase shift of molecular k in radian.

Below are two specific examples of the outputs (Full data are available as Movies 1–4 and Supplementary Figs. 1–4 in the online Supplementary Materials.).

Output example 1: Restricted/hindered diffusion between intra- and extracellular spaces. The simulation was performed in 2D mode. One thousand molecules were distributed randomly within a 20-m2 as the initial state, with red molecules inside the cells and blue ones outside the cells (Fig. 2a; also see Movie 1 and Supplementary Fig. 1. Some of the major parameters are listed in Supplementary Table 1). The blue (extracellular) dots were then removed from the initial state and the simulation was performed (Fig. 2b). The permeability of the cell membranes was set to zero, giving pure restricted diffusion. Therefore, the mean square displacement of the molecules reached a plateau in this setting (Fig. 2, graph b-2) The signal intensity (S/S0) and the apparent diffusion coefficient (ADC) were obtained by simulating a pair of MPGs to give a b-value = 1000 [s/mm2] and Δ/δ = 28/25 [ms] (not shown in the figure). The simulation was then repeated, but the red dots were removed instead of the blue ones (Fig. 2c). This configuration is called hindered diffusion, and is characterized by approximately linearly increasing mean square displacement over time (Fig. 2, graph c-2.) Therefore, the ADC is larger than with restricted diffusion, and this difference increases as the diffusion time increases. We also performed a simulation similar to those shown in Figs. 1b and 1c but using both red and blue dots (Fig. 2d). The mean square displacement graph has characteristics similar to those of Fig. 1b at the start, but it does not reach a plateau (Fig. 2, graph d-2.) We then performed a simulation similar to that in Fig. 1d, but when the dots encountered a cell membrane, they penetrated it instead of being reflected, with an average 3 ± 0.3% probability of penetration (Fig. 2e). The dots diffused more than in Fig. 1d, and this is reflected in the S/S0 and ADC values.

Fig. 2

Output example 1: Restricted/hindered diffusion between intra- and extracellular spaces.

Output example 2: Measuring water exchange by using DDE. These simulations (Fig. 3; also see Movies 2–4 and Supplementary Figs. 2–4. Some of the major parameters are listed in Supplementary Table 1) confirmed a basic idea introduced in an earlier report, namely filter exchange imaging.[3] The simulated field was similar to those in Figs. 1d and 1e, but DDE was performed to visualize molecule exchange between the intra- and extracellular spaces. In the first simulation (Fig. 3a) DDE was performed by using two different time intervals between the MPG pairs, 5 ms and 25 ms. The two MPG pairs were the same; each was set to have a b-value = 800 [s/mm2] and Δ/δ = 14/13 [ms]. The first MPG pair is used to filter the signals of molecules with high diffusion coefficient that exist mainly in the extracellular space. Then, the second MPG pair is used to observe how molecules in the intra- and extracellular spaces are exchanged during the mixing time. When there is no exchange, the different mixing times do not alter the final signal intensity, because the conditions of the field do not change during the mixing time (Fig. 3b) (The effect of T2 relaxation is not simulated, so this signal intensity corresponds to the signal intensity normalized by that without MPGs.). When there was 3% exchange, the molecules that existed mainly inside the cells at the first MPG pair, which relatively maintain phase coincidence and can therefore provide higher signal intensity than those outside the cells, permeated to the extracellular space during the mixing time, where they could move more freely (Fig. 3c). The observed signal intensity after the second MPG pair may become lower if more of such molecules exist in the extracellular space at the time of the second MPG. Therefore, the signal will decrease with the longer mixing time, as simulated in the graph.

Fig. 3

Output example 2: Measuring water exchange by using DDE.

Comparison with the Previously Published Monte Carlo Simulator

Previously, Yeh et al.[7] published a Monte Carlo simulation software which the basic models for simulating the Brownian motion and interactions with cell membrane were similar with our proposed one. In addition, their software can simulate more complicated diffusing field than ours; however, the proposed software has several merits instead. First, MPGs can be designed more freely. Bipolar MPG for example, and even other unknown MPG patterns based on new ideas are available, because the intensity and direction of the gradient can be defined individually at each time step. Second, the dynamics of dephasing and rephasing are visualized which may help the user understanding the simulated phenomenon comprehensively. Third, in combination with the 2D axon packing software,[4] the users can easily and quickly simulate multiple different diffusing fields that have similar cell radii (mean and variance) and distance between the cells, that might be useful to confirm the reproducibility of the results.

Limitations of the Simulation

So far, the obstacles in the simulation field are only circles or cylinders arranged in parallel. This point would be resolved in the future version.

Supplementary files below are available online.

Supplementary Table 1

Values of parameters used in Monte Carlo simulations.

Supplementary Fig. 1

Restricted and hindered diffusion.

Supplementary Fig. 2

Water exchange or no water exchange between the intra- and extracellular spaces.

Supplementary Fig. 3

Double diffusion encoding (DDE) with water exchange.

Supplementary Fig. 4

Double diffusion encoding (DDE) without water exchange.