API

CRBSourceTerms(dCRB_dxf, dCRB_dyf)

Callable struct storing the source terms ∂CRB/∂m for adjoint backpropagation. Call as source(iSeq, t, g) to get an 11-component SVector with non-zero entries at positions 6 (∂CRB/∂(∂xf/∂θ)) and 7 (∂CRB/∂(∂yf/∂θ)).

R2slInterpolants

Struct holding interpolated R2sl and its partial derivatives. Returned by precompute_R2sl. All fields are callable as f(TRF, α, B1, T2s).

Fields

• R2sl: R2sl(TRF, α, B1, T2s)
• dR2sl_dT2s: ∂R2sl/∂T2s
• dR2sl_dB1: ∂R2sl/∂B1
• dR2sl_dω1: ∂R2sl/∂ω1
• dR2sl_dTRF: ∂R2sl/∂TRF
• d2R2sl_dT2s_dω1: ∂²R2sl/(∂T2s ∂ω1)
• d2R2sl_dB1_dω1: ∂²R2sl/(∂B1 ∂ω1)
• d2R2sl_dT2s_dTRF: ∂²R2sl/(∂T2s ∂TRF)
• d2R2sl_dB1_dTRF: ∂²R2sl/(∂B1 ∂TRF)

costate = adjoint_backpropagate(crb_source, cycle_ops, propagators)

Compute the adjoint (co-state) variables by backpropagating the CRB source terms through the pulse sequence.

The adjoint equation is the time-reversed analog of the forward propagation, driven by the source terms ∂CRB/∂m(t,g). The computation has two phases:

Phase 1 — Accumulate boundary condition: Sum all source contributions backward through the sequence to find λ = Σ_t E[T→t]ᵀ · source(t), then solve the periodic boundary condition P[end] = Q⁻ᵀ · λ.

Phase 2 — Backward propagation: Starting from P[end], propagate the co-state backward: P[t-1] = E[t+1]ᵀ · P[t] + source(t), accumulating source terms at each step.

apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p, t)

Apply the generalized Bloch Hamiltonian to m and write the resulting derivative wrt. time into ∂m∂t.

Arguments

• ∂m∂t::Vector{Real}: Vector describing to derivative of m wrt. time; this vector has to be of the same size as m, but can contain any value, which is replaced by H * m
• m::Vector{Real}: Vector the spin ensemble state of the form [xf, yf, zf, zs, 1] if now gradient is calculated or of the form [xf, yf, zf, zs, 1, ∂xf/∂θ1, ∂yf/∂θ1, ∂zf/∂θ1, ∂zs/∂θ1, 0, ..., ∂xf/∂θn, ∂yf/∂θn, ∂zf/∂θn, ∂zs/∂θn, 0] if n derivatives wrt. θn are calculated
• mfun: History function; can be initialized with mfun(p, t; idxs=nothing) = typeof(idxs) <: Real ? 0.0 : zeros(5n + 5) for n gradients, and is then updated by the delay differential equation solvers
• p::NTuple{6,Any}: (ω1, B1, ω0, R1s, T2s, g) or
• p::NTuple{6,Any}: (ω1, B1, φ, R1s, T2s, g) or
• p::NTuple{10,Any}: (ω1, B1, ω0, m0s, R1f, R2f, Rex, R1s, T2s, g) or
• p::NTuple{10,Any}: (ω1, B1, φ, m0s, R1f, R2f, Rex, R1s, T2s, g) or
• p::NTuple{12,Any}: (ω1, B1, ω0, m0s, R1f, R2f, Rex, R1s, T2s, g, dG_o_dT2s_x_T2s, grad_list) or
• p::NTuple{12,Any}: (ω1, B1, ω0, m0s, R1f, R2f, Rex, R1s, T2s, g, dG_o_dT2s_x_T2s, grad_list) with the following entries
  • ω1::Real: Rabi frequency in rad/s (rotation about the y-axis) or
  • ω1(t)::Function: Rabi frequency in rad/s as a function of time for shaped RF-pulses
  • B1::Real: B1 scaling normalized so that B1=1 corresponds to a perfectly calibrated RF field
  • ω0::Real: Larmor or off-resonance frequency in rad/s or
  • φ::Function: RF-phase in rad as a function of time for frequency/phase-sweep pulses (works only in combination with ω1(t)::Function)
  • m0s::Real: Fractional semi-solid spin pool size in the range of 0 to 1
  • R1f::Real: Longitudinal spin relaxation rate of the free pool in 1/seconds
  • R2f::Real: Transversal spin relaxation rate of the free pool in 1/seconds
  • Rex::Real: Exchange rate between the two pools in 1/seconds
  • R1s::Real: Longitudinal spin relaxation rate of the semi-solid pool in 1/seconds
  • T2s::Real: Transversal spin relaxation time of the semi-solid pool in seconds
  • g::Function: Green's function of the form G(κ) = G((t-τ)/T2s)
  • dG_o_dT2s_x_T2s::Function: Derivative of the Green's function wrt. T2s, multiplied by T2s; of the form dG_o_dT2s_x_T2s(κ) = dG_o_dT2s_x_T2s((t-τ)/T2s)
  • grad_list::Vector{grad_param}: List of gradients to be calculated, i.e., any subset of [grad_m0s(), grad_R1f(), grad_R2f(), grad_Rex(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]; length of the vector must be n (cf. arguments m and ∂m∂t); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a()
• t::Real: Time in seconds

Optional:

• pulsetype=:normal: Use default for a regular RF-pulse; the option pulsetype=:inversion should be handled with care as it is only intended to calculate the saturation of the semi-solid pool and its derivative.

Examples

julia> using DelayDiffEq

julia> using DifferentialEquations

julia> α = π/2;

julia> TRF = 100e-6;

julia> ω1 = α/TRF;

julia> B1 = 1;

julia> ω0 = 0;

julia> m0s = 0.2;

julia> R1f = 1/3;

julia> R2f = 15;

julia> R1s = 2;

julia> T2s = 10e-6;

julia> Rex = 30;

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

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

julia> mfun(p, t; idxs=nothing) = typeof(idxs) <: Real ? 0.0 : zeros(5);

julia> sol = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), (ω1, B1, ω0, m0s, R1f, R2f, Rex, R1s, T2s, G)), MethodOfSteps(Tsit5()));

