Balanced Hybrid-State Free Precession Pulse Sequence

This section explains the interface for calculating the spin evolution during a train of RF pulses, assuming balanced gradient moments (cf. Hybrid-state free precession in nuclear magnetic resonance). For this, we need the following packages:

using MRIgeneralizedBloch
using MAT
using Plots

and we use the pulse train described in the paper Rapid quantitative magnetization transfer imaging: utilizing the hybrid state and the generalized Bloch model:

control = matread(normpath(joinpath(pathof(MRIgeneralizedBloch), "../../docs/control_MT_v3p2_TR3p5ms_discretized.mat")))
α   = control["alpha"]
TRF = control["TRF"]
grad_moment = [i == 1 ? :crusher : :balanced for i ∈ eachindex(α)];

The vector grad_moment defines the gradient moments for each TR. Further, we define

TR = 3.5e-3; # s

The control has the following shape:

t = TR .* (1:length(TRF))
p1 = plot(t, α/π, ylabel="α/π", label=:none)
p2 = plot(t, TRF, xlabel="t (s)", ylabel="TRF (s)", label=:none)
p = plot(p1, p2, layout=(2,1))

We simulate the signal for the following biophysical parameters:

M0 = 1
m0s = 0.15
R1f = 0.5   # 1/s
R2f = 15    # 1/s
Rex = 30    # 1/s
R1s = 3     # 1/s
T2s = 10e-6 # s
ω0 = 0      # rad/s
B1 = 1;     # in units of B1_nominal

For speed purposes, it is advisable to use the linear approximation of the generalized Bloch model, which requires a precomputed $R_2^{s,l}$

R2slT = precompute_R2sl();

Now we have everything set up to calculate the signal:

s_linapp, _ = simulate_linearapprox(α, TRF, TR, ω0, B1, M0, m0s, R1f, R2f, Rex, R1s, T2s, R2slT; grad_moment)

(ComplexF64[-0.355900208848632 - 0.0im, -0.2587967919638072 + 0.0im, -0.1419918348574934 + 0.0im, -0.06842108345163278 + 0.0im, -0.025059481062928022 + 0.0im, -0.024838885366256218 + 0.0im, -0.007136448918897868 + 0.0im, -0.01533693523425317 + 0.0im, -0.0034896625579282727 + 0.0im, -0.012787743351199645 + 0.0im  …  0.009571162811461562 + 0.0im, 0.01015200266899922 + 0.0im, 0.011462673073508452 + 0.0im, 0.014501245088328337 + 0.0im, 0.021284663250329967 + 0.0im, 0.035816661379313164 + 0.0im, 0.07337672133106851 + 0.0im, 0.17203648560756177 + 0.0im, 0.2838939689166648 + 0.0im, 0.3750841419498616 + 0.0im], nothing)

By default, the output is a complex valued array where each element describes the transversal magnetization $x^f + i y^f$ of the free spin pool in each $T_\text{R}$. With $ω_0 = 0$, however, the imaginary part of the signal vanishes:

p = plot(xlabel="t (s)", ylabel="signal (normalized)"; legend=:topleft)
plot!(p, t, real.(vec(s_linapp)), label="Re(s); lin. approx.")
plot!(p, t, imag.(vec(s_linapp)), label="Im(s); lin. approx.")

For comparison, we can also solve the full integro-differential equation (IDE) for each RF pulse, which is more accurate, but much slower:

s_ide, _ = simulate_gbloch_ide(α, TRF, TR, ω0, B1, M0, m0s, R1f, R2f, Rex, R1s, T2s)

plot!(p, t, real.(vec(s_ide)), label="Re(s); IDE")
plot!(p, t, imag.(vec(s_ide)), label="Im(s); IDE")

Clicking on the legend entries allows to select and de-select individual graphs.

Real-valued magnetization vector

As an alternative to the complex-valued signal, we can also calculate the full magnetization vector $(x^f, y^f, z^f, x^s, z^s, 1)$ by supplying the keyword argument output=:realmagnetization. Here, $x$, $y$, $z$ denote the dimensions in space, the superscripts $f$ and $s$ denote the free and the semi-solid spin pool, respectively. We neglect the $y^s$ component, assuming (without loss of generality) ωₓ = 0 and given that $R_2^{s,l} \gg ω_0$.

