Tutorial

This tutorial will walk through an example of modelling the MRI signal evolution for a diffusion-weighted sequence. The spins in this simulation will be constrained by regularly packed cylinders. This tutorial will use the programatic Julia interface, which you can run in the Julia REPL or in a Jupyter notebook. If you would prefer to use the command line interface, you can find a tutorial doing the same analysis here.

After installation we can load MCMRSimulator.jl using

using MCMRSimulator
using CairoMakie  # used for plotting; use GLMakie or WGLMakie for interactive plots
update_theme!(Theme(Axis=(xgridvisible=false, ygridvisible=false))) # hide grid lines

In general, running a simulation will consist of the following steps:

• Defining the microstructure and one or more sequences by creating an appropriate Simulation object.
• Simulating a random walk of the spins through the microstructure and the MR signal produced by those spins. For signal simulations, the recommended approach is adaptive readout using a target SNR.
• Plotting the MR signal or storing it to disk.

We will look through each of these steps below.

Defining the simulation

The first step is to define the environment through which the spins will evolve. We will do so by creating an appropriate Simulation object. This Simulation will contain information on the microstructure, how spins interact with that microstructure, and the enabled sequence(s).

These different steps are described in more detail in other sections of this documentation:

• How to define the microstrutural geometry
• Controlling spin behaviour

For a discussion on sequence generation see MRIBuilder.jl.

First we will define a geometry formed of regularly packed axons. This is represented by a single cylinder with a radius of 1 micrometer that repeats itself every 2.5 micrometer (in both the x-, and y-direction).

geometry = Cylinders(radius=1., repeats=[2.5, 2.5])

f = plot(PlotPlane(size=5), geometry)
xlims!(f.axis, -2.5, 2.5)
ylims!(f.axis, -2.5, 2.5)
f

More complicated geometries can be generated as described here. More details on plotting geometries can be found in the plot_geometry documentation.

The next step is to load a sequence. Here we load a pre-defined diffusion-weighted MRI sequence.

The next step is to define a sequence using MRIBuilder.jl. Here we will adopt a diffusion-weighted MRI sequence defined in pulseq.

sequence = read_pulseq("dwi_te_80_bval_2.seq", TR=300)
f = plot(sequence)
f

Once we have both a geometry and one or more sequences, we can put them together in a Simulation object:

simulation = Simulation(sequence, R2=0.012, R1=3e-3, diffusivity=2., off_resonance=0.1, geometry=geometry)

┌ Info: # Timestep determination
│ Maximum timestep set by tortuosity constraint based on size of geometry to 0.015 ms.
│ Size scale of smallest object in the simulation was automatically determined to be 1.0 um.
│ If this value is too small, you can set the size scale explicitly by passing on `size_scale=<new_value>` to the `Simulator` constructor.
└ The actual timestep will be reduced based on the MR sequence(s).

By default there is no T1 or T2 relaxation and a diffusivity of 3 um^2/ms. Enabling spin relaxation and diffusion requires setting the appropriate parameters in the Simulation as seen here. The spin relaxation rates can be different in different compartments as described here.

Simulating a cell from an SWC file

The SWC specification defines a standard text format for digital reconstructions. SWC files can be loaded directly as connected sphere-and-cylinder geometries. For example, the tutorial includes a small branched cell morphology:

cell_geometry = read_swc("cell.swc")
cell_simulation = Simulation(sequence, diffusivity=2., geometry=cell_geometry)

┌ Info: # Timestep determination
│ Maximum timestep set by tortuosity constraint based on size of geometry to 0.96 ms.
│ Size scale of smallest object in the simulation was automatically determined to be 8.0 um.
│ If this value is too small, you can set the size scale explicitly by passing on `size_scale=<new_value>` to the `Simulator` constructor.
└ The actual timestep will be reduced based on the MR sequence(s).

The default sampling box for this finite, non-repeating geometry is calculated from the cell itself. To simulate only initially intracellular spins, use a Subset as the adaptive readout filter:

intracellular_signal = readout(
    cell_simulation;
    target_snr=1,
    max_spins=100,
    filter=Subset(inside=true),
)
transverse(intracellular_signal; mean=true)

