Simulations of spin and gradient echo sequences¶

This notebook illustrates the use of MCMRSimulator.jl to model the interplay between diffusion, magnetisation transfer, and magnetic susceptibility in explaining the T2 relaxation in white matter tissue.

Installation instructions and documentation of the simulator can be found at https://open.win.ox.ac.uk/pages/ndcn0236/mcmrsimulator.jl/dev/. MCMRSimulator.jl v0.3 was used to run this notebook.

This notebook was used to prepare figures for an abstract submitted to ISMRM 2023 in Toronto.

View this notebook online at https://users.fmrib.ox.ac.uk/~ndcn0236/notebooks/gradient_spin_echo.html

Download this notebook from https://users.fmrib.ox.ac.uk/~ndcn0236/notebooks/gradient_spin_echo.ipynb

In [1]:
using MCMRSimulator 
using CairoMakie
import Random

Simulations of diffusion-T2 spectra¶

Choosing a scanner¶

First, we will select a scanner on which to run the sequence. This will set the B0 field strength as well as the maximum slew rate and gradient strength. For these simulations, only the B0 strength really matters (as it affects the size of the myelin-induced off-resonance field).

In [2]:
scanner = MCMRSimulator.Siemens_Prisma
Out[2]:
Scanner(3.0, 80.0, 200.0)

Creating a substrate¶

Axons will be simulated as infinitely long cylinders all pointing in the y-direction. The configuration of axons will repeat every 30 micrometers. Outer radii are drawn from a gamma distribution with a mean of 1 micrometer and a variance of 1 micrometer squared. The final distribution will have a density of around 70% and be non-overlapping.

In [3]:
Random.seed!(123)
(positions, radii) = random_positions_radii((30., 30.), 0.7, 2, mean=1., variance=1.);

Based on this set of non-overlapping circles we create myelinated axons as a set of annuli, where the inner radius is 0.7 of the outer radius (i.e., all axons have a g-ratio of 0.7).

In [4]:
function get_tissue(; MT=0., myelin=false)
    annuli(0.7 .* radii, radii; positions=positions, rotation=:y, repeats=(30., 30.), myelin=myelin, MT_fraction=MT)
end
Out[4]:
get_tissue (generic function with 1 method)

Here we plot the cylinders with the off-resonance field generated by the myelin:

In [5]:
pp = PlotPlane(:y, size=30.)
f, ax, p = plot_off_resonance(pp, get_tissue(myelin=true), ngrid=200)
plot!(ax, pp, get_tissue(myelin=true))
f
Out[5]:
In [6]:
pp_small = PlotPlane(:y, [7., 0., 7.], size=7.)
f, ax, p = plot_off_resonance(pp_small, get_tissue(myelin=true), ngrid=200)
plot!(ax, pp_small, get_tissue())
#Colorbar(f, p, label="off-resonance (ppm)")

# test particle
Random.seed!(1)
path = trajectory([7., 0., 7.], Simulation([], diffusivity=3., geometry=get_tissue(), timestep=1e-3), 0:1e-3:0.5)
plot!(pp_small, path)

f
Out[6]:

The distribution of the off-resonance field in the three compartments (intra-axonal, myelin and extra-axonal) is shown below. The myelin magnetic susceptibility has both an isotropic ($\chi_I=100 {\rm ppb}$) and anisotropic ($\chi_A=-100 {\rm ppb}$) components. Because of the anisotropic component there is a net field offset in the intra-axonal space with an opposite sign to the magnetic field offset within the myelin water. In the extra-axonal space the dipole-like field does not produce an net offset.

In [7]:
snap = Snapshot(10000)
o = [MCMRSimulator.off_resonance(get_tissue(myelin=true), p) for p in MCMRSimulator.position.(snap)]
l = isinside(get_tissue(), snap)
f = hist(o[l .== 0], label="extra")
hist!(o[l .== 1], label="myelin")
hist!(o[l .== 2], label="intra")
axislegend()
f
Out[7]:

Helper functions to plot the signal profiles¶