m_linapp, _ = simulate_linearapprox(α, TRF, TR, ω0, B1, M0, m0s, R1f, R2f, Rex, R1s, T2s, R2slT;
    output=:realmagnetization, grad_moment)

p = plot(xlabel="t (s)", ylabel="magnetization (normalized)"; legend=:topleft)
plot!(p, t, [m_linapp[i][1] for i ∈ axes(m_linapp,1)] ./ (1 - m0s), label="xᶠ / m₀ᶠ")
plot!(p, t, [m_linapp[i][2] for i ∈ axes(m_linapp,1)] ./ (1 - m0s), label="yᶠ / m₀ᶠ")
plot!(p, t, [m_linapp[i][3] for i ∈ axes(m_linapp,1)] ./ (1 - m0s), label="zᶠ / m₀ᶠ")
plot!(p, t, [m_linapp[i][4] for i ∈ axes(m_linapp,1)] ./      m0s , label="xˢ / m₀ˢ")
plot!(p, t, [m_linapp[i][5] for i ∈ axes(m_linapp,1)] ./      m0s , label="zˢ / m₀ˢ")

Gradients

The same interface can also be used to calculate the derivatives of the signal wrt. the biophysical parameters. One can specify any subset of derivatives in any order with a vector of identifier objects:

grad_list = (grad_M0(), grad_m0s(), grad_R1f(), grad_R2f(), grad_Rex(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1());

Calling the function simulate_linearapprox with the keyword argument grad_list and this vector

_, g_linapp = simulate_linearapprox(α, TRF, TR, ω0, B1, 1, m0s, R1f, R2f, Rex, R1s, T2s, R2slT; grad_list, grad_moment);

returns the derivatives in the specified order:

p = plot(xlabel="t (s)", ylabel="signal (normalized)"; legend=:topleft)
plot!(p, t, real.(g_linapp[:,1] .* M0 ), label="Re(∂s/∂M₀ )*M₀")
plot!(p, t, real.(g_linapp[:,2] .* m0s), label="Re(∂s/∂m₀ˢ)*m₀ˢ")
plot!(p, t, real.(g_linapp[:,3] .* R1f), label="Re(∂s/∂R₁ᶠ)*R₁ᶠ")
plot!(p, t, real.(g_linapp[:,4] .* R2f), label="Re(∂s/∂R₂ᶠ)*R₂ᶠ")
plot!(p, t, real.(g_linapp[:,5] .* Rex), label="Re(∂s/∂Rₓ )*Rₓ ")
plot!(p, t, real.(g_linapp[:,6] .* R1s), label="Re(∂s/∂R₁ˢ)*R₁ˢ")
plot!(p, t, real.(g_linapp[:,7] .* T2s), label="Re(∂s/∂T₂ˢ)*T₂ˢ")
plot!(p, t, real.(g_linapp[:,8] .* ω0 ), label="Re(∂s/∂ω₀ )*ω₀ ")
plot!(p, t, real.(g_linapp[:,9] .* B1 ), label="Re(∂s/∂B₁ )*B₁ ")

Note that the first row is always the signal itself, which is equivalent to ∂s/∂M₀, as this toolbox always assumes M₀ = 1.

Apparent R₁

Above code calculates separate derivatives for $R_1^f$ and $R_1^s$. Yet, many publications, including our own paper "Rapid quantitative magnetization transfer imaging: utilizing the hybrid state and the generalized Bloch model" assumes an apparent longitudinal relaxation rate $R_1^a = R_1^f = R_1^f$. The derivatives wrt. this apparent relaxation rate can be calculated with

R1a = 1 # 1/s
grad_list = (grad_M0(), grad_R1a())
_, g_linapp = simulate_linearapprox(α, TRF, TR, ω0, B1, M0, m0s, R1a, R2f, Rex, R1a, T2s, R2slT; grad_list, grad_moment)

p = plot(xlabel="t (s)", ylabel="signal (normalized)"; legend=:topleft)
plot!(p, t, real.(g_linapp[:,1] .* M0 ), label="Re(∂s/∂M₀)/M₀")
plot!(p, t, real.(g_linapp[:,2] .* R1a), label="Re(∂s/∂R₁ᵃ)*R₁ᵃ")

Note that R1a appears here twice in the arguments of the simulate_linearapprox in place of R1f and R1s.

This page was generated using Literate.jl.