0.10831796614958912

The filter is applied before the random walks start, so spins outside the cell are not simulated. This differs from subset, which selects spins from the signal after they have been simulated.

Initialising the simulation

The current state of the simulation at any time is given by a Snapshot object. This is essentially a vector of Spin objects with a time stamp. Each Spin represents a single diffusing particle. Besides containing its current position, it also contains its contribution to the MR signal for each of the sequences in the simulation and whether it is stuck on any surfaces.

When a full spin state is needed, the recommended way to initialise is to call Snapshot(<number of spins>, <simulation>, [bounding_box]). This will create randomly distributed spins within some BoundingBox. By default this bounding box is the finite geometry bounds for non-repeating geometries, or an isotropic voxel with a size of 1 mm centered on the origin otherwise.

After initialisation or after running the simulation, the Snapshot can be later filtered to include only spins inside/outside specific compartments or only include free/bound spins using [get_subset].

The simulation can also be initialised explicitly using a sequence of positions (i.e., length-3 vectors) with the initial spin positions. Note that such a simulation will start with all spins free and not necessarily randomly distributed, which means it might take some time to reach an equilibrium.

For each of these initialisations the initial magnetisation can be explicitly set using the transverse, longitudinal, and phase flags. The default is for spins to start in equilibrium (i.e., transverse magnetisation of 0 and longitudinal magnetisation of 1).

Finally, one could start a simulation using a Snapshot from a previous simulation.

Running the simulation

The main way to run a signal simulation is by calling readout with a target SNR. The simulator then generates independent batches of spins until the requested precision is reached. Fixed-spin readouts are also available when an explicit number of spins is required. There are various ways to define when the output will be read out (as described in the readout documentation).

Here, we will illustrate various examples of using this function:

Adaptive signal readouts

Most sequences will contain one or more ADC events, which define when the sequence will be read out during each repetition time (TR). To get the signal at this time, we can simply call:

signal = readout(simulation; target_snr=1, max_spins=100)

SpinOrientationSum(longitudinal=1.2786987632228022, transverse=22.86136024598594, phase=94.93985651228182°, nspins=100, snr=[12.118069061531715, 9.356141801032033, 3.295505035590032e15])

The returned MCMRSimulator.SpinOrientationSum stores the total signal and the number of contributing spins. For values that are easier to compare between simulations with different number of spins, request the per-spin mean:

transverse(signal; mean=true), longitudinal(signal; mean=true)

(0.2286136024598594, 0.012786987632228021)

This signal is not truely representative from what we expect in a true diffusion-weighted MRI sequence, because the longitudinal signal has not had a chance to relax across multiple repetition times. To see what the signal will look like after such equilibriation, we can delay our readout with several TRs:

readout(simulation; target_snr=1, max_spins=100, skip_TR=2)

SpinOrientationSum(longitudinal=5.075264340325728, transverse=2.5310157340205492, phase=-153.6677898024463°, nspins=100, snr=[0.7844515538147622, 0.8696824821480094, 39.10262064432811])

In addition, to the total signal, we can also get the signal associated with individual compartments:

signals_by_subset = readout(simulation; target_snr=1, max_spins=100, subset=[Subset(inside=true), Subset(inside=false)])

2-element Vector{MCMRSimulator.Spins.SpinOrientationSum}:
 SpinOrientationSum(longitudinal=0.6649233568758605, transverse=19.90650341003997, phase=89.9653845853568°, nspins=52, snr=[359.55167670745493, 32340.412822572172, 2.339414434409847e15])
 SpinOrientationSum(longitudinal=0.6137754063469311, transverse=1.2027993136445798, phase=-155.06026527019097°, nspins=48, snr=[0.6138589395653459, 0.6604398008769443, 2.282174993884324e15])

Note that we now get two signal outputs. The first respresents the signal within the cylinders, which is very close to number of spins, indicating that there has been very little dephasing due to the diffusion weighting inside the cylinders. On the other hand, we did lose most of the signal outside of the cylinders. All the spins are either inside or outside the cylinders, so in this case the first row is simply the sum of the next two.