julia> dG_o_dT2s_x_T2s = interpolate_greens_function(dG_o_dT2s_x_T2s_superlorentzian, 0, TRF / T2s);

julia> grad_list = (grad_R2f(), grad_m0s());

julia> m0 = [0; 0; 1-m0s; m0s; 1; zeros(5*length(grad_list))];

julia> mfun2(p, t; idxs=nothing) = typeof(idxs) <: Real ? 0.0 : zeros(5 + 5*length(grad_list));

julia> sol = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun2, (0, TRF), (ω1, B1, ω0, m0s, R1f, R2f, Rex, R1s, T2s, G, dG_o_dT2s_x_T2s, grad_list)), MethodOfSteps(Tsit5()));

apply_hamiltonian_sled!(∂m∂t, m, p, t)

Apply Sled's Hamiltonian to m and write the resulting derivative wrt. time into ∂m∂t.

Arguments

- `∂m∂t::Vector{<:Real}`: Vector of length 1 describing to derivative of `m` wrt. time; this vector can contain any value, which is replaced by `H * m`
- `m::Vector{<:Real}`: Vector of length 1 describing the `zs` magnetization
- `p::NTuple{6 or 10, Any}`: `(ω1, B1, ω0, R1s, T2s, g)` for a simulating an isolated semi-solid pool or `(ω1, B1, ω0, m0s, R1f, R2f, Rex, R1s, T2s, g)` for simulating a coupled spin system; with
- `ω1::Real`: Rabi frequency in rad/s (rotation about the y-axis) or
- `ω1(t)::Function`: Rabi frequency in rad/s as a function of time for shaped RF-pulses
- `B1::Real`: B1 scaling normalized so that `B1=1` corresponds to a perfectly calibrated RF field
- `ω0::Real`: Larmor or off-resonance frequency in rad/s (is only used for the free spin pool)
- `R1f::Real`: Longitudinal spin relaxation rate of the free pool in 1/seconds
- `R2f::Real`: Transversal spin relaxation rate of the free pool in 1/seconds
- `R1s::Real`: Longitudinal spin relaxation rate of the semi-solid in 1/seconds
- `Rex::Real`: Exchange rate between the two pools in 1/seconds
- `T2s::Real`: Transversal spin relaxation time in seconds
- `g::Function`: Green's function of the form `G(κ) = G((t-τ)/T2s)`

• t::Real: Time in seconds

Examples

julia> using DifferentialEquations

julia> α = π/2;

julia> TRF = 100e-6;

julia> ω1 = α/TRF;

julia> B1 = 1;

julia> ω0 = 0;

julia> R1s = 2;

julia> T2s = 10e-6;

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

julia> m0 = [1];

julia> sol = solve(ODEProblem(apply_hamiltonian_sled!, m0, (0, TRF), (ω1, 1, ω0, R1s, T2s, G)), Tsit5());

fisher = assemble_fisher_matrix(magnetization)

Assemble the Fisher information matrix from the steady-state magnetization.

The signal derivative wrt. parameter g at time t has real part m[t,g][6] (= ∂xf/∂θg) and imaginary part m[t,g][7] (= ∂yf/∂θg). The FIM is the real symmetric matrix:

F[g1, g2] = Σ_{t,seq} (∂xf/∂θ_g1 · ∂xf/∂θ_g2 + ∂yf/∂θ_g1 · ∂yf/∂θ_g2)

x = bound_omega1_TRF!(ω1, TRF; ω1_min=zeros(size(ω1)), ω1_max=fill(2e3π, size(ω1)), TRF_min=fill(100e-6, size(TRF)), TRF_max=fill(500e-6, size(TRF)))

Bound the controls ω1 and TRF (over-written in place) and return a vector of length 2Npulses * NSeq with values in the range [-Inf, Inf] that relate to the bounded ω1 and TRF via the tanh function.

Arguments

• ω1: Control vector of length Npulses or matrix with the number of sequences in the second dimension
• TRF: Control vector of length Npulses or matrix with the number of sequences in the second dimension

Optional Keyword Arguments:

• ω1_min: elementwise lower bound for ω1 in rad/s
• ω1_max: elementwise upper bound for ω1 in rad/s
• TRF_min: elementwise lower bound for TRF in s
• TRF_max: elementwise bound for TRF in s

Examples

julia> ω1 = collect(range(0, 2000π, 100));

julia> TRF = collect(range(100e-6, 500e-6, 100));

julia> x = bound_omega1_TRF!(ω1, TRF)
200-element Vector{Float64}:
 -Inf
  -2.2924837393352853
  -1.9407818989717183
  -1.7328679513998637
  -1.5837912652403254
  -1.4669284349179519
  -1.3704200119626004
  -1.2879392139968635
  -1.2157089824185072
  -1.151292546497023
   ⋮
   1.2157089824185072
   1.2879392139968635
   1.3704200119626009
   1.4669284349179519
   1.5837912652403245
   1.7328679513998637
   1.9407818989717212
   2.2924837393352826
  Inf

build_crushed_propagator!(E, dEdω1, dEdTRF, t, g, ...)

