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 command line interface, which we assume is available through the mcmr command (see installation instructions). If you would prefer to interact with MCMRSimulator in Julia, you can find a tutorial doing the same analysis here.
The command line tool requires one or more MRI sequence definitions as input. These can be supplied as pulseq files, such as those produced by MRIBuilder.jl.
In general, running a simulation will consist of the following two steps:
- Creating a geometry using one or more calls to
mcmr geometry(full description). - Running the actual simulation using
mcmr run.
We will look through each of these steps below.
Defining the geometry
First we will define a geometry formed of regularly packed axons. This is represented by a single cylinder pointing in the z-direction with a radius of 1 micrometer that repeats itself every 2.5 micrometer (in both the x-, and y-direction).
mcmr geometry create cylinders 1 geometry.json --radius 1 --repeats 2.5,2.5This will create a JSON file with the full information on the geometry:
{
"type": "Cylinders",
"number": 1,
"#radius_description": "Radius of the cylinder. Field is required. Expected type: Float64.",
"radius": 1.0,
"#rotation_description": "Rotation applied to all obstructions in group. Can be set to a matrix or one of :x, :y, or, :z (see [`MCMRSimulator.get_rotation`](@ref MCMRSimulator.Methods.get_rotation)). Field is required. Expected type: StaticArraysCore.SMatrix{3, 2, Float64, 6}.",
"rotation": [[1.0,0.0,0.0],[0.0,1.0,0.0]],
"#R1_surface_description": "Additional longitudinal relaxation rate (kHz). Surface property. Field is required. Expected type: Float64.",
"R1_surface": 0.0,
"#R1_inside_description": "Additional longitudinal relaxation rate (kHz). Inside property. Field is required. Expected type: Float64.",
"R1_inside": 0.0,
"#R2_surface_description": "Additional transverse relaxation rate (kHz). Surface property. Field is required. Expected type: Float64.",
"R2_surface": 0.0,
"#R2_inside_description": "Additional transverse relaxation rate (kHz). Inside property. Field is required. Expected type: Float64.",
"R2_inside": 0.0,
"#off_resonance_surface_description": "Additional off-resonance field offset (kHz). Surface property. Field is required. Expected type: Float64.",
"off_resonance_surface": 0.0,
"#off_resonance_inside_description": "Additional off-resonance field offset (kHz). Inside property. Field is required. Expected type: Float64.",
"off_resonance_inside": 0.0,
"#position_description": "Spatial offset of obstruction from origin. Field is required. Expected type: StaticArraysCore.MVector{2, Float64}.",
"position": [0.0,0.0],
"#g_ratio_description": "Inner/outer radius used for susceptibility calculation Field can be null. Expected type: Float64.",
"g_ratio": 1.0,
"#susceptibility_iso_description": "Isotropic component of the susceptibility (in ppm um). It can be calculated by multiplying the isotropic susceptibility of the simulated tissue by the thickness. Field can be null. Expected type: Float64.",
"susceptibility_iso": -0.1,
"#susceptibility_aniso_description": "Anisotropic component of the susceptibility (in ppm um). It can be calculated by multiplying the anisotropic susceptibility of the simulated tissue by the thickness Field can be null. Expected type: Float64.",
"susceptibility_aniso": -0.1,
"#lorentz_radius_description": "Only compute field explicitly for cylinders with this Lorentz radius. Field can be null. Expected type: Float64.",
"lorentz_radius": 5.0,
"#repeats_description": "Length scale on which the obstructions are repeated (um). Field can be null. Expected type: StaticArraysCore.MVector{2, Float64}.",
"repeats": [2.5,2.5],
"#use_boundingbox_description": "Use bounding boxes for an initial filtering of possible intersections. Field can be null. Expected type: Bool.",
"use_boundingbox": true,
"#grid_resolution_description": "Resolution of the grid that the volume is split up into (um). Defaults to roughly one grid element per obstruction. Field can be null. Expected type: Float64.",
"grid_resolution": null,
"#dwell_time_description": "Average time a particle stays stuck to the surface (ms). Surface property. Field can be null. Expected type: Float64.",
"dwell_time": null,
"#density_description": "Surface density of stuck particles relative to the volume density (um). Surface property. Field can be null. Expected type: Float64.",
"density": null,
"#permeability_description": "Rate of particle passing through the obstruction in arbitrary units. Surface property. Field can be null. Expected type: Float64.",
"permeability": null,
"#relaxation_description": "Rate of signal loss at each collision. The actual signal loss at each collision is e^(-x * sqrt(t)), where x is this rate and t is the timestep. Surface property. Field can be null. Expected type: Float64.",
"relaxation": null,
"#size_scale_description": "Size of the smallest obstructions. If not set explicitly, this will be determined by the minimum radius or distance between objects (see `size_scale`). Field can be null. Expected type: Float64.",
"size_scale": null
}
You can see how that the repeats and radius keywords have been set to our predefined values. You can alter these and other geometry properties by editing this JSON directly or using the flags when creating the geometry. For a full overview of these flags, you can run:
mcmr geometry create cylinders --helpusage: mcmr geometry create cylinders --radius RADIUS
[--rotation ROTATION]
[--R1_surface R1_SURFACE]
[--R1_inside R1_INSIDE]
[--R2_surface R2_SURFACE]
[--R2_inside R2_INSIDE]
[--off_resonance_surface OFF_RESONANCE_SURFACE]
[--off_resonance_inside OFF_RESONANCE_INSIDE]
[--position POSITION] [--g_ratio G_RATIO]
[--susceptibility_iso SUSCEPTIBILITY_ISO]
[--susceptibility_aniso SUSCEPTIBILITY_ANISO]
[--lorentz_radius LORENTZ_RADIUS]
[--repeats REPEATS] [--no-use_boundingbox]
[--grid_resolution GRID_RESOLUTION]
[--dwell_time DWELL_TIME] [--density DENSITY]
[--permeability PERMEABILITY]
[--relaxation RELAXATION]
[--size_scale SIZE_SCALE] [-h] number
output_file
Create a geometry JSON file filled with only cylinders with any
properties defined by the flags.
positional arguments:
number Number of obstructions to create. (type:
Int64)
output_file Geometry JSON output filename.
optional arguments:
--radius RADIUS Radius of the cylinder. Field is required.
Expected type: Float64. (type:
FieldParser{Float64})
--rotation ROTATION Rotation applied to all obstructions in group.
Can be set to a matrix or one of :x, :y, or,
:z (see [`MCMRSimulator.get_rotation`](@ref
MCMRSimulator.Methods.get_rotation)). Field is
required. Expected type:
StaticArraysCore.SMatrix{3, 2, Float64, 6}.
default value: [1.0 0.0; 0.0 1.0; 0.0 0.0]
(type: RotationParser{2}, default:
RotationParser{2}([1.0 0.0; 0.0 1.0; 0.0
0.0]))
--R1_surface R1_SURFACE
Additional longitudinal relaxation rate (kHz).
Surface property. Field is required. Expected
type: Float64. default value: 0.0 (type:
FieldParser{Float64}, default:
FieldParser{Float64}(0.0))
--R1_inside R1_INSIDE
Additional longitudinal relaxation rate (kHz).
Inside property. Field is required. Expected
type: Float64. default value: 0.0 (type:
FieldParser{Float64}, default:
FieldParser{Float64}(0.0))
--R2_surface R2_SURFACE
Additional transverse relaxation rate (kHz).
Surface property. Field is required. Expected
type: Float64. default value: 0.0 (type:
FieldParser{Float64}, default:
FieldParser{Float64}(0.0))
--R2_inside R2_INSIDE
Additional transverse relaxation rate (kHz).
Inside property. Field is required. Expected
type: Float64. default value: 0.0 (type:
FieldParser{Float64}, default:
FieldParser{Float64}(0.0))
--off_resonance_surface OFF_RESONANCE_SURFACE
Additional off-resonance field offset (kHz).
Surface property. Field is required. Expected
type: Float64. default value: 0.0 (type:
FieldParser{Float64}, default:
FieldParser{Float64}(0.0))
--off_resonance_inside OFF_RESONANCE_INSIDE
Additional off-resonance field offset (kHz).
Inside property. Field is required. Expected
type: Float64. default value: 0.0 (type:
FieldParser{Float64}, default:
FieldParser{Float64}(0.0))
--position POSITION Spatial offset of obstruction from origin.
Field is required. Expected type:
StaticArraysCore.MVector{2, Float64}. default
value: [0.0, 0.0] (type:
FieldParser{MVector{2, Float64}}, default:
FieldParser{MVector{2, Float64}}([0.0, 0.0]))
--g_ratio G_RATIO Inner/outer radius used for susceptibility
calculation Field can be null. Expected type:
Float64. default value: 1.0 (type:
FieldParser{Float64}, default:
FieldParser{Float64}(1.0))
--susceptibility_iso SUSCEPTIBILITY_ISO
Isotropic component of the susceptibility (in
ppm um). It can be calculated by multiplying
the isotropic susceptibility of the simulated
tissue by the thickness. Field can be null.
Expected type: Float64. default value: -0.1
(type: FieldParser{Float64}, default:
FieldParser{Float64}(-0.1))
--susceptibility_aniso SUSCEPTIBILITY_ANISO
Anisotropic component of the susceptibility
(in ppm um). It can be calculated by
multiplying the anisotropic susceptibility of
the simulated tissue by the thickness Field
can be null. Expected type: Float64. default
value: -0.1 (type: FieldParser{Float64},
default: FieldParser{Float64}(-0.1))
--lorentz_radius LORENTZ_RADIUS
Only compute field explicitly for cylinders
with this Lorentz radius. Field can be null.
Expected type: Float64. default value: 5.0
(type: FieldParser{Float64}, default:
FieldParser{Float64}(5.0))
--repeats REPEATS Length scale on which the obstructions are
repeated (um). Field can be null. Expected
type: StaticArraysCore.MVector{2, Float64}.
(type: FieldParser{MVector{2, Float64}})
--no-use_boundingbox Use bounding boxes for an initial filtering of
possible intersections. Field can be null.
Expected type: Bool. default value: true
--grid_resolution GRID_RESOLUTION
Resolution of the grid that the volume is
split up into (um). Defaults to roughly one
grid element per obstruction. Field can be
null. Expected type: Float64. (type:
FieldParser{Float64})
--dwell_time DWELL_TIME
Average time a particle stays stuck to the
surface (ms). Surface property. Field can be
null. Expected type: Float64. (type:
FieldParser{Float64})
--density DENSITY Surface density of stuck particles relative to
the volume density (um). Surface property.
Field can be null. Expected type: Float64.
(type: FieldParser{Float64})
--permeability PERMEABILITY
Rate of particle passing through the
obstruction in arbitrary units. Surface
property. Field can be null. Expected type:
Float64. (type: FieldParser{Float64})
--relaxation RELAXATION
Rate of signal loss at each collision. The
actual signal loss at each collision is e^(-x
* sqrt(t)), where x is this rate and t is the
timestep. Surface property. Field can be null.
Expected type: Float64. (type:
FieldParser{Float64})
--size_scale SIZE_SCALE
Size of the smallest obstructions. If not set
explicitly, this will be determined by the
minimum radius or distance between objects
(see `size_scale`). Field can be null.
Expected type: Float64. (type:
FieldParser{Float64})
-h, --help show this help message and exit
How these various properties affect the simulation is described here.
The procedure to create Walls, Spheres, or Annuli is very similar as for the Cylinders illustrated above. Randomly distributed cylinders, annuli, and spheres can be created using mcmr geometry create-random.
Simulating a cell from an SWC file
The SWC specification defines a standard text format for digital reconstructions. An SWC morphology can be passed directly to mcmr run; no intermediate JSON geometry file is needed. The tutorial includes a small branched cell morphology in cell.swc. SWC nodes are loaded as connected spheres and cylinders by default, and the adaptive --filter-inside option restricts the simulation to spins that start inside the cell:
mcmr run cell.swc dwi_te_80_bval_2.seq --target-snr 10 --filter-inside -o intracellular_signal.csvdf = 1×18 DataFrame
Row │ sequence sequence_index TR readout subset nspins longitudinal transverse phase Sx Sy Sz SE_Sx SE_Sy SE_Sz SNR_Sx SNR_Sy SNR_Sz
│ SubStrin… Int64 Int64 Int64 Int64 Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ dwi_te_80_bval_2.seq 1 1 1 0 33 -2.01859e-16 0.136862 -49.3797 0.089103 -0.103884 -2.01859e-16 4.17069 3.99989 5.84846e-16 1.0829 1.12914 5.64251e16
┌ Info: # Timestep determination
│ Maximum timestep set by tortuosity constraint based on size of geometry to 0.64 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).
For finite, non-repeating geometries, adaptive spin initialisation uses the geometry bounds by default. The --filter-inside option is only available for adaptive simulations, because it filters each generated batch before it is simulated.
Running the simulation
To get instructions on running the simulations, we can check the help message of mcmr run:
mcmr run --helpusage: mcmr run [-D DIFFUSIVITY] [--R1 R1] [--R2 R2]
[--swc-as-spheres] [-o OUTPUT-SIGNAL]
[--output-snapshot OUTPUT-SNAPSHOT] [--nTR NTR]
[-T TIMES [TIMES...]] [--skip-TR SKIP-TR]
[--target-snr TARGET-SNR] [--max-spins MAX-SPINS]
[--batch-size BATCH-SIZE] [--filter-inside]
[--subset SUBSET [SUBSET...]] [-N NSPINS]
[--voxel-size VOXEL-SIZE]
[--longitudinal LONGITUDINAL]
[--transverse TRANSVERSE] [--phase PHASE]
[--seed SEED] [-h] geometry [sequence...]
Runs a Monte Carlo simulation of the MRI signal evolution for spins
interacting with the geometry.
optional arguments:
-h, --help show this help message and exit
Define the simulation parameters:
geometry Geometry input file. JSON files describe
obstructions and their biophysical properties;
PLY files provide meshes; SWC files are loaded
as connected spheres and cylinders by default.
Use --swc-as-spheres to load overlapping
spheres without connecting cylinders.
sequence One of more pulseq .seq files describing the
sequences to be run.
-D, --diffusivity DIFFUSIVITY
Diffusivity of free water (um^2/ms). (type:
Float64, default: 3.0)
--R1 R1 Longitudinal relaxation in 1/ms. This
relaxation rate will at the very least be
applied to free, extra-cellular spins. It
might be overriden in the 'geometry' for bound
spins or spins inside any obstructions. (type:
Float64, default: 0.0)
--R2 R2 Transverse relaxation in 1/ms. This relaxation
rate will at the very least be applied to
free, extra-cellular spins. It might be
overriden in the 'geometry' for bound spins or
spins inside any obstructions. (type: Float64,
default: 0.0)
--swc-as-spheres If set, SWC files will be read as overlapping
spheres instead of as a graph of spheres and
cylinders.
Output flags. At least one is required:
-o, --output-signal OUTPUT-SIGNAL
Writes the total signal at the readouts to
this file as a comma-separated value (CSV)
table.
--output-snapshot OUTPUT-SNAPSHOT
Writes the state of all the spins at the
readouts to this file as a comma-separated
value (CSV) table.
Readout flags. These control when the signal/spin states will be read out:
--nTR NTR Acquire the signal provided at the sequence
readouts for this many repetition times (TRs).
Output will be stored as a CSV file. (type:
Int64, default: 1)
-T, --times TIMES [TIMES...]
Acquire the signal at the given times within
each TR (in ms). Multiple values can be
provided (e.g., '-T 0 10 15.3'). By default,
the Readout markers in the sequence will be
used instead. (type: Float64)
--skip-TR SKIP-TR The number of repetition times the simulation
will run before starting to acquire data.
(type: Int64, default: 0)
--target-snr TARGET-SNR
Adaptively simulate spins until the target
signal-to-noise ratio is reached for every
subset and magnetisation component. (type:
Float64)
--max-spins MAX-SPINS
Maximum number of spins for adaptive
simulations. (type: Int64)
--batch-size BATCH-SIZE
Number of spins added per batch in adaptive
simulations. (type: Int64, default: 1000)
--filter-inside In adaptive simulations, only simulate spins
that are initially inside the geometry.
--subset SUBSET [SUBSET...]
Can be provided multiple times. For each time
it is provided, the signal will be computed at
each readout for a specific subset of spins.
This subset is defined by one or two values
from bound/free/inside/outside. Afterwards
they can include an integer value to select a
specific geometry, and an additional integer
value to select a specific obstruction within
that geometry. The selected user geometry is
passed to the subset. For example: - `--subset
free`: include any free spins - `--subset
inside`: include any spins inside any geometry
- `--subset outside 2`: include any spins
outside of the second obstruction group in the
geometry - `--subset inside bound 2 3`:
include any spins inside and bound to the 3rd
obstruction in the second obstruction group of
the geometry.
Initialisation flags. These control the spins initial state:
-N, --Nspins NSPINS Number of spins to simulate. Ignored if --init
is set. (type: Int64)
--voxel-size VOXEL-SIZE
Size of the voxel (in mm) over which the
initial spins are spread. (type: Float64,
default: 1.0)
--longitudinal LONGITUDINAL
Initial value of the longitudinal
magnetisation for each spin. Note the the
equilibrium longitudinal magnetisation for
each spin is 1. (type: Float64, default: 1.0)
--transverse TRANSVERSE
Initial value of the magnitude of the
transverse magnetisation for each spin. (type:
Float64, default: 0.0)
--phase PHASE Initial value of the phase of the transverse
magnetisation for each spin in degrees. (type:
Float64, default: 0.0)
--seed SEED Initialisation for random number seed. Supply
this to get reproducible results. If --init is
also set, this flag will override the seed
stored in this initialisation file. (type:
Int64)
We can see that in addition to defining the geometry and the sequence, we can also control the simulation properties such as the --diffusivity, --R1, and --R2.
The recommended approach is to set a target SNR and let the simulator determine how many spins are needed. For adaptive simulations, spins are randomly distributed within the geometry bounds for finite, non-repeating geometries, or within the default 1 mm voxel otherwise. Fixed-spin simulations use the voxel size given by --voxel-size. The initial state might also contain bound spins (if the --density flag was set to a non-zero value during the geometry generation).
We will use a pre-defined diffusion-weighted MRI sequence with an echo time of 80 ms and a b-value of 2 mm^2/s.
The DWI sequence defined above contains a single ADC event at the echo time (80 ms). By default, this is used for readout:
mcmr run geometry.json dwi_te_80_bval_2.seq --target-snr 10 -o signal.csvdf = 1×18 DataFrame
Row │ sequence sequence_index TR readout subset nspins longitudinal transverse phase Sx Sy Sz SE_Sx SE_Sy SE_Sz SNR_Sx SNR_Sy SNR_Sz
│ SubStrin… Int64 Int64 Int64 Int64 Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ dwi_te_80_bval_2.seq 1 1 1 0 100 -2.10942e-16 0.474785 82.1518 0.0648312 0.470338 -2.10942e-16 5.18407 7.16699 1.31553e-15 9.15855 6.62461 7.60151e16
┌ Info: # Timestep determination
│ Maximum timestep set by tortuosity constraint based on size of geometry to 0.01 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).
This produces the CSV file, which looks like
sequence,sequence_index,TR,readout,subset,nspins,longitudinal,transverse,phase,Sx,Sy,Sz,SE_Sx,SE_Sy,SE_Sz,SNR_Sx,SNR_Sy,SNR_Sz
dwi_te_80_bval_2.seq,1,1,1,0,100,-2.1094237467877973e-16,0.47478521084356556,82.1518236545574,0.0648312466602298,0.4703380761667606,-2.1094237467877973e-16,5.184065062397645,7.166991894732021,1.3155278076321379e-15,9.158550387173886,6.624609289603754,7.601511683739573e16
The columns in this file store the following information:
- "sequence": string; input filename of the sequence
- "sequence_index": integer; index of the sequence (always 1 if only a single sequence is used)
- "TR": integer; index of the repetition time that this data was acquired (between 1 and the value of
--nTR) - "readout": integer; index of the readout within a TR.
- "subset": integer; index of the subset of the total signal (e.g., intra-axonal) that has been output (see the
--subsetflag). The total signal will always be included with "subset" of 0. - "nspins": total number of spins contributing to the signal (might change for certain subsets of spins)
- "longitudinal": average longitudinal signal
- "transverse": average transverse signal
- "phase": average phase of the signal (in degrees)
- "Sx": signal strength in the x-direction
- "Sy": signal strength in the y-direction
- "Sz": signal strength in the z-direction
- "SESx", "SESy", "SE_Sz": estimated standard errors of the signal components
- "SNRSx", "SNRSy", "SNR_Sz": estimated signal-to-noise ratios of the signal components
The signal and magnetisation values in adaptive CLI output are reported per contributing spin. The nspins column reports the number of spins in the corresponding subset. If a finite --max-spins limit is reached before the target SNR, the program emits a warning.
The CLI has an --Nspins option for fixed-spin simulations as an alternative to defining the --target-snr. These options cannot be combined; use --max-spins to cap an adaptive run.
We can also output the signal of specific subsets of spins. For example, in the following we request to separately the output for just the spins inside the cylinders and just the spins outside of the cylinders.
mcmr run geometry.json dwi_te_80_bval_2.seq --target-snr 10 -o signal.csv --subset inside --subset outsidedf = 3×18 DataFrame
Row │ sequence sequence_index TR readout subset nspins longitudinal transverse phase Sx Sy Sz SE_Sx SE_Sy SE_Sz SNR_Sx SNR_Sy SNR_Sz
│ SubStrin… Int64 Int64 Int64 Int64 Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ dwi_te_80_bval_2.seq 1 1 1 0 100 -2.13163e-16 0.590382 90.7133 -0.00734991 0.590336 -2.13163e-16 4.94797 6.4281 1.17749e-15 11.9318 9.18439 8.49264e16
2 │ dwi_te_80_bval_2.seq 1 1 1 1 54 -2.17933e-16 0.999885 90.0847 -0.00147865 0.999884 -2.17933e-16 0.112371 0.00114376 8.67514e-16 480.494 47207.2 6.22468e16
3 │ dwi_te_80_bval_2.seq 1 1 1 2 46 -2.07563e-16 0.110484 97.4065 -0.0142423 0.109562 -2.07563e-16 4.97588 4.65713 8.03326e-16 1.02138 1.09129 5.7262e16
┌ Info: # Timestep determination
│ Maximum timestep set by tortuosity constraint based on size of geometry to 0.01 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).
We can see two additional rows in the output. These new rows are the in same order as the --subset flags provided to mcmr run and can be distinguished based on the "subset" column.
sequence,sequence_index,TR,readout,subset,nspins,longitudinal,transverse,phase,Sx,Sy,Sz,SE_Sx,SE_Sy,SE_Sz,SNR_Sx,SNR_Sy,SNR_Sz
dwi_te_80_bval_2.seq,1,1,1,0,100,-2.1316282072803005e-16,0.590381902490443,90.71331786859588,-0.007349914515239769,0.5903361496171936,-2.1316282072803005e-16,4.947972050304589,6.428104154012709,1.1774900437437752e-15,11.93179542018836,9.18438607006547,8.492640810962157e16
dwi_te_80_bval_2.seq,1,1,1,1,54,-2.1793266779679e-16,0.9998852270867836,90.08473011995912,-0.0014786493015213722,0.999884133758823,-2.1793266779679e-16,0.11237134227051725,0.0011437612376951182,8.675141184239348e-16,480.49441407137675,47207.23214181781,6.224682555957136e16
dwi_te_80_bval_2.seq,1,1,1,2,46,-2.0756343503861623e-16,0.11048424370325331,97.40647683164175,-0.014242269331343833,0.10956242910310694,-2.0756343503861623e-16,4.975882940636219,4.657125918501693,8.033256164698305e-16,1.0213815861391284,1.0912900572773734,5.726196084987857e16
We can see in the second row that inside the cylinders the transverse signal is very close to the 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 (i.e., the transverse signal is much lower in the third row). The reported values are per-spin means, so the first row is not the sum of the next two. To recover total signals, multiply each row by its nspins value.
A more complete state of all the spins can be produced using the --output-snapshot flag. For example, the command
mcmr run geometry.json dwi_te_80_bval_2.seq --output-snapshot snapshot.csv┌ Info: # Timestep determination
│ Maximum timestep set by tortuosity constraint based on size of geometry to 0.01 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).
will produce a file named "snapshot.csv" with:
sequence,sequence_index,TR,readout,spin,x,y,z,longitudinal,transverse,phase,Sx,Sy
dwi_te_80_bval_2.seq,1,1,1,1,-95.43492248597597,-110.07693499870821,174.93225962814844,-2.220446049250313e-16,1.0,88.32535756017796,0.029223863072789114,0.9995728917028027
dwi_te_80_bval_2.seq,1,1,1,2,-353.20003521145736,103.31937288472888,20.662280542099253,-2.220446049250313e-16,1.0,127.9184466902334,-0.6145392182034071,0.7888862714548562
dwi_te_80_bval_2.seq,1,1,1,3,-461.69264006555454,394.77883477900133,-536.191907581498,-2.220446049250313e-16,1.0,90.44935795923084,-0.007842695508475514,0.9999692455906638
dwi_te_80_bval_2.seq,1,1,1,4,123.51340678190198,146.71984534518168,70.10240794604688,-2.220446049250313e-16,1.0,86.19758147277093,0.06631601880134028,0.9977986698980613
...Each row corresponds to the state of a single spin. In addition to the sequence, sequence index, TR, and readout columns listed above, snapshot output contains:
- "spin": integer; index of the spin
- "x"/"y"/"z": floats; position of the spin at the time of the readout
- "longitudinal": longitudinal magnetisation
- "transverse": magnitude of the transverse magnetisation
- "phase": phase of the transverse magnetisation in degrees
- "Sx"/"Sy": transverse magnetisation components
The readout times can be adjusted using the --nTR, --times, and --skip-TR flags.