In [8]:
function plot_all_signal(ax, times_plot, signal_plot; func=transverse, norm=false, ast=false, tissue=get_tissue())
    label = isinside(tissue, signal_plot[1])
    
    transform(use, signal_trans) = begin
        x = func.(s[use] for s in signal_trans)
        norm_by = sum(use)
        if ast
            @. -times_plot / log(x / norm_by)
        elseif norm
            x / norm_by
        else
            norm_alt = func === phase ? 1. : length(use)
            x / norm_alt
        end
    end
    plot_line!(use, label; linestyle=:solid) = lines!(ax, times_plot, transform(use, signal_plot), label=label, linestyle=linestyle)
    plot_line!(label .== label, "total")
    for (l_ind, l_text) in [(0, "extra"), (1, "myelin"), (2, "intra")]
        use = label .== l_ind
        plot_line!(use, l_text, linestyle=:dash)
    end
    plot_line!(label .!= 1, "no myelin")
end
Out[8]:
plot_all_signal (generic function with 1 method)
In [9]:
function plot_all_signal(times_plot, signal_plot; xlog=false, ylog=false, kwargs...)
    f = Figure(resolution=(400, 300))
    ax = Axis(f[1, 1], xscale=xlog ? log10 : identity, yscale=ylog ? log10 : identity, )
    plot_all_signal(ax, times_plot, signal_plot; kwargs...)
    f
end
Out[9]:
plot_all_signal (generic function with 2 methods)

Gradient echo sequence¶

The gradient echo sequence has a single 90 degree excitation pulse. The repetition time (TR) does not matter for these simulations.

In [10]:
ge_sequence = Sequence(TR=1000, pulses=[RFPulse(flip_angle=90)], scanner=scanner)
Out[10]:
Sequence{1, RFPulse, LinearGradients}(Scanner(3.0, 80.0, 200.0), LinearGradients([0.0, 1000.0], StaticArraysCore.SVector{3, Float64}[[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]], [0.0, 0.0, 0.0]), RFPulse[RFPulse(0.0, 1.5707963267948966, 6.123233995736766e-17, 1.0, 0.0, 1.0, 0.0)], 1000.0)

We create a simulation based on the gradient echo sequence setting the following parameters:

  • the diffusivity is set to the value for free water (3 micrometer^2/millisecond)
  • this random diffusion is constrained by the geometry set by get_tissue. The myelin in the tissue with generate an off-resonance field (myelin=true), but there will be no magnetisation transfer (MT=0.)
  • The relaxation rates are set to values appropriate for the CSF at 3T
    • T2 = 2000 milliseconds (Spijkerman, J.M. et al. (2018) ‘T2 mapping of cerebrospinal fluid: 3 T versus 7 T’, Magma (New York, N.y.), 31(3), pp. 415–424. doi:10.1007/s10334-017-0659-3.)
    • T1 = 4000 milliseconds (Yamashiro, A., Kobayashi, M. and Saito, T. (2019) ‘Cerebrospinal fluid T1 value phantom reproduction at scan room temperature’, Journal of Applied Clinical Medical Physics, 20(7), pp. 166–175. doi:10.1002/acm2.12659.) <- like the sequence TR, this will not affect the simulations shown here
In [11]:
gradient_echo = Simulation(ge_sequence, diffusivity=3., geometry=get_tissue(myelin=true, MT=5e-3), R1=1/4000, R2=1/2000, timestep=0.1);

The trajectory function returns the spin positions and orientations at the requested times (as defined in times_ge).

In [12]:
times_ge = 0.1:0.1:50
@time signal_ge = trajectory(1000, gradient_echo, times_ge);

  3.036572 seconds (3.06 M allocations: 228.077 MiB, 1.95% gc time, 28.68% compilation time)

We check that there is no exchange between compartments between the first and last timestep

In [13]:
label = isinside(get_tissue(), signal_ge[1])
label2 = isinside(get_tissue(), signal_ge[end])
all(label .== label2)
Out[13]:
true

The resulting signal profile. While the signal loss is small in individual compartments (dashed lines; magenta: intra; orange: extra; green:myelin), there is a large signal loss when adding the compartments together (light blue without myelin water; dark blue with myelin water).

In [14]:
f = plot_all_signal(times_ge, signal_ge, func=transverse)
plot!(ge_sequence)
f
Out[14]:
In [15]:
plot_all_signal(times_ge, signal_ge, func=transverse, ast=true, ylog=true)
Out[15]:

This rapid signal loss when the signal is summed across compartments is due to the dephasing of the various compartments.

In [16]:
plot_all_signal(times_ge, signal_ge, func=phase)
Out[16]:

Here we produce a grid of T2* estimates with readouts at various echo times with the myelin off-resonance field turned off (left) or on (right) and without (top) or with (bottom) magnetisation transfer. The dominant effect on the estimated T2* is the myelin off-resonance field, which leads to T2* estimates of around 20 ms. Note that this T2* is not the T2* of the myelin water, but rather the T2* due to the difference in the off-resonance field of the myelin water and the remaining tissue water.

In [17]:
function ge_grid_plot(loc; kwargs...)
    gradient_echo = Simulation(ge_sequence, diffusivity=3., geometry=get_tissue(; kwargs...), R1=1/4000, R2=1/2000, timestep=0.1)
    @time signal_ge = trajectory(1000, gradient_echo, times_ge)
    ax = Axis(loc, yscale=log10)
    plot_all_signal(ax, times_ge, signal_ge, func=transverse, ast=true)
    ylims!(ax, 8, 3000)
    ticks = [10, 30, 100, 300, 1000, 2000]
    ax.yticks = (ticks, string.(ticks))
end
Out[17]:
ge_grid_plot (generic function with 1 method)
In [18]:
f = Figure()
ge_grid_plot(f[1, 1], myelin=false, MT=0.)
ge_grid_plot(f[1, 2], myelin=true, MT=0.)
ge_grid_plot(f[2, 1], myelin=false, MT=5e-3)
ge_grid_plot(f[2, 2], myelin=true, MT=5e-3)
f

  1.527592 seconds (1.05 M allocations: 118.009 MiB, 2.70% gc time)
  1.957466 seconds (1.05 M allocations: 118.027 MiB)
  1.193477 seconds (1.05 M allocations: 118.016 MiB)
  1.970595 seconds (1.05 M allocations: 118.022 MiB)
Out[18]:

Spin echo sequence¶

The spin echo sequence is identical to the gradient echo sequence, but with the addition of an 180 degree refocus pulse.

In [19]:
se_sequence = Sequence(TR=1000, pulses=[RFPulse(flip_angle=90), RFPulse(time=25, flip_angle=180)], scanner=scanner)
spin_echo = Simulation(se_sequence, diffusivity=3., geometry=get_tissue(myelin=true, MT=5e-3), timestep=0.1, R1=1/4000, R2=1/2000);
In [20]:
times_se = 0.1:0.1:50
@time signal_se = trajectory(1000, spin_echo, times_se);

  2.640547 seconds (1.89 M allocations: 162.956 MiB, 1.83% gc time, 20.14% compilation time)
In [21]:
f = plot_all_signal(times_se, signal_se, func=transverse)
plot!(se_sequence)
f
Out[21]:
In [22]:
plot_all_signal(times_se, signal_se, func=phase)
Out[22]:

Modeling the signal evolution for a range of different spin echo sequence is a bit more difficult than for the gradient echo. For the spin echo, we need to simulate multiple sequences as for each echo time the refocus pulse will be at a different time. We now also add a Readout to the sequence to indicate at what time the signal from that particular sequence should be read out.

In [23]:
function se_grid_plot(loc; kwargs...)
    times_se = 1:50
    se_sequence = [Sequence(TR=1000, pulses=[RFPulse(flip_angle=90), RFPulse(time=t/2, flip_angle=180), Readout(time=t)]) for t in times_se]
    spin_echo = Simulation(se_sequence, diffusivity=3., geometry=get_tissue(; kwargs...), R1=1/4000, R2=1/2000)
    @time signal_se = [s[1] for s in readout(1000, spin_echo)]
    ax = Axis(loc, yscale=log10)
    plot_all_signal(ax, times_se, signal_se, func=transverse, ast=true)
    ylims!(ax, 8, 3000)
    ticks = [10, 30, 100, 300, 1000, 2000]
    ax.yticks = (ticks, string.(ticks))