Compute the propagator for a crushed (inversion) pulse at time index t and gradient parameter index g. Crushed pulses have zero derivatives wrt. ω₁ and TRF (they are not optimized).

E, dEdω1, dEdTRF = build_propagators(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rex, R1s, T2s, R2slT, grad_list; grad_moment)

Compute the per-pulse propagator E[t, g] and its derivatives dE/dω₁[t, g] and dE/dTRF[t, g] for each time step t and gradient parameter index g.

Each propagator maps the 11-component state vector from one TR to the next, incorporating free precession, RF pulse, phase cycling, and gradient spoiling.

Derivatives are computed via augmented 22×22 matrix exponentials: exp([H 0; dH/dθ H] · τ) yields both exp(H·τ) and d(exp(H·τ))/dθ in a single matrix exponential.

build_pulse_propagator!(E, dEdω1, dEdTRF, t, g, ...)

Compute the propagator and its ω₁/TRF derivatives for a single non-crushed RF pulse at time index t and gradient parameter index g, writing results into E[t,g], dEdω1[t,g], and dEdTRF[t,g].

The full TR propagator is: prop_freeprec · prop_pulse · prop_phasecycle · prop_freeprec, where prop_freeprec includes gradient spoiling and half-TR free precession.

grad_ω1, grad_TRF = control_gradients(costate, magnetization, propagators, dprop_dω1, dprop_dTRF)

Compute the gradients of the CRB wrt. the RF controls ω1 and TRF by evaluating the inner product ⟨P[t-1], dE[t]/dθ · m[t-1]⟩ at each time step, where P is the co-state (adjoint variable) and m is the magnetization.

CRB, crb_source = crb_and_derivatives(magnetization, weights)

Compute the weighted Cramér-Rao bound CRB = wᵀ diag(F⁻¹) and the source terms ∂CRB/∂m needed to drive the adjoint backpropagation.

Returns crb_source as a CRBSourceTerms callable struct, where crb_source(iSeq, t, g) returns an 11-component SVector with non-zero entries only at positions 6 (∂CRB/∂(∂xf/∂θ)) and 7 (∂CRB/∂(∂yf/∂θ)). These source terms use the chain rule through CRB = wᵀ diag(F⁻¹) and ∂F⁻¹/∂x = -F⁻¹ (∂F/∂x) F⁻¹.

CRB, grad_ω1, grad_TRF = crb_gradient(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rex, R1s, T2s, R2slT, grad_list, weights; grad_moment=...)

Calculate the Cramér-Rao bound (CRB) and its gradients wrt. the RF controls ω1 (amplitude) and TRF (duration) using the adjoint state method.

The CRB is computed as wᵀ diag(F⁻¹), where F is the Fisher information matrix assembled from the steady-state signal derivatives, and w are user-supplied weights for each fitted parameter.

Arguments

