Continuous Wave Simulation

The following code replicates the continuous wave simulation of Fig. 2 and is slightly more comprehensive in the sense that all discussed models are simulated.

For these simulations we need the following packages:

using MRIgeneralizedBloch
using DifferentialEquations
using QuadGK
using Plots

and we simulate an isolated semi-solid spin pool with the following parameters:

R₁ = 1.0 # 1/s
T₂ˢ = 10e-6 # s

Tʳᶠ = 2e-3 # s
ω₁ = 2000π # rad/s
ω₀ = 200π # rad/s

t = range(0, Tʳᶠ, length=1001) # time points for plotting
tspan = (0.0, Tʳᶠ); # simulation range

These parameters correspond to Fig. 2b, the parameters for replicating Fig. 2a are ω₁ = 200π and Tʳᶠ = 1s.

Lorentzian Lineshape

In this script, we simulate the three lineshapes separately, starting with the Lorentzian lineshape for which the Bloch model provides a ground truth.

\[\partial_t \begin{pmatrix} x \\ y \\ z \\ 1 \end{pmatrix} = \begin{pmatrix} -R_2 & -ω_0 & ω_1 & 0 \\ ω_0 & -R_2 & 0 & 0 \\ -ω_1 & 0 & -R_1 & R_1 \\ 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \\ 1 \end{pmatrix} ,\]

where the matrix is the Hamiltonian of the Bloch model. For a constant $ω_0$ and $ω_1$, we can evaluate the Bloch model by taking the matrix exponential of its Hamiltonian:

H(ω₁, ω₀, R₂, R₁) = [-R₂  -ω₀  ω₁  0;
                       ω₀ -R₂   0  0;
                      -ω₁   0 -R₁ R₁;
                        0   0   0  0]

z_Bloch = similar(t)
for i = 1:length(t)
    (_, _, z_Bloch[i], _) = exp(H(ω₁, ω₀, 1 / T₂ˢ, R₁) * t[i]) * [0; 0; 1; 1]
end

Henkelman's Steady-State Solution

When assuming an isolated semi-solid pool, Eq. (9) in Henkelman, R. Mark, et al. "Quantitative interpretation of magnetization transfer." Magnetic resonance in medicine 29.6 (1993): 759-766 reduces to

g_Lorentzian(ω₀) = T₂ˢ / π / (1 + (T₂ˢ * ω₀)^2)
z_steady_state_Lorentzian = R₁ / (R₁ + π * ω₁^2 * g_Lorentzian(ω₀))

0.0025267290589107306

where g_Lorentzian(ω₀) is the Lorentzian lineshape.

Graham's Single Frequency Approximation

The lineshape is also used to calculate Graham's single frequency approximation, which describes an exponential decay with the RF-induced saturation rate Rʳᶠ:

Rʳᶠ = π * ω₁^2 * g_Lorentzian(ω₀)
z_Graham_Lorentzian = @. (Rʳᶠ * exp(-t * (R₁ + Rʳᶠ)) + R₁) / (R₁ + Rʳᶠ);

Sled's Model

Sled's model is given by the ordinary differential equation (ODE)

\[\partial_t z(t) = \left(-\pi \int_0^t G(t-τ) \omega_1(τ)^2 dτ \right) z(t) + R_1 (1-z),\]

where $G(t-τ)$ is the Green's function. The Hamiltonian of this ODE is implemented in apply_hamiltonian_sled! and we solve this ODE with the DifferentialEquations.jl package:

z₀ = [1.0] # initial z-magnetization
param = (ω₁, 1, ω₀, R₁, T₂ˢ, greens_lorentzian) # defined by apply_hamiltonian_sled!
prob = ODEProblem(apply_hamiltonian_sled!, z₀, tspan, param)
z_Sled_Lorentzian = solve(prob);

Generalized Bloch Model

The generalized Bloch model is an integro-differential equation (IDE) as it depends on z(τ) instead of z(t):

\[\partial_t z(t) = - ω_1(t) \int_0^t G(t,τ) ω_1(τ) z(τ) dτ + R_1 (1 - z(t)) .\]

For off-resonant RF-pulses with $ω_1 = ω_x + i ω_y$, it is given by

\[\partial_t z(t) = - ω_y(t) \int_0^t G(t,τ) ω_y(τ) z(τ) dτ - ω_x(t) \int_0^t G(t,τ) ω_x(τ) z(τ) dτ + R_1 (1 - z(t)) .\]

The Hamiltonian of the IDE is implemented in apply_hamiltonian_gbloch! and we can solve this IDE with the delay-differential equation (DDE) solver of the DifferentialEquations.jl package:

zfun(p, t) = [1.0] # initialize history function (will be populated with an interpolation by the DDE solver)