end
Out[23]:
se_grid_plot (generic function with 1 method)
In [24]:
f = Figure()
se_grid_plot(f[1, 1], myelin=false, MT=0.)
se_grid_plot(f[1, 2], myelin=true, MT=0.)
se_grid_plot(f[2, 1], myelin=false, MT=5e-3)
se_grid_plot(f[2, 2], myelin=true, MT=5e-3)
f

  7.500952 seconds (5.77 M allocations: 445.426 MiB, 1.53% gc time, 92.32% compilation time)
  0.938407 seconds (938.04 k allocations: 184.965 MiB, 4.53% gc time)
  0.866871 seconds (937.98 k allocations: 184.963 MiB, 5.72% gc time)
  0.758667 seconds (937.99 k allocations: 184.963 MiB, 5.79% gc time)
Out[24]:

For the spin echo the myelin has a very small effect, because the dephasing between compartments is fully recoverable with only a minor dephasing within compartments remaining. Instead to explain the T2 differences between compartments we need to add a signal loss due to cross-relaxation or magnetisation transfer whenever the spin hits any membrane.

This can explain the T2 in intra-axonal water ($\sim$80 ms) as well as for the intra-axonal water ($\sim$50-60 ms), although the T2 of myelin water is significantly overestimated (should be $\sim$20 ms). This overestimation might simply reflect that the myelin configuration is actually better represented by a spiral than the two cylinders modeled here

In [25]:
times_se = 1:50
se_sequence = [Sequence(TR=1000, pulses=[RFPulse(flip_angle=90), RFPulse(time=t/2, flip_angle=180), Readout(time=t)], scanner=scanner) for t in times_se]
spin_echo = Simulation(se_sequence, diffusivity=3., geometry=get_tissue(; myelin=true, MT=5e-3), R1=1/4000, R2=1/2000)
@time signal_se = [s[1] for s in readout(1000, spin_echo)]
plot_all_signal(times_se, signal_se, func=transverse, ast=true)

  1.075943 seconds (998.78 k allocations: 188.279 MiB, 5.44% gc time, 3.21% compilation time)
Out[25]:

All axons having the same size¶

The increase of the T2 with echo time is due to the relatively larger contribution of the longer-T2 water in larger axons at the longer echo times.

In [26]:
Random.seed!(1)
(constant_positions, constant_radii) = random_positions_radii((30., 30.), 0.7, 2, variance=1e-8);
In [27]:
function constant_tissue(; MT=0., myelin=false)
    annuli(0.7 .* constant_radii, constant_radii; positions=constant_positions, rotation=:y, repeats=(30., 30.), myelin=myelin, MT_fraction=MT)
end
Out[27]:
constant_tissue (generic function with 1 method)
In [28]:
pp = PlotPlane(:y, size=30.)
f, ax, p = plot_off_resonance(pp, constant_tissue(myelin=true), ngrid=200)
plot!(ax, pp, constant_tissue(myelin=true))
f
Out[28]:
In [29]:
times_se = 1:50
se_sequence = [Sequence(TR=1000, pulses=[RFPulse(flip_angle=90), RFPulse(time=t/2, flip_angle=180), Readout(time=t)], scanner=scanner) for t in times_se]
spin_echo = Simulation(se_sequence, diffusivity=3., geometry=constant_tissue(; myelin=true, MT=5e-3), R1=1/4000, R2=1/2000)
@time signal_se = [s[1] for s in readout(1000, spin_echo)]
plot_all_signal(times_se, signal_se, func=transverse, ast=true, tissue=constant_tissue())

  3.332547 seconds (1.85 M allocations: 233.531 MiB, 1.17% gc time, 65.99% compilation time)
Out[29]:

Varying parameters¶

