Simulation of a Single RF Pulse

The core of generalized Bloch model is implemented in the function apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p, t), which calculates the derivative ∂m/∂t for a given magnetization vector m and stores it in-place in the variable ∂m∂t. The function interface is written in a format that can be fed directly into a differential equation solver of the DifferentialEquations.jl package.

We need the following packages for this tutorial:

using MRIgeneralizedBloch
using DifferentialEquations
using SpecialFunctions
using QuadGK
using Plots

and we define the properties of a coupled spin system:

m0s = 0.15
R1f = 0.5 # 1/s
R2f = 13 # 1/s
R1s = 3 # 1/s
T2s = 12e-6 # s
Rx = 17; # 1/s

For most parts of this tutorial, we assume a perfectly calibrated, on-resonant RF-pulse:

B1 = 1
ω0 = 0; # rad/s

as well as a super-Lorentzian lineshape. We interpolate the corresponding Green's function in the range TRF ∈ [0, 1000 ⋅ T2s] to improve performance:

G = interpolate_greens_function(greens_superlorentzian, 0, 1000);

Rectangular RF-Pulses

First, we simulate a rectangular RF-pulse with a constant ω1:

α = π # rad
TRF = 500e-6 # s
ω1 = α/TRF # rad/s

6283.185307179586

Isolated Semi-Solid Spin Pool

The first example demonstrates how to simulate an isolated semi-solid spin pool for which the magnetization vector is defined by m = [zs; 1]. The appended 1 facilitates a more compact implementation of longitudinal relaxation to a non-zero thermal equilibrium. Here, we initialize the magnetization with the thermal equilibrium:

m0 = [1; 1];

The generalized Bloch model is a so-called integro-differential equation where the derivative $∂m/∂t$ at the time $t_1$ does not just depend on $m(t_1)$, but on $m(t)$ with $t \in [0, t_1]$. This is solved with a delay differential equation (DDE) solver that stores an interpolated history function mfun(p, t), which we use in the apply_hamiltonian_gbloch! function to evaluate the integral. This history function has to be initialized with

mfun(p, t) = m0;

For slight performance improvements, we could also initialize the history function with mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? m0[idxs] : m0. This syntax allows for direct indexing of the history function in apply_hamiltonian_gbloch!, which improves performance. Following the syntax of the DifferentialEquations.jl package, we can define and solve the differential equation:

param = (ω1, B1, ω0, R1s, T2s, G) # defined by apply_hamiltonian_gbloch!
prob = DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), param)
z_gBloch = solve(prob)

retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 6-element Vector{Float64}:
 0.0
 9.999999999999999e-5
 0.0001723322143021384
 0.0003139133745598435
 0.00048120310025674554
 0.0005
u: 6-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [0.8787349423147068, 1.0]
 [0.743794709072904, 1.0]
 [0.48522531688005804, 1.0]
 [0.24137682860240595, 1.0]
 [0.21901321820163858, 1.0]

The function apply_hamiltonian_gbloch! is implemented such that it infers from param = (ω1, B1, ω0, R1s, T2s, G) that you are only supplying the relaxation properties of the semi-solid spin pool and hence it simulates the spin dynamics of an isolated semi-solid spin pool. The DifferentialEquations.jl package also implements a plot function for the solution object, which can be used to display the result:

p = plot(z_gBloch, xlabel="t [s]", ylabel="zˢ(t)", idxs=1, label="g. Bloch")

For comparison, we can simulate the signal with Graham's spectral model, which describes an exponential saturation with the rate

f_ω1(t) = ω1
Rʳᶠ = graham_saturation_rate_spectral(ω0 -> lineshape_superlorentzian(ω0, T2s), f_ω1, TRF, ω0) # 1/s

2333.750703345905

The function graham_saturation_rate_spectral calculates the spectral power density S(ω₀,Δω,ω₁(t)) of the RF-pulse with an off-resonance frequency Δω and a pulse shape ω₁(t). The spectral power density is the squared absolute value of the pulse's Fourier transform, divided by the pulse duration. Thereafter, the function calculates the integral

\[<Rʳᶠ> = \int_{-∞}^{+∞} dω₀ S(ω₀,Δω,ω₁(t)) g(ω₀, T₂ˢ) .\]

where g(ω₀, T₂ˢ) is the spectral lineshape of the spin pool. Given this saturation rate, we can simply solve the ordinary differential equation ∂z/∂t = R₁ (1 - z) - Rʳᶠ z, which has the following analytical solution:

z_Graham(t) = (Rʳᶠ * exp(-t * (R1s + Rʳᶠ)) + R1s) / (R1s + Rʳᶠ)
plot!(p, z_Graham, 0, TRF, label="Graham")

This plot reveals a substantial difference between Graham's spectral model and the generalized Bloch model at the example of a 500μs inversion pulse, as, e.g., used in our in vivo experiments. This also entails a substantial difference in the z-magnetization of the semi-solid spin pool at the end of the RF-pulse:

