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 the variable ∂m∂t. The function interface is written in a way that we can directly feed it into a differential equation solver of the DifferentialEquations.jl package.

For this example, we need the following packages:

using MRIgeneralizedBloch
using DifferentialEquations
using Plots

We simulate the dynamics of a coupled spin system with the following parameters:

m0s = 0.15
R1f = 0.3 # 1/s
R2f = 15 # 1/s
R1s = 2 # 1/s
T2s = 10e-6 # s
Rx = 30; # 1/s

and the thermal equilibrium of the magnetization m = [xf; yf; zf; zs; 1]:

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

during a rectangular RF-pulse with the flip angle and pulse duration

α = π
TRF = 100e-6; # s

Further, 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 to improve performance:

G = interpolate_greens_function(greens_superlorentzian, 0, TRF / T2s);

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)$ for $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, 0) = m0. Here, we use a slightly more complicated initialization that allows us to index the history function in apply_hamiltonian_gbloch!, which improves performance:

mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? m0[idxs] : m0;

With this, we are ready to formulate and solve the differential equation:

param = (α/TRF, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G) # defined by apply_hamiltonian_gbloch!
prob = DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), param)
sol = solve(prob)

retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 12-element Vector{Float64}:
 0.0
 6.467225550500944e-8
 7.113948105551039e-7
 3.889778302843712e-6
 1.0349991481452107e-5
 1.9482147235515678e-5
 3.113999790802603e-5
 4.558766309671583e-5
 6.165620037675622e-5
 7.85495587736603e-5
 9.933867018921202e-5
 0.0001
u: 12-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.85, 0.15, 1.0]
 [0.0017269759778720242, 0.0, 0.8499982456170496, 0.14999969040313774, 1.0]
 [0.018995075304559135, 0.0, 0.8497877291258531, 0.1499625434748019, 1.0]
 [0.10360949274847489, 0.0, 0.8436614598386818, 0.14888421807808372, 1.0]
 [0.2715159022384952, 0.0, 0.805463505995969, 0.14227471329634592, 1.0]
 [0.48829308063593013, 0.0, 0.6957166284325105, 0.12426853976802493, 1.0]
 [0.7049489432221873, 0.0, 0.4747340290822344, 0.09160570833204726, 1.0]
 [0.8415621387574719, 0.0, 0.1176178685562197, 0.04646861650187269, 1.0]
 [0.7932878482861954, 0.0, -0.3039964508672744, 0.0023665474216662333, 1.0]
 [0.5301405392764524, 0.0, -0.6636236071528878, -0.028988822278229506, 1.0]
 [0.01771090612228332, 0.0, -0.8491007617900579, -0.044530765659317795, 1.0]
 [6.712405280141429e-5, 0.0, -0.8492833364736676, -0.04465901778432846, 1.0]

The plot function is implemented for such solution objects and we can plot the solution simply with

p = plot(sol, labels=["xᶠ" "yᶠ" "zᶠ" "zˢ" "1"], xlabel="t [s]", ylabel="m(t)")

More details on 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.