Adiabatic RF Pulse

In this tutorial, we will build a hyperbolic-secant (HS) adiabatic inversion pulse. The key point is that KomaMRI can keep the RF frequency modulation explicit, then simulate the RF-frame dynamics and the resulting / robustness.

Defining a frequency-modulated pulse

Parameters follow the default pulse in Figure 1(c) of this JMRI adiabatic inversion example.

b1max = 13.5e-6
Trf = 18.3e-3
β̂ = 4
μ = 6
β = 2 * β̂ / Trf;

For an HS pulse, the hardest point for adiabatic following is near the center of the sweep, where is maximal and  . Requiring the RF precession rate to dominate that sweep rate gives the threshold used in the paper:

b1_threshold = β * sqrt(μ) / (2π * γ)
b1max > b1_threshold

true

First, we define the RF amplitude as a hyperbolic secant and the frequency sweep as a hyperbolic tangent.

t = range(-Trf / 2, Trf / 2, 201)
B1 = b1max .* sech.(β .* t)
Δf = -μ * β .* tanh.(β .* t) ./ (2π);

The argument of RF can be a scalar, for a constant RF offset such as a slice offset, or a waveform, as in this frequency-modulated pulse.

f = range(-2e3, 2e3, 161) |> collect
seq = Sequence()
@addblock seq += RF(B1, Trf, Δf, 0);

This representation is native to KomaMRI. Pulseq RF events do not store frequency-modulation waveforms, so exporting this sequence with write_seq requires first converting the sweep into RF phase samples. KomaMRI may do this conversion automatically in the future.

Plotting an adiabatic pulse

plot_seq can show both views. The default plot keeps explicit; the freq_in_phase keyword shows the same pulse after moving the frequency sweep into the RF phase.

p_freq = plot_seq(seq; max_rf_samples=Inf, slider=false, height=360, title="Frequency-modulated RF", showlegend=false)
p_phase = plot_seq(seq; freq_in_phase=true, max_rf_samples=Inf, slider=false, height=360, title="Phase-modulated RF", showlegend=false)

Keeping the frequency sweep explicit is useful because KomaMRI simulates a frequency-modulated RF pulse in the RF rotating frame. If the same pulse is first converted to sampled phase modulation, the simulation is more sensitive to the RF time sampling.

Comparing RF and rotating frames

Using a callback, we can record the magnetization during the RF pulse.

trajectory = NamedTuple[]
call_every_N_blocks = 1

record_traj = Callback(
    call_every_N_blocks,
    (progress_info, sim_blocks_info, device_data, sim_params) -> begin
        j = last(sim_blocks_info.parts[progress_info.block])
        push!(trajectory, (;
            Mxy=device_data.Xt.xy[1], Mz=device_data.Xt.z[1],
            ψ=device_data.seqd.ψ[j], B1=device_data.seqd.B1[j], Δf=device_data.seqd.Δf[j],
        ))
    end,
)

sim_params = KomaMRICore.default_sim_params()
sim_params["return_type"] = "state"
sim_params["max_rf_block_length"] = 1; # very inefficient; just for plots
obj0 = Phantom(; x=[0.0], Δw=[0.0]);
simulate(obj0, seq, sys; sim_params, callbacks=(record_traj,), verbose=false);

The blue vector is the normalized rotation axis , where   . With this sign convention, the magnetization precesses right-handed around the plotted axis.

HS adiabatic pulse B0 and B1 robustness

To showcase the off-resonance and robustness of this type of pulse, we can show its effect in a heatmap with     and   :

obj = Phantom(; x=zeros(length(f)), Δw=2π .* f);
b1_scales = range(0.05, 1.2, 47) |> collect;
sim_params = Dict{String, Any}("return_type" => "state");

Mz = map(b1_scales) do scale
    seq_scale = Sequence()
    @addblock seq_scale += RF(scale .* B1, Trf, Δf, 0)
    simulate(obj, seq_scale, sys; sim_params, verbose=false).z
end;

The dashed line is the analytic adiabatic threshold for this HS pulse, where we can see that the inversion is achieved after the threshold is surpassed.

---

This page was generated using Literate.jl.