z_gBloch(TRF)[1]

0.2190132182016387

z_Graham(TRF)

0.31175631470797394

Coupled Spin System

For a coupled spin system, the magnetization vector is defined as m = [xf; yf; zf; zs; 1] and the thermal equilibrium magnetization is given by:

m0 = [0; 0; 1-m0s; m0s; 1];

To indicate to the apply_hamiltonian_gbloch! function that we would like to simulate a coupled spin system, we simple provide it with the properties of both pools in the following format:

param = (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G);

Thereafter, we can use the same function calls as above to simulate the spin dynamics:

prob = DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), param)
m_gBloch = solve(prob)
p = plot(m_gBloch, xlabel="t [s]", ylabel="m(t)", idxs=1:4, labels=["xᶠ" "yᶠ" "zᶠ" "zˢ"])

Shaped RF-Pulses

The function apply_hamiltonian_gbloch! also allows for the simulation of RF-pulses with arbitrary shapes. To this end, ω₁(t) has to be defined as a function that takes time in seconds as an input and returns ω₁ at this particular point in time. For example, we can define a 1ms sinc-pulse:

TRF = 1e-3 # s
NSideLobes = 1
f_ω1(t) = sinc(2(NSideLobes+1) * t/TRF - (NSideLobes+1)) * α / (sinint((NSideLobes+1)π) * TRF/π / (NSideLobes+1))
p = plot(f_ω1, 0, TRF, xlabel="t [s]", ylabel="ω₁(t)", labels=:none)

NSideLobes defines here the number of side lobes on each side as can be seen in the plot. With numerical integration we can check if the RF-pulse has the correct flip angle:

quadgk(f_ω1, 0, TRF)[1] / α

1.0000000000000002

Isolated Semi-Solid Spin Pool

In order to calculate the spin dynamics of an isolated semi-solid spin pool during this shaped RF-pulse, we define the same tuple param as we did in Section Rectangular RF-Pulses with the only difference that the first element is a subtype of the abstract type Function:

param = (f_ω1, B1, ω0, R1s, T2s, G)
typeof(f_ω1) <: Function

true

With this definition of param, we can use the same function call as we did before:

m0 = [1; 1]
prob = DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), param)
z_gBloch = solve(prob)
p = plot(z_gBloch, xlabel="t [s]", ylabel="zˢ(t)", idxs=1, label="g. Bloch")

For comparison, we can simulate the signal with Graham's spectral model, which describes an exponential saturation with the rate

Rʳᶠ = graham_saturation_rate_spectral(ω0 -> lineshape_superlorentzian(ω0, T2s), f_ω1, TRF, ω0)

2386.0951117624327

z_Graham(t) = (Rʳᶠ * exp(-t * (R1s + Rʳᶠ)) + R1s) / (R1s + Rʳᶠ)
plot!(p, z_Graham, 0, TRF, label="Graham")

z_gBloch(TRF)[1]

-0.08069805953984102

z_Graham(TRF)

0.09285317751344668

The difference between Graham's model and the generalized Bloch model is more pronounced for this 1ms sinc-inversion pulse compared to the 500μs rectangular inversion pulse. We note, however, that the peak pulse amplitude of the sinc-pulse is higher and potentially too high clinical MRI systems.

Coupled Spin System

We can perform the same change to param to simulate a coupled spin system:

param = (f_ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G)
m0 = [0; 0; 1-m0s; m0s; 1]
prob = DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), param)
m_gBloch = solve(prob)
p = plot(m_gBloch, xlabel="t [s]", ylabel="m(t)", idxs=1:4, labels=["xᶠ" "yᶠ" "zᶠ" "zˢ"])

Double click on zˢ in the legend to isolate the semi-solid spin pool in the plot and compare the simulation to the last section.

ω₀-Sweep or Adiabatic RF-Pulses

The apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p, t) function is also implemented for RF-pulses with a varying RF frequency ω₀(t) as, e.g., used in adiabatic pulses. In order to simulate such pulses, the first element of param has to be ω₁(t)::Function like for Shaped RF-Pulses and, additionally, the third element has to be φ::Function instead of ω₀::Number. Notice two differences here: first, the (abstract) type Function instead of Number will tell the compiler to use the adiabatic-pulse implementation. Second, this implementation requires the phase of the RF-pulse as function of time φ(t) instead of the frequency because φ(t) ≠ ω₀(t) ⋅ t, if ω₀ is a function of time.

To demonstrate the interface at a practical example, we can defined a hyperbolic secant adiabatic inversion pulse:

TRF = 10.24e-3 # s
γ = 267.522e6 # gyromagnetic ratio in rad/s/T
ω₁ᵐᵃˣ = 13e-6 * γ # rad/s
μ = 5 # shape parameter in rad
β = 674.1 # shape parameter in 1/s

f_ω1(t) = ω₁ᵐᵃˣ * sech(β * (t - TRF/2)) # rad/s
f_ω0(t) = -μ * β * tanh(β * (t - TRF/2)) # rad/s
f_φ(t)  = -μ * log(cosh(β * t) - sinh(β*t) * tanh(β*TRF/2)); # rad