• ω1: RF amplitude per pulse in rad/s. Vector (single sequence) or matrix (columns = sequences).
• TRF: RF pulse duration per pulse in seconds. Same shape as ω1.
• TR::Real: Repetition time in seconds
• ω0::Real: Off-resonance frequency in rad/s
• B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
• m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
• R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
• R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
• Rex::Real: Exchange rate between the two spin pools in 1/seconds
• R1s::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
• T2s::Real: Transversal relaxation time of the semi-solid pool in seconds
• R2slT::NTuple{3, Function}: Tuple of three functions: R2sl(TRF, ω1, B1, T2s), dR2sldB1(TRF, ω1, B1, T2s), and R2sldT2s(TRF, ω1, B1, T2s). Can be generated with precompute_R2sl
• grad_list::Tuple{<:grad_param}: Tuple that specifies the gradients that are calculated; any subset/order of (grad_M0(), grad_m0s(), grad_R1f(), grad_R2f(), grad_Rex(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a(). Including grad_M0() is recommended to account for the equilibrium magnetization in the CRB.
• weights::transpose(Vector{Real}): Row vector of weights applied to the Cramér-Rao bounds of the individual parameters, matching grad_list in order and length.

Optional Keyword Arguments:

• grad_moment: Gradient spoiling scheme per pulse (:balanced, :crusher, :spoiler_dual, :spoiler_prepulse). :balanced simulates a TR with all gradient moments nulled. :crusher assumes equivalent (non-zero) gradient moments before and simulates the refocussing path of the extended phase graph. :spoiler_prepulse nulls all transverse magnetization before the RF pulse, emulating an idealized FLASH. :spoiler_dual nulls all transverse magnetization before and after the RF pulse.

Examples

julia> CRB, grad_ω1, grad_TRF = crb_gradient(range(pi/2, π, 100), range(100e-6, 400e-6, 100), 3.5e-3, 0, 1, 0.15, 0.5, 15, 30, 4, 10e-6, precompute_R2sl(), [grad_M0(), grad_m0s(), grad_R2f()], transpose([0, 1, 1]));

See also: Optimal Control

dG_o_dT2s_x_T2s_gaussian(κ)

Evaluate the derivative of Green's function, corresponding to a Gaussian lineshape, wrt. T2s at κ = (t-τ)/T2s and multiply it by T2s.

The multiplication is added so that the function merely depends on κ = (t-τ)/T2s. The actual derivative is given by dG_o_dT2s_x_T2s_gaussian((t-τ)/T2s)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> dGdT2s = dG_o_dT2s_x_T2s_gaussian((t-τ)/T2s)/T2s
1.9287498479639177e-15

dG_o_dT2s_x_T2s_lorentzian(κ)

Evaluate the derivative of Green's function, corresponding to a Lorentzian lineshape, wrt. T2s at κ = (t-τ)/T2s and multiply it by T2s.

The multiplication is added so that the function merely depends on κ = (t-τ)/T2s. The actual derivative is given by dG_o_dT2s_x_T2s_lorentzian((t-τ)/T2s)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> dGdT2s = dG_o_dT2s_x_T2s_lorentzian((t-τ)/T2s)/T2s
45.39992976248485

dG_o_dT2s_x_T2s_superlorentzian(κ)

Evaluate the derivative of Green's function, corresponding to a super-Lorentzian lineshape, wrt. T2s at κ = (t-τ)/T2s and multiply it by T2s.

The multiplication is added so that the function merely depends on κ = (t-τ)/T2s. The actual derivative is given by dG_o_dT2s_x_T2s_superlorentzian((t-τ)/T2s)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> dGdT2s = dG_o_dT2s_x_T2s_superlorentzian((t-τ)/T2s)/T2s
15253.09503367097

qM = fit_gBloch(data, α, TRF, TR;
    reM0 = (-Inf,   1,  Inf),
    imM0 = (-Inf,   0,  Inf),
    m0s  = (   0, 0.2,    1),
    R1f  = (   0, 0.3,  Inf),
    R2f  = (   0,  15,  Inf),
    Rex  = (   0,  20,  Inf),
    R1s  = (   0,   3,  Inf),
    T2s  = (8e-6,1e-5,12e-6),
    ω0   = (-Inf,   0,  Inf),
    B1   = (   0,   1,  1.5),
    R1a  = (   0, 0.7,  Inf),
    u=1,
    fit_apparentR1=false,
    show_trace=false,
    maxIter=100,
    R2slT = precompute_R2sl(TRF_min=minimum(TRF), TRF_max=maximum(TRF), T2s_min=minimum(T2s), T2s_max=maximum(T2s), ω1_max=maximum(α ./ TRF), B1_max=maximum(B1)),
    )

Fit the generalized Bloch model for a train of RF pulses and balanced gradient moments to data.

Arguments

• data::Vector{Number}: Array of measured data points, either in the time or a compressed domain (cf. u)
• α::Vector{Real}: Array of flip angles in radians; can also be a Vector{Vector{Real}} which simulates each RF pattern and concatenates the signals of each simulation
• TRF::Vector{Real}: Array of the RF-pulse durations in seconds (or Vector{Vector{Real}} if α::Vector{Vector{Real}}`)
• TR::Real: Repetition time in seconds
• ω0::Real: Off-resonance frequency in rad/s

Optional Keyword Arguments:

• reM0::Union{Real, Tuple{Real, Real, Real}}: Real part of M0; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• imM0::Union{Real, Tuple{Real, Real, Real}}: Imaginary part of M0; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• m0s::Union{Real, Tuple{Real, Real, Real}}: Fractional size of the semi-solid pool (should be in range of 0 to 1); either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• R1f::Union{Real, Tuple{Real, Real, Real}}: Longitudinal relaxation rate of the free pool in 1/s; only used in combination with fit_apparentR1=false; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• R2f::Union{Real, Tuple{Real, Real, Real}}: Transversal relaxation rate of the free pool in 1/s; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• Rex::Union{Real, Tuple{Real, Real, Real}}: Exchange rate between the two spin pools in 1/s; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• R1s::Union{Real, Tuple{Real, Real, Real}}: Longitudinal relaxation rate of the semi-solid pool in 1/s; only used in combination with fit_apparentR1=false; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• T2s::Union{Real, Tuple{Real, Real, Real}}: Transversal relaxation time of the semi-solid pool in s; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• ω0::Union{Real, Tuple{Real, Real, Real}}: Off-resonance frequency in rad/s; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• B1::Union{Real, Tuple{Real, Real, Real}}: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• R1a::Union{Real, Tuple{Real, Real, Real}}: Apparent longitudinal relaxation rate in 1/s; only used in combination with fit_apparentR1=true; either fixed value as a Real or fit limits thereof as a Tuple with the elements (min, start, max)
• u::Union{Number, Matrix}: Compression matrix that transform the simulated time series to a series of coefficients. Set to 1 by default to enable the fitting in the time domain
• fit_apparentR1::Bool: Switch between fitting R1f and R1s separately (false; default) and an apparent R1a = R1f = R1s (true)
• show_trace::Bool: print output during the optimization; default=false
• maxIter::Int: Maximum number of iteration; default=100
• R2slT::NTuple{3, Function}: Tuple of three functions: R2sl(TRF, ω1, B1, T2s), dR2sldB1(TRF, ω1, B1, T2s), and R2sldT2s(TRF, ω1, B1, T2s). By default generated with precompute_R2sl

Examples

c.f. Non-Linear Least Square Fitting

ω1, TRF = get_bounded_omega1_TRF(x)

Transform a vector of length 2Npulses with values in the range [-Inf, Inf] into two vectors of length Npulses, which describe the bounded controls ω1 and TRF.

Arguments

• x::Vector{Real}: Control vector of length = 2Npulses with values in the range [-Inf, Inf]

Optional Keyword Arguments:

• ω1_min::Vector{Real}: elementwise lower bound for ω1 in rad/s
• ω1_max::Vector{Real}: elementwise upper bound for ω1 in rad/s
• TRF_min::Vector{Real}: elementwise lower bound for TRF in s
• TRF_max::Vector{Real}: elementwise bound for TRF in s

Examples

julia> x = repeat(range(-1000.0, 1000.0, 100), 2);

julia> ω1, TRF = get_bounded_omega1_TRF(x)
([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  6283.185307179586, 6283.185307179586, 6283.185307179586, 6283.185307179586, 6283.185307179586, 6283.185307179586, 6283.185307179586, 6283.185307179586, 6283.185307179586, 6283.185307179586], [0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001  …  0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.0005, 0.0005])

graham_saturation_rate_single_frequency(lineshape, ω1, TRF, Δω)

Calculate saturation rate (in units of 1/s) according to Graham's single frequency approximation.

Arguments

• lineshape::Function: as a function of ω₀ (in rad/s). Supply, e.g., the anonymous function ω₀ -> lineshape_superlorentzian(ω₀, T2s). Note that the integral over the lineshape has to be 1.
• ω1::Function: ω1 in rad/s as a function of time (in units of s) where the puls shape is defined for t ∈ [0,TRF]
• TRF::Real: duration of the RF pulse in s
• Δω::Real: offset frequency in rad/s

Examples

julia> using SpecialFunctions

julia> T2s = 10e-6;

julia> α = π;

julia> TRF = 100e-6;

julia> NSideLobes = 1;

julia> ω1(t) = sinc(2(NSideLobes+1) * t/TRF - (NSideLobes+1)) * α / (sinint((NSideLobes+1)π) * TRF/π / (NSideLobes+1));

julia> Δω = 200;

julia> graham_saturation_rate_single_frequency(ω₀ -> lineshape_superlorentzian(ω₀, T2s), ω1, TRF, Δω)
419969.3376658947

graham_saturation_rate_spectral(lineshape, ω1, TRF, Δω)

Calculate saturation rate (in units of 1/s) according to Graham's spectral model.

Arguments

• lineshape::Function: as a function of ω₀ (in rad/s). Supply, e.g., the anonymous function ω₀ -> lineshape_superlorentzian(ω₀, T2s). Note that the integral over the lineshape has to be 1.
• ω1::Function: ω1 in rad/s as a function of time (in units of s) where the puls shape is defined for t ∈ [0,TRF]
• TRF::Real: duration of the RF pulse in s
• Δω::Real: offset frequency in rad/s

Examples

julia> using SpecialFunctions

julia> T2s = 10e-6;

julia> α = π;

julia> TRF = 100e-6;

julia> NSideLobes = 1;

julia> ω1(t) = sinc(2(NSideLobes+1) * t/TRF - (NSideLobes+1)) * α / (sinint((NSideLobes+1)π) * TRF/π / (NSideLobes+1));

julia> Δω = 200;

julia> graham_saturation_rate_spectral(ω₀ -> lineshape_superlorentzian(ω₀, T2s), ω1, TRF, Δω)
56135.388046022905

greens_gaussian(κ)

Evaluate the Green's function corresponding to a Gaussian lineshape at κ = (t-τ)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> greens_gaussian((t-τ)/T2s)
1.9287498479639178e-22

greens_lorentzian(κ)

Evaluate the Green's function corresponding to a Lorentzian lineshape at κ = (t-τ)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> greens_lorentzian((t-τ)/T2s)
4.5399929762484854e-5

greens_superlorentzian(κ)

Evaluate the Green's function corresponding to a super-Lorentzian lineshape at κ = (t-τ)/T2s.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> greens_superlorentzian((t-τ)/T2s)
0.14712468680944424

hamiltonian_linear(ω1, B1, ω0, T, M0, m0s, R1f, R2f, Rex, R1s, R2s[, dR2sdT2s, dR2sdB1, grad_type])

Calculate the hamiltonian of the linear approximation of the generalized Bloch model.

If no gradient is supplied, it returns a 6x6 (static) matrix with the dimensions (in this order) [xf, yf, zf, xs, zs, 1]; the attached 1 is a mathematical trick to allow for $T_1$ relaxation to a non-zero thermal equilibrium. If a gradient is supplied, it returns a 11x11 (static) matrix with the dimensions (in this order) [xf, yf, zf, xs, zs, dxf/dθ, dyf/dθ, dzf/dθ, dxs/dθ, dzs/dθ, 1], where θ is the parameter specified by grad_type

Arguments

• ω1::Real: Rabi frequency in rad/s (rotation about the y-axis)
• B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
• ω0::Real: Larmor (or off-resonance) frequency in rad/s (rotation about the z-axis)
• T::Real: Time in seconds; this can, e.g., be the RF-pulse duration, or the time of free precession with ω1=0
• M0::Real: Equilibrium magnetization (proton density) scaling factor
• m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
• R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
• R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
• Rex::Real: Exchange rate between the two spin pools in 1/seconds
• R1s::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
• R2s::Real: Transversal relaxation rate of the semi-solid pool in 1/seconds; this number can be calculated with the first function returned by precompute_R2sl to implement the linear approximation described in the generalized Bloch paper

Optional Arguments:

• dR2sdT2s::Real: Derivative of linearized R2sl wrt. the actual T2s; only required if grad_type = grad_T2s(); this number can be calculated with the second function returned by precompute_R2sl
• dR2sdB1::Real: Derivative of linearized R2sl wrt. B1; only required if grad_type = grad_B1(); this number can be calculated with the third function returned by precompute_R2sl
• grad_type::grad_param: grad_m0s(), grad_R1f(), grad_R1s(), grad_R2f(), grad_Rex(), grad_T2s(), grad_ω0(), or grad_B1(); create one hamiltonian for each desired gradient

Examples

julia> α = π;

julia> T = 500e-6;

julia> ω1 = α/T;

julia> B1 = 1;

julia> ω0 = 0;

julia> m0s = 0.1;

julia> R1f = 1;

julia> R2f = 15;

julia> Rex = 30;

julia> R1s = 6.5;

julia> R2s = 1e5;

julia> M0 = 1;

julia> m0 = [0, 0, M0*(1-m0s), 0, M0*m0s, 1];

julia> (xf, yf, zf, xs, zs, _) = exp(hamiltonian_linear(ω1, B1, ω0, T, M0, m0s, R1f, R2f, Rex, R1s, R2s)) * m0
6-element StaticArraysCore.SVector{6, Float64} with indices SOneTo(6):
  0.0010647535813058293
  0.0
 -0.8957848274535014
  0.005126529591877105
  0.08122007142111888
  1.0

interpolate_greens_function(f, κmin, κmax)

Interpolate the Green's function f in the range between κmin and κmax.

The interpolation uses the ApproxFun.jl package that incorporates Chebyshev polynomials and ensures an approximation to machine precision.

Examples

julia> t = 100e-6;

julia> τ = 0;

julia> T2s = 10e-6;

julia> greens_superlorentzian((t-τ)/T2s)
0.14712468680944424

julia> Gint = interpolate_greens_function(greens_superlorentzian, 0, 20);

julia> Gint((t-τ)/T2s)
0.14712468680944407

F = penalty_RF_power!(grad_ω1, grad_TRF, ω1, TRF; λ=1, Pmax=3e6, TR=3.5e-3)

Calculate RF power penalty and add the gradients in place.

Arguments

• grad_ω1::Vector{Real}: Gradient of control, which will be added in place (matrix if more than 1 sequence are optimized)
• grad_TRF::Vector{Real}: Gradient of control, which will be added in place (matrix if more than 1 sequence are optimized)
• ω1::Vector{Real}: Control vector (matrix if more than 1 sequence are optimized)
• TRF::Vector{Real}: Control vector (matrix if more than 1 sequence are optimized)

Optional Keyword Arguments:

• λ::Real: regularization parameter
• Pmax::Real: Maximum average power deposition in (rad/s)²; everything above this value will be penalized and with an appropriate λ, the resulting power will be equal to or less than this value.
• TR::Real: Repetition time of the pulse sequence

Examples

julia> ω1 = range(0, 4000π, 100);

julia> TRF = range(100e-6, 500e-6, 100);

julia> grad_ω1 = similar(ω1);

julia> grad_TRF = similar(ω1);

julia> F = penalty_RF_power!(grad_ω1, grad_TRF, ω1, TRF; λ=1e3, Pmax=3e6, TR=3.5e-3)
9.418321886730644e15

F = penalty_TRF_variation!(grad_TRF, TRF; λ=1, grad_moment=[i[1] == 1 ? :spoiler_dual : :balanced for i ∈ CartesianIndices(TRF)])

Calculate the total variation penalty of TRF and add to grad_TRF in place.

Arguments

• grad_TRF::Vector{Real}: Gradient of control, which will added in place (matrix if more than 1 sequence are optimized)
• TRF::Vector{Real}: Control vector (matrix if more than 1 sequence are optimized)

Optional Keyword Arguments:

• λ::Real: regularization parameter
• grad_moment = [i[1] == 1 ? :spoiler_dual : :balanced for i ∈ CartesianIndices(TRF)]: Different types of gradient moments of each TR are possible (:balanced, :spoilerdual, :crusher). Skip :crusher and :spoilerdual TRs for the TRF TV penalty

Examples

julia> TRF = range(100e-6, 500e-6, 100);

julia> grad_TRF = similar(TRF);

julia> F = penalty_TRF_variation!(grad_TRF, TRF; λ = 1e-3)
3.9595959595959597e-7

F = penalty_alpha_curvature!(grad_ω1, grad_TRF, ω1, TRF; λ=1, grad_moment=[i[1] == 1 ? :spoiler_dual : :balanced for i ∈ CartesianIndices(ω1)])

Calculate second order penalty of variations of the flip angle α and adds in place to the gradients.

Arguments

• grad_ω1::Vector{Real}: Gradient of control, which will be added in place (matrix if more than 1 sequence are optimized)
• grad_TRF::Vector{Real}: Gradient of control, which will be added in place (matrix if more than 1 sequence are optimized)
• ω1::Vector{Real}: Control vector (matrix if more than 1 sequence are optimized)
• TRF::Vector{Real}: Control vector (matrix if more than 1 sequence are optimized)

Optional Keyword Arguments:

• λ::Real: regularization parameter
• grad_moment = [i[1] == 1 ? :spoiler_dual : :balanced for i ∈ CartesianIndices(ω1)]: Different types of gradient moments of each TR are possible (:balanced, :spoilerdual, :crusher). Skip :crusher and :spoilerdual TRs second order penalty

Examples

julia> ω1 = range(0, 2000π, 100);

julia> TRF = range(100e-6, 500e-6, 100);

julia> grad_ω1 = similar(ω1);

julia> grad_TRF = similar(ω1);

julia> F = penalty_alpha_curvature!(grad_ω1, grad_TRF, ω1, TRF; λ = 1e-3)
0.005015194549476384

precompute_R2sl([;TRF_min=100e-6, TRF_max=500e-6, T2s_min=5e-6, T2s_max=21e-6, ω1_max=π/TRF_max, B1_max=1.4, greens=greens_superlorentzian])

Pre-compute and interpolate the linearized R2sl(TRF, α, B1, T2s) and its derivatives dR2sldB1(TRF, α, B1, T2s), R2sldT2s(TRF, α, B1, T2s) etc. in the range specified by the arguments.

The function solves the generalized Bloch equations of an isolated semi-solid pool for values in the specified range, calculates the linearized R2sl that minimizes the error of zs at the end of the RF-pulse, and interpolates between the different samples.

Optional Arguments:

• TRF_min::Real: lower bound of the RF-pulse duration range in seconds
• TRF_max::Real: upper bound of the RF-pulse duration range in seconds
• T2s_min::Real: lower bound of the T2s range in seconds
• T2s_max::Real: upper bound of the T2s range in seconds
• ω1_max::Real: upper bound of the Rabi frequency ω1, the default is the frequency of a 500μs long π-pulse
• B1_max::Real: upper bound of the B1 range, normalized so that B1 = 1 corresponds to a perfectly calibrated RF field
• greens=greens_superlorentzian: Greens function in the form G(κ) = G((t-τ)/T2s). This package supplies the three Greens functions greens=greens_superlorentzian (default), greens=greens_lorentzian, and greens=greens_gaussian

Examples

julia> R2slT = precompute_R2sl();


julia> R2sl, dR2sldB1, R2sldT2s, _ = precompute_R2sl(TRF_min=100e-6, TRF_max=500e-6, T2s_min=5e-6, T2s_max=15e-6, ω1_max=π/500e-6, B1_max=1.4, greens=greens_gaussian);

simulate_gbloch_ide(α, TRF, TR, ω0, B1, M0, m0s, R1f, R2f, Rex, R1s, T2s[; grad_list=nothing, Ncyc=2, output=:complexsignal])

Calculate the signal or magnetization evolution with the full generalized Bloch model assuming a super-Lorentzian lineshape (slow).

The simulation assumes a sequence of rectangular RF-pulses with varying flip angles α and RF-pulse durations TRF, but a fixed repetition time TR. Further, it assumes balanced gradient moments.

Always returns a tuple (signal, gradients) where signal is a vector (for output=:complexsignal) or matrix (for output=:realmagnetization) scaled by M0, and gradients is a matrix with one column per entry in grad_list, or nothing if no gradients are requested. M0 must be real-valued; for complex-valued M0 (e.g. to account for receive coil phase), simulate with M0=1 and multiply the complex M0 afterward.

Arguments

• α::Vector{Real}: Array of flip angles in radians
• TRF::Vector{Real}: Array of the RF-pulse durations in seconds
• TR::Real: Repetition time in seconds
• ω0::Real: Off-resonance frequency in rad/s
• B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
• M0::Real: Equilibrium magnetization (proton density) scaling factor
• m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
• R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
• R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
• Rex::Real: Exchange rate between the two spin pools in 1/seconds
• R1s::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
• T2s::Real: Transversal relaxation time of the semi-solid pool in seconds

Optional Arguments:

• grad_list=nothing: Tuple that specifies the gradients that are calculated; the default nothing means no gradient, or pass any subset/order of grad_list=(grad_M0(), grad_m0s(), grad_R1f(), grad_R2f(), grad_Rex(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a().
• Ncyc=2: The magnetization is initialized with thermal equilibrium and then performed Ncyc times and only the last cycle is stored. The default value is usually a good approximation for antiperiodic boundary conditions. Increase the number for higher precision at the cost of computation time.
• output=:complexsignal: The default keywords triggers the function to output a complex-valued signal (xf + 1im yf); the keyword output=:realmagnetization triggers an output of the entire (real valued) vector [xf, yf, zf, xs, zs]
• greens=(greens_superlorentzian, dG_o_dT2s_x_T2s_superlorentzian): Tuple of a Greens function G(κ) = G((t-τ)/T2s) and its partial derivative wrt. T2s, multiplied by T2s ∂G((t-τ)/T2s)/∂T2s * T2s. This package supplies the three Greens functions greens=(greens_superlorentzian, dG_o_dT2s_x_T2s_superlorentzian) (default), greens=(greens_lorentzian, dG_o_dT2s_x_T2s_lorentzian), and greens=(greens_gaussian, dG_o_dT2s_x_T2s_gaussian)

Examples

julia> s, g = simulate_gbloch_ide(fill(π/2, 100), fill(5e-4, 100), 4e-3, 0, 1, 1, 0.2, 0.3, 15, 20, 2, 10e-6);

julia> typeof(s)
Vector{ComplexF64} (alias for Array{Complex{Float64}, 1})

julia> size(s)
(100,)

julia> typeof(g)
Nothing

simulate_graham_ode(α, TRF, TR, ω0, B1, M0, m0s, R1f, R2f, Rex, R1s, T2s[; grad_list=nothing, Ncyc=2, output=:complexsignal])

Calculate the signal or magnetization evolution with Graham's spectral model assuming a super-Lorentzian lineshape.

The simulation assumes a sequence of rectangular RF-pulses with varying flip angles α and RF-pulse durations TRF, but a fixed repetition time TR. Further, it assumes balanced gradient moments.

Always returns a tuple (signal, gradients) where signal is a vector (for output=:complexsignal) or matrix (for output=:realmagnetization) scaled by M0, and gradients is a matrix with one column per entry in grad_list, or nothing if no gradients are requested. M0 must be real-valued; for complex-valued M0 (e.g. to account for receive coil phase), simulate with M0=1 and multiply the complex M0 afterward.

Arguments

• α::Vector{Real}: Array of flip angles in radians
• TRF::Vector{Real}: Array of the RF-pulse durations in seconds
• TR::Real: Repetition time in seconds
• ω0::Real: Off-resonance frequency in rad/s
• B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
• M0::Real: Equilibrium magnetization (proton density) scaling factor
• m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
• R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
• R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
• Rex::Real: Exchange rate between the two spin pools in 1/seconds
• R1s::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
• T2s::Real: Transversal relaxation time of the semi-solid pool in seconds

Optional:

• grad_list=nothing: Tuple that specifies the gradients that are calculated; the default nothing means no gradient, or pass any subset/order of grad_list=(grad_M0(), grad_m0s(), grad_R1f(), grad_R2f(), grad_Rex(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a().
• Ncyc=2: The magnetization is initialized with thermal equilibrium and then performed Ncyc times and only the last cycle is stored. The default value is usually a good approximation for antiperiodic boundary conditions. Increase the number for higher precision at the cost of computation time.
• output=:complexsignal: The default keywords triggers the function to output a complex-valued signal (xf + 1im yf); the keyword output=:realmagnetization triggers an output of the entire (real valued) vector [xf, yf, zf, xs, zs]

Examples

julia> s, g = simulate_graham_ode(fill(π/2, 100), fill(5e-4, 100), 4e-3, 0, 1, 1, 0.2, 0.3, 15, 20, 2, 10e-6);

julia> typeof(s)
Vector{ComplexF64} (alias for Array{Complex{Float64}, 1})

julia> size(s)
(100,)

julia> typeof(g)
Nothing

simulate_linearapprox(α, TRF, TR, ω0, B1, M0, m0s, R1f, R2f, Rex, R1s, T2s, R2slT[; grad_list=nothing, rfphase_increment=π, m0=:periodic, output=:complexsignal, grad_moment = [i == 1 ? :spoiler_dual : :balanced for i ∈ eachindex(α)]])

Calculate the signal or magnetization evolution with the linear approximation of the generalized Bloch model assuming a super-Lorentzian lineshape.

The simulation assumes a sequence of rectangular RF-pulses with varying flip angles α and RF-pulse durations TRF, but a fixed repetition time TR. Further, it assumes balanced gradient moments.

Always returns a tuple (signal, gradients) where signal is a vector scaled by M0 and gradients is a matrix with one column per entry in grad_list, or nothing if no gradients are requested. M0 must be real-valued; for complex-valued M0 (e.g. to account for receive coil phase), simulate with M0=1 and multiply the complex M0 afterward.

Arguments

• α::Vector{Real}: Array of flip angles in radians
• TRF::Vector{Real}: Array of the RF-pulse durations in seconds
• TR::Real: Repetition time in seconds
• ω0::Real: Off-resonance frequency in rad/s
• B1::Real: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
• M0::Real: Equilibrium magnetization (proton density) scaling factor
• m0s::Real: Fractional size of the semi-solid pool; should be in range of 0 to 1
• R1f::Real: Longitudinal relaxation rate of the free pool in 1/seconds
• R2f::Real: Transversal relaxation rate of the free pool in 1/seconds
• Rex::Real: Exchange rate between the two spin pools in 1/seconds
• R1s::Real: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
• T2s::Real: Transversal relaxation time of the semi-solid pool in seconds
• R2slT::NTuple{3, Function}: Tuple of three functions: R2sl(TRF, ω1, B1, T2s), dR2sldB1(TRF, ω1, B1, T2s), and R2sldT2s(TRF, ω1, B1, T2s). Can be generated with precompute_R2sl

Optional Arguments:

• grad_list=nothing: Tuple that specifies the gradients that are calculated; the default nothing means no gradient, or pass any subset/order of grad_list=(grad_M0(), grad_m0s(), grad_R1f(), grad_R2f(), grad_Rex(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()); the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a().
• rfphase_increment=π::Real: Increment of the RF phase between consecutive pulses. The default value π, together with $ω0=0$ corresponds to the on-resonance condition.
• m0=:periodic: With the default keyword :periodic, the signal and their derivatives are calculated assuming $m(0) = -m(T)$, where T is the duration of the RF-train. With the keyword :thermal, the magnetization $m(0)$ is initialized with thermal equilibrium [xf, yf, zf, xs, zs] = [0, 0, 1-m0s, 0, m0s], followed by a α[1]/2 - TR/2 prep pulse; and with the keyword :IR, this initialization is followed an inversion pulse of duration TRF[1], (set α[1]=π) and a α[2]/2 - TR/2 prep pulse.
• preppulse=false: if true, a α/2 - TR/2 preparation is applied. In the case of m0=:IR, it is applied after the inversion pulse based on α[2], otherwise it is based on α[1]
• output=:complexsignal: The default keywords triggers the function to output a complex-valued signal (xf + 1im yf); the keyword output=:realmagnetization triggers an output of the entire (real valued) vector [xf, yf, zf, xs, zs, 1]
• grad_moment = [i == 1 ? :spoiler_dual : :balanced for i ∈ eachindex(α)]: Different types of gradient moments of each TR are possible (:balanced, :crusher, :spoiler_dual, :spoiler_prepulse). :balanced simulates a TR with all gradient moments nulled. :crusher assumes equivalent (non-zero) gradient moments before and simulates the refocussing path of the extended phase graph. :spoiler_prepulse nulls all transverse magnetization before the RF pulse, emulating an idealized FLASH. :spoiler_dual nulls all transverse magnetization before and after the RF pulse.

Examples

julia> R2slT = precompute_R2sl();

julia> s, g = simulate_linearapprox(range(0, π/2, 100), fill(5e-4, 100), 4e-3, 0, 1, 1, 0.1, 1, 15, 30, 6.5, 10e-6, R2slT);

julia> typeof(s)
Vector{ComplexF64} (alias for Array{Complex{Float64}, 1})

julia> size(s)
(100,)

julia> typeof(g)
Nothing

magnetization = steady_state_magnetization(cycle_ops, propagators)

Solve for the periodic steady-state magnetization and propagate it through all time steps. Returns magnetization[t, g] as an 11-component SVector for each time step t and gradient parameter g.

cycle_ops = steady_state_operator(propagators)

Compute the cycle operator Q[g] = A₀(A) - A for each gradient parameter, where A = E[1] · E[end] · ... · E[2] is the full-cycle propagator. The steady-state magnetization satisfies Q · m = c, where c encodes the thermal equilibrium source and A₀ zeroes the affine (thermal equilibrium) row.