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"];

and

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:

m0s = 0.15
R1f = 0.5   # 1/s
R2f = 15    # 1/s
Rx = 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 = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT)

1142×1×1 Array{ComplexF64, 3}:
[:, :, 1] =
    -0.3559004931993566 - 0.0im
    -0.2587969946905679 + 0.0im
   -0.14199195026591527 + 0.0im
   -0.06842113472542963 + 0.0im
  -0.025059503752018035 + 0.0im
  -0.024838902407844275 + 0.0im
  -0.007136457160852868 + 0.0im
  -0.015336945033751357 + 0.0im
 -0.0034896676844255983 + 0.0im
   -0.01278775140470097 + 0.0im
                        ⋮
   0.010152010471859227 + 0.0im
   0.011462681836658495 + 0.0im
   0.014501256117609435 + 0.0im
    0.02128467935644647 + 0.0im
    0.03581668834759564 + 0.0im
    0.07337677636041184 + 0.0im
    0.17203661446590307 + 0.0im
     0.2838941845821245 + 0.0im
    0.37508444162781684 + 0.0im

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 = calculatesignal_gbloch_ide(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, 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 = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT;
    output=:realmagnetization)

p = plot(xlabel="t (s)", ylabel="magnetization (normalized)"; legend=:topleft)
plot!(p, t, [m_linapp[i][1] for i=1:size(m_linapp,1)] ./ (1 - m0s), label="xᶠ / m₀ᶠ")
plot!(p, t, [m_linapp[i][2] for i=1:size(m_linapp,1)] ./ (1 - m0s), label="yᶠ / m₀ᶠ")
plot!(p, t, [m_linapp[i][3] for i=1:size(m_linapp,1)] ./ (1 - m0s), label="zᶠ / m₀ᶠ")
plot!(p, t, [m_linapp[i][4] for i=1:size(m_linapp,1)] ./      m0s , label="xˢ / m₀ˢ")
plot!(p, t, [m_linapp[i][5] for i=1:size(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_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1());

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

s_linapp = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT; grad_list);

returns the derivatives in the specified order:

p = plot(xlabel="t (s)", ylabel="signal (normalized)"; legend=:topleft)
plot!(p, t, real.(s_linapp[:,1,1]       ), label="Re(∂s/∂M₀ )*M₀")
plot!(p, t, real.(s_linapp[:,1,2] .* m0s), label="Re(∂s/∂m₀ˢ)*m₀ˢ")
plot!(p, t, real.(s_linapp[:,1,3] .* R1f), label="Re(∂s/∂R₁ᶠ)*R₁ᶠ")
plot!(p, t, real.(s_linapp[:,1,4] .* R2f), label="Re(∂s/∂R₂ᶠ)*R₂ᶠ")
plot!(p, t, real.(s_linapp[:,1,5] .* Rx ), label="Re(∂s/∂Rₓ )*Rₓ ")
plot!(p, t, real.(s_linapp[:,1,6] .* R1s), label="Re(∂s/∂R₁ˢ)*R₁ˢ")
plot!(p, t, real.(s_linapp[:,1,7] .* T2s), label="Re(∂s/∂T₂ˢ)*T₂ˢ")
plot!(p, t, real.(s_linapp[:,1,8] .* ω0 ), label="Re(∂s/∂ω₀ )*ω₀ ")
plot!(p, t, real.(s_linapp[:,1,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_R1a(),)
s_linapp = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1a, R2f, Rx, R1a, T2s, R2slT; grad_list)

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

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

This page was generated using Literate.jl.