For this example, we analytically solved the integral $φ(t) = \int_0^t ω₀(τ) dτ$.

This pulse a hyperbolic secant amplitude:

p = plot(f_ω1, 0, TRF, xlabel="t [s]", ylabel="ω₁(t) [rad/s]", labels=:none)

and hyperbolic tangent frequency sweep

p = plot(f_ω0, 0, TRF, xlabel="t [s]", ylabel="ω₀(t) [rad/s]", labels=:none)

As explained above, we actually don't use the frequency in the implementation. Instead, we use the RF-phase:

p = plot(f_φ, 0, TRF, xlabel="t [s]", ylabel="φ(t) [rad]", labels=:none)

This interface, of course, also allows for the simulation of an isolated semi-solid spin pool with above described modifications to param. For brevity, however, we here directly simulate a coupled spin pool starting from thermal equilibrium:

m0 = [0, 0, 1-m0s, m0s, 1]
p = (f_ω1, B1, f_φ, m0s, R1f, R2f, Rx, R1s, T2s, G)
m_gBloch = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), p))
p = plot(m_gBloch, xlabel="t [s]", ylabel="m(t)", idxs=1:4, labels=["xᶠ" "yᶠ" "zᶠ" "zˢ"])

This simulation shows the intended inversion of the free spin pool and a near complete saturation of the semi-solid spin pool (double click on the corresponding legend entry). The transversal magnetization of the free pool exhibits some oscillations and at this point I should highlight a distinct difference between the implementation for adiabatic RF-pulses the above described implementation for constant ω₀: the latter case uses a frame of references that rotates with the RF-frequency about the z-axis, i.e. the RF-pulses rotate the magnetization with ω₁ about the y-axis and, additionally, the magnetizations rotates with ω₀ about the z-axis. The implementation of the adiabatic pulses uses a rotating frame of reference that is on resonance with the Larmor frequency of the spin isochromat and the phase of the RF-pulse rotates with ω₀(t). To simulate off-resonance, we can simply add a static value to above function or, more precisely, add a phase slope to φ:

Δω0 = 1000 # rad/s
f_φ_or(t) = f_φ(t) + Δω0 * t; # rad

We can, additionally, change B1 to demonstrate the robustness of adiabatic pulses:

B1 = 1.2 # 20% miss-calibration
p = (f_ω1, B1, f_φ_or, m0s, R1f, R2f, Rx, R1s, T2s, G)
m_gBloch = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), p))
p = plot(m_gBloch, xlabel="t [s]", ylabel="m(t)", idxs=1:4, labels=["xᶠ" "yᶠ" "zᶠ" "zˢ"])

While the spin dynamics during the pulse is changed, the final magnetization of the free pool is approximately the same compared to the on-resonant isochromat with B1 = 1. The final magnetization of the semi-solid spin pool is, like before, close to zero, but a close look reveals a small negative zˢ-magnetization (double click on zˢ in the plot's legend). The difference to the previous simulation highlights that the dynamics of semi-solid spin pool is not adiabatic and does depend on B₁.

Transversal Magnetization of the Semi-Solid Spin Pool

Throughout this tutorial, we only ever calculated and plotted the longitudinal magnetization of the semi-solid spin pool. This is foremost a result of way we formulate and solve the generalized Bloch equations (cf. Eq. (9) in the paper). But this implementation is also reflective of the standard use-case in magnetization transfer imaging, where we are foremost interested in the longitudinal magnetization of the semi-solid spin pool and its effect on the free spin pool. If required, however, it is easily possible to calculate the transversal magnetization with Eqs. (4-5) from the paper:

ωx(t) = -B1 * f_ω1(t) * sin(f_φ_or(t))
ωy(t) =  B1 * f_ω1(t) * cos(f_φ_or(t))
zs_gBloch(t) = m_gBloch(t)[4]
xs_gBloch(t) = quadgk(τ -> G((t - τ) / T2s) * ωx(τ) * zs_gBloch(τ), 0, t)[1]
ys_gBloch(t) = quadgk(τ -> G((t - τ) / T2s) * ωy(τ) * zs_gBloch(τ), 0, t)[1];

The last two lines calculate the numerical integral of the Green's function multiplied by the oscillating RF-fields. Similar code is also used in the implementation of apply_hamiltonian_gbloch!. Plotting these functions reveals the spin dynamics of the semi-solid spin pool:

p = plot(xs_gBloch, 0, TRF, xlabel="t [s]", ylabel="m(t)", label="xˢ")
plot!(p, ys_gBloch, 0, TRF, label="yˢ")
plot!(p, zs_gBloch, 0, TRF, label="zˢ")

More details about the interface, including the linear approximation of the generalized Bloch model can found in the following scripts that replicate all simulations, data analyses, and figures of the generalized Bloch paper.

This page was generated using Literate.jl.