Instead of just running the simulation for multiple TRs without readouts, we could also visualise the equilibriation process by outputting the signal for multiple TRs:

signals = readout(simulation; target_snr=1, max_spins=100, nTR=6)
f = lines(longitudinal.(signals))
lines!(transverse.(signals))
f

┌ Warning: Adaptive readout did not reach target SNR 1 after 100 spins.
└ @ MCMRSimulator.Evolve /builds/ndcn0236/mcmrsimulator.jl/src/evolve.jl:475

At each timepoint readout returns the total MR signal (for each sequence) as a MCMRSimulator.SpinOrientationSum object. The snr field contains the independently estimated SNR for the Sx, Sy, and Sz components when using adaptive readout. From this one can estimate the transverse component, the longitudinal component, and the phase. The longitudinal and transverse functions are used above to get those respective components.

We can also override, when the signal will be read out, by passing on the readout times as a third argument to readout. Here we use this to plot the actual transverse signal evolution.

times = 0:0.1:100
# simulate 3000 spins for a single repetition time
average_signals = readout(simulation; target_snr=1, max_spins=100, readout_times=times)
f = lines(times, transverse.(average_signals; mean=true))

┌ Warning: Adaptive readout did not reach target SNR 1 after 100 spins.
└ @ MCMRSimulator.Evolve /builds/ndcn0236/mcmrsimulator.jl/src/evolve.jl:475

Fixed-spin simulations

When an exact number of spins is needed, pass that number explicitly:

fixed_signal = readout(1000, simulation)
transverse(fixed_signal; mean=true)

0.21131332278689735

This is useful for controlled comparisons with a more predictable runtime or benchmarking. Adaptive readout is preferable for ordinary signal simulations because it controls the Monte Carlo noise directly.

Reading out the full snapshot

Instead of returning just the total signal readout can also return the full spin state in a Snapshot object by setting the return_snapshot keyword to true. Note that this is very memory intensive, so is only recommended when you only output a small number of timepoints or a small number of spins.

Here, we use this to visualise the trajectory of spins through the geometry. To plot the trajectory we first need to output the state of the all spins at a high temporal resolution:

# Simulate 2 spins with given starting positions for 3 ms
snapshots = readout([[0, 0, 0], [1, 1, 0]], simulation, 0:0.01:3., return_snapshot=true)

pp = PlotPlane(size=5.)
f = plot(pp, geometry)
plot(pp, snapshots)
f

In this plot the color at each timepoint encodes the spin orientation. The brightness of the spin indicates the size of the transverse component with purely longitudinal spins being in black. The color of the spin encodes the phase of the MR signal in the transverse plane.

The trajectories can also be plotted in 3D:

f = plot(snapshots)

More details on plotting trajectories can be found in the plot_trajectory documentation.

We can also use this future to plot the complete Snapshot at a specific time. In this example we do not set this time explicitly, so it will default to the time of the sequence MRIBuilder.SingleReadout/MRIBuilder.ADC as discussed above:

snapshot = readout(3000, simulation, return_snapshot=true)
pp = PlotPlane(size=2.5)
f = plot(pp, snapshot)
plot!(pp, geometry)
f

We can see that the spins within the cylinder stay far more in-sync with each other (i.e., they have a more uniform colour), than those outside of the cylinder.

Note that object returned when return_snapshot=true is the same Snapshot object as was used to initialise this simulation. This means that it can be used as an initialisation for future simulations.

One complication with this occurs when running a Simulation with multiple sequences. Often, these different sequences will readout at different times. So, each Snapshot object returned by readout will only contain the spin magnetisation of the sequence that is readout at that particular time. To get a Snapshot object with the spin states for all sequences, you can use the evolve function instead.

More details on plotting snapshots can be found in the plot_snapshot documentation. For example, the same snapshot can also be plotted using dyads or as an image:

f = plot(pp, snapshot, kind=:dyad)
plot!(pp, geometry)
f

f = plot(pp, snapshot, kind=:image)
plot!(pp, geometry)
f