param = (ω₁, 1, ω₀, R₁, T₂ˢ, greens_lorentzian) # defined by apply_hamiltonian_gbloch!
prob = DDEProblem(apply_hamiltonian_gbloch!, z₀, zfun, tspan, param)
z_gBloch_Lorentzian = solve(prob);

Now that we have solved all five models, we can plot the solutions for comparison:

p = plot(xlabel="t [ms]", ylabel="zˢ(t)")
plot!(p, 1e3t, zero(similar(t)) .+ z_steady_state_Lorentzian, label="Henkelman's steady-state")
plot!(p, 1e3t, z_Graham_Lorentzian, label="Graham's model")
plot!(p, 1e3t, (hcat(z_Sled_Lorentzian(t).u...)'), label="Sled's model")
plot!(p, 1e3t, (hcat(z_gBloch_Lorentzian(t).u...)'), label="generalized Bloch model")
plot!(p, 1e3t, z_Bloch, label="Bloch model")

Zooming into the plot, reveals virtually perfect (besides numerical differences) agreement between Bloch and generalized Bloch model and subtle, but existing differences when compared to the other models. Choosing a longer T₂ˢ amplifies these differences.

Gaussian Lineshape

We can repeat these simulations (with the exception of the Bloch model) for the Gaussian lineshape:

g_Gaussian(ω₀) = T₂ˢ / sqrt(2π) * exp(-(T₂ˢ * ω₀)^2 / 2)
z_steady_state_Gaussian = R₁ / (R₁ + π * ω₁^2 * g_Gaussian(ω₀))

Rʳᶠ = π * ω₁^2 * g_Gaussian(ω₀)
z_Graham_Gaussian = @. (Rʳᶠ * exp(-t * (R₁ + Rʳᶠ)) + R₁) / (R₁ + Rʳᶠ)

param = (ω₁, 1, ω₀, R₁, T₂ˢ, greens_gaussian) # defined by apply_hamiltonian_sled!
prob = ODEProblem(apply_hamiltonian_sled!, z₀, tspan, param)
z_Sled_Gaussian = solve(prob)

prob = DDEProblem(apply_hamiltonian_gbloch!, z₀, zfun, tspan, param)
z_gBloch_Gaussian = solve(prob)

p = plot(xlabel="t [ms]", ylabel="zˢ(t)")
plot!(p, 1e3t, zero(similar(t)) .+ z_steady_state_Gaussian, label="Henkelman's steady-state")
plot!(p, 1e3t, z_Graham_Gaussian, label="Graham' model")
plot!(p, 1e3t, (hcat(z_Sled_Gaussian(t).u...)'), label="Sled's model")
plot!(p, 1e3t, (hcat(z_gBloch_Gaussian(t).u...)'), label="generalized Bloch model")

Super-Lorentzian Lineshape

And we can repeat these simulations (with the exception of the Bloch model) for the super-Lorentzian lineshape, which reveals the most pronounced deviations between the models due to the substantially slower decay of the Green's function:

g_superLorentzian(ω₀) = sqrt(2 / π) * T₂ˢ * quadgk(ct -> exp(-2 * (T₂ˢ * ω₀ / abs(3 * ct^2 - 1))^2) / abs(3 * ct^2 - 1), 0.0, sqrt(1 / 3), 1)[1]
z_steady_state_superLorentzian = R₁ / (R₁ + π * ω₁^2 * g_superLorentzian(ω₀))

Rʳᶠ = π * ω₁^2 * g_superLorentzian(ω₀)
z_Graham_superLorentzian = @. (Rʳᶠ * exp(-t * (R₁ + Rʳᶠ)) + R₁) / (R₁ + Rʳᶠ)

G_superLorentzian = interpolate_greens_function(greens_superlorentzian, 0, Tʳᶠ/T₂ˢ)

param = (ω₁, 1, ω₀, R₁, T₂ˢ, G_superLorentzian)
prob = ODEProblem(apply_hamiltonian_sled!, z₀, tspan, param)
z_Sled_superLorentzian = solve(prob)

prob = DDEProblem(apply_hamiltonian_gbloch!, z₀, zfun, tspan, param)
z_gBloch_superLorentzian = solve(prob)


p = plot(xlabel="t [ms]", ylabel="zˢ(t)")
plot!(p, 1e3t, zero(similar(t)) .+ z_steady_state_superLorentzian, label="Henkelman's steady-state")
plot!(p, 1e3t, z_Graham_superLorentzian, label="Graham's model")
plot!(p, 1e3t, (hcat(z_Sled_superLorentzian(t).u...)'), label="Sled's model")
plot!(p, 1e3t, (hcat(z_gBloch_superLorentzian(t).u...)'), label="generalized Bloch model")

This page was generated using Literate.jl.