In [30]:
"Plot gradient or spin echo signal at TE=30 ms for a range of geometries"
function plot_varying(x, geometries, spin_echo=false)
    sequence = (
        spin_echo ?
        Sequence(TR=1000, pulses=[RFPulse(flip_angle=90), RFPulse(time=15, flip_angle=180), Readout(time=30.)], scanner=scanner) :
        Sequence(TR=1000, pulses=[RFPulse(flip_angle=90), Readout(time=30.)], scanner=scanner)
    )

    @time signal = [evolve(1000,
            Simulation(sequence, geometry=geom, diffusivity=3., R1=1/4000, R2=1/2000),
            30.,
            ) for geom in geometries];
    function get_compartment(geometry, signal, index)
        label = isinside(geometry, signal)
        if index == -1
            if maximum(label) == 1
                return signal
            else
                return signal[label .!= 1]
            end
        end
        if maximum(label) == 1
            # no myelin
            if index == 2
                index = 1
            elseif index == 1
                index = 2
            end
        end
        signal[label .== index]
    end
    function mytrans(s)
        length(s) == 0 ? NaN : transverse(s) / length(s)
    end

    f = Figure(resolution=(400, 300))
    ax = Axis(f[1, 1])
    for (compartment, ls) in [(nothing, :solid), (0, :dash), (1, :dash), (2, :dash), (-1, :solid)]
        part_signal = isnothing(compartment) ? signal : get_compartment.(geometries, signal, compartment)
        attenuation = mytrans.(part_signal)
        lines!(ax, x, -30 ./ log.(attenuation), linestyle=ls)
    end
    ylims!(0, 250)
    f
end
Out[30]:
plot_varying

Myelination (g-ratio)¶

In [31]:
grat = 0.6:0.1:1
geometries = TransformObstruction[annuli(0.6 .* radii, (0.6 / g) .* radii; positions=positions, rotation=:y, repeats=(30., 30.), myelin=true, MT_fraction=5e-3) for g in grat[1:end-1]]
# g-ratio of 1 is represented by a set of unmyelinated cylinders rather than annuli
push!(geometries, cylinders(0.6 .* radii; positions=positions, rotation=:y, repeats=(30., 30.), MT_fraction=5e-3));
plot_varying(grat, geometries, false)

  2.511681 seconds (3.29 M allocations: 193.408 MiB, 1.40% gc time, 68.52% compilation time)
Out[31]:
In [32]:
plot_varying(grat, geometries, true)

  1.723285 seconds (1.99 M allocations: 121.078 MiB, 3.02% gc time, 61.09% compilation time)
Out[32]:

Axon diameter¶

Changing the axonal diameter also scales the size of the extra-axonal space (for a constant density)

In [33]:
cradii = [0.5, 1, 2, 3, 5]
constant_positions, constant_radii
geometries = [
    annuli(0.7 .* constant_radii .* r, constant_radii .* r; positions=constant_positions .* r, rotation=:y, repeats=(30. * r, 30. * r), myelin=true, MT_fraction=5e-3)
    for r in cradii
]
plot_varying(cradii, geometries, false)

  2.824155 seconds (859.84 k allocations: 56.056 MiB, 13.87% compilation time)
Out[33]:
In [34]:
plot_varying(cradii, geometries, true)

  2.054231 seconds (832.51 k allocations: 54.883 MiB, 2.51% gc time, 18.46% compilation time)
Out[34]:

Angle with B0¶

In [35]:
angle = 0:10:90
orientation = [[0., sin(deg2rad(a)), cos(deg2rad(a))] for a in angle]
geometries = [annuli(0.7 .* radii, radii; positions=positions, rotation=o, repeats=(30., 30.), myelin=true, MT_fraction=5e-3) for o in orientation]
plot_varying(angle, geometries, false)

  1.619376 seconds (871.93 k allocations: 64.913 MiB, 3.00% compilation time)
Out[35]:
In [36]:
plot_varying(angle, geometries, true)

  1.648471 seconds (932.85 k allocations: 68.787 MiB, 2.85% compilation time)
Out[36]:

Axonal density¶

In [37]:
densities = 0.5:0.1:0.8
geometries = TransformObstruction[]
for d in densities
    (p, r) = random_positions_radii((30., 30.), d, 2, mean=1., variance=1.);
    push!(geometries,
            annuli(0.7 .* r, r; positions=p, rotation=:y, repeats=(30., 30.), myelin=true, MT_fraction=5e-3)
         )
end
plot_varying(densities, geometries, false)

  2.064110 seconds (2.60 M allocations: 148.125 MiB, 76.30% compilation time)
Out[37]:
In [38]:
plot_varying(densities, geometries, true)

  1.309036 seconds (1.38 M allocations: 83.909 MiB, 61.59% compilation time)
Out[38]: