API

CRB, grad_ω1, grad_TRF = CRB_gradient_OCT(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, weights; isInversionPulse = [true; falses(length(ω1)-1)])

Calculate the Cramer-Rao bound of a pulse sequence along with the derivatives wrt. ω1 and TRF.

Arguments

• ω1::Vector{<:Number}: Control vector of length = Npulses
• TRF::Vector{<:Number}: Control vector of length = Npulses
• TR::Number: Repetition time in seconds
• ω0::Number: Off-resonance frequency in rad/s
• B1::Number: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
• m0s::Number: Fractional size of the semi-solid pool; should be in range of 0 to 1
• R1f::Number: Longitudinal relaxation rate of the free pool in 1/seconds
• R2f::Number: Transversal relaxation rate of the free pool in 1/seconds
• Rx::Number: Exchange rate between the two spin pools in 1/seconds
• R1f::Number: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
• T2s::Number: 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::Vector{<:grad_param}: Vector that specifies the gradients that are calculated; the vector elements can either be any subset/order of grad_list=[grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]; the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a()
• weights::transpose(Vector{<:Number}): Row vector of weights applied to the Cramer-Rao bounds (CRB) of the individual parameters. The first entry always refers to the CRB of M0, followed by the values defined in grad_list in the order defined therein. Hence, the vector weights has to have one more entry than grad_list

Optional Keyword Arguments:

• isInversionPulse::Vector{Bool}: Indicates all inversion pulses; must have the same length as α; the default = [true; falses(length(ω1)-1)] indicates that the first pulse is an inversion pulse and all others are not

Examples

julia> CRB, grad_ω1, grad_TRF = MRIgeneralizedBloch.CRB_gradient_OCT(rand(100) .* π, rand(100) .* 400e-6 .+ 100e-6, 3.5e-3, 0, 1, 0.15, 0.5, 15, 30, 2, 10e-6, precompute_R2sl(), [grad_m0s(), grad_R1f()], transpose([0, 1, 1]); isInversionPulse = [true, falses(99)...])
(4.326345257791475e20, [0.0, 3.8149138162765024e20, -1.4073464193072628e20, 1.196651310804305e20, -1.9989638164742383e20, 3.6493209787269874e20, -1.5536324770044346e20, 1.7284798914645443e20, -2.0772522228321904e20, 2.7044669618035792e20  …  -6.3928347250126905e19, 1.5671768655962505e20, -1.1080723702078613e20, 5.7454290837437145e19, -7.845659845122268e19, 1.3447202941410777e20, -6.377328799158878e19, 5.1500654841007415e19, -3.97683716404465e19, 1.2714666624472279e20], [0.0, 2.9969441368202493e24, -1.4650231537214497e24, 1.8226585318981564e24, -4.7462030797597407e23, 1.97549455130354e24, -1.5270500210926142e24, 9.099568864141483e23, -1.4811032902868818e24, 1.2220830933233262e24  …  -5.258326682397155e23, 6.918585992930617e23, -2.8293355719514886e23, 4.015546075884712e23, -4.1767064624316516e23, 5.061776753372734e23, -4.1049901738345103e23, 9.311794792265391e23, -4.131984395429261e23, 7.610379029805631e23])

c.f. Optimal Control

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

Calculate RF power penalty and over-write the gradients in place.

Arguments

• grad_ω1::Vector{<:Number}: Gradient of control, which will be over-written in place
• grad_TRF::Vector{<:Number}: Gradient of control, which will be over-written in place
• ω1::Vector{<:Number}: Control vector
• TRF::Vector{<:Number}: Control vector

Optional Keyword Arguments:

• λ::Number: regularization parameter
• Pmax::Number: 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::Number: Repetition time of the pulse sequence

Examples

julia> ω1 = 4000π * rand(100);

julia> TRF = 500e-6 * rand(100);

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

julia> grad_TRF = similar(ω1);

julia> F = MRIgeneralizedBloch.RF_power!(grad_ω1, grad_TRF, ω1, TRF; λ=1e3, Pmax=3e5)
9.54432950518758e15

F = TRF_TV!(grad_TRF, ω1, TRF; λ = 1)

Calculate the total variation penalty of TRF and over-write grad_TRF in place.

Arguments

• grad_TRF::Vector{<:Number}: Gradient of control, which will be over-written in place
• ω1::Vector{<:Number}: Control vector
• TRF::Vector{<:Number}: Control vector

Optional Keyword Arguments:

• λ::Number: regularization parameter

Examples

julia> ω1 = 4000π * rand(100);

julia> TRF = 500e-6 * rand(100);

julia> grad_TRF = similar(ω1);

julia> F = MRIgeneralizedBloch.TRF_TV!(grad_TRF, ω1, TRF; λ = 1e-3)
1.826335488333797e-5

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{<:Number}: 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{<:Number}: 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) <: Number ? 0.0 : zeros(5n + 5) for n gradients, and is then updated by the delay differential equation solvers
• p::NTuple{9,10, or 11, Any}: (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g), with -ω1::Number: Rabi frequency in rad/s (rotation about the y-axis) -B1::Number: B1 scaling normalized so that B1=1 corresponds to a perfectly calibrated RF field -ω0::Number: Larmor or off-resonance frequency in rad/s -m0s::Number: Fractional semi-solid spin pool size in the range of 0 to 1 -R1f::Number: Longitudinal spin relaxation rate of the free pool in 1/seconds -R2f::Number: Transversal spin relaxation rate of the free pool in 1/seconds -Rx::Number: Exchange rate between the two pools in 1/seconds -R1s::Number: Longitudinal spin relaxation rate of the semi-solid pool in 1/seconds -T2s::Number: Transversal spin relaxation time of the semi-solid pool in seconds -g::Function: Green's function of the form G(κ) = G((t-τ)/T2s) or (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, zs_idx, g) with
  • zs_idx::Integer: Index to be used history function to be used in the Green's function; Default is 4 (zs), and for derivatives 9, 14, ... are used
  or (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g, dG_o_dT2s_x_T2s, grad_list) with
  • 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; any subset of [grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), 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::Number: 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 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> Rx = 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) <: Number ? 0.0 : zeros(5);

julia> sol = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G)))
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 9-element Vector{Float64}:
 0.0
 1.375006182301112e-7
 1.512506800531223e-6
 8.042561696923577e-6
 2.107848894861101e-5
 3.911414415070652e-5
 6.26879093553081e-5
 9.147705752659822e-5
 0.0001
u: 9-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.8, 0.2, 1.0]
 [0.0017278806030763402, 0.0, 0.7999981340131751, 0.19999953350448, 1.0]
 [0.019004717382235078, 0.0, 0.7997742277135814, 0.19994357804868362, 1.0]
 [0.10079111348917136, 0.0, 0.7936248122939504, 0.19842287240365722, 1.0]
 [0.26002578672576243, 0.0, 0.7565529666157937, 0.18981913039644657, 1.0]
 [0.4610419882566734, 0.0, 0.6537242214798688, 0.16937688382096108, 1.0]
 [0.6661738538876186, 0.0, 0.44261236945975563, 0.1358931514238721, 1.0]
 [0.7923116826717905, 0.0, 0.10713144280454787, 0.09390268562369869, 1.0]
 [0.7994211188440815, 0.0, 0.0004403374355099447, 0.08214809683848684, 1.0]

julia> using Plots

julia> plot(sol, labels=["xf" "yf" "zf" "zs" "1"], xlabel="t (s)", ylabel="m(t)");




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> mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? 0.0 : zeros(5 + 5*length(grad_list));

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




julia> plot(sol);

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{<:Number}: 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{<:Number}: 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, Rx, R1s, T2s, g) for simulating a coupled spin system; with -ω1::Number: Rabi frequency in rad/s (rotation about the y-axis) -B1::Number: B1 scaling normalized so that B1=1 corresponds to a perfectly calibrated RF field -ω0::Number: Larmor or off-resonance frequency in rad/s -R1f::Number: Longitudinal spin relaxation rate of the free pool in 1/seconds -R2f::Number: Transversal spin relaxation rate of the free pool in 1/seconds -R1s::Number: Longitudinal spin relaxation rate of the semi-solid in 1/seconds -Rx::Number: Exchange rate between the two pools in 1/seconds -T2s::Number: Transversal spin relaxation time in seconds -g::Function: Green's function of the form G(κ) = G((t-τ)/T2s)
• t::Number: 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())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 3-element Vector{Float64}:
 0.0
 7.475414666720001e-5
 0.0001
u: 3-element Vector{Vector{Float64}}:
 [1.0]
 [0.631392823181197]
 [0.48953654496619187]

julia> using Plots

julia> plot(sol, labels=["zs"], xlabel="t (s)", ylabel="m(t)");

x = bound_ω1_TRF!(ω1, TRF; ω1_min = 0, ω1_max = 2e3π, TRF_min = 100e-6, TRF_max = 500e-6)

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

Arguments

• ω1::Vector{<:Number}: Control vector of length Npulses
• TRF::Vector{<:Number}: Control vector of length Npulses

Optional Keyword Arguments (see above for defaults):

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

Examples

julia> ω1 = 4000π * rand(100);

julia> TRF = 500e-6 * rand(100);

julia> x = MRIgeneralizedBloch.bound_ω1_TRF!(ω1, TRF)
200-element Vector{Float64}:
  Inf
   0.16123919298957837
  Inf
  Inf
  Inf
  Inf
   0.17719147195925186
   0.33617937553604216
  -1.1288911659500056
  Inf
   ⋮
 -Inf
  -0.5475311433614531
  -0.6572991695812251
 -Inf
  -0.5487931795341404
  -2.097250369053942
 -Inf
  -1.1536520473794025
   0.07161170293075586

calculatesignal_gbloch_ide(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s[; grad_list, 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.

Arguments

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

Optional:

• grad_list=[]: Vector that specifies the gradients that are calculated; the vector can either be empty [] for no gradient, or contain any subset/order of grad_list=[grad_m0s(), grad_R1s(), grad_R2f(), grad_Rx(), 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> calculatesignal_gbloch_ide(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6)
100×1 Matrix{ComplexF64}:
  -0.007966316445016917 + 0.0im
  0.0012590590419314762 - 0.0im
  -0.006088855588122765 + 0.0im
  0.0024187389404174905 - 0.0im
  -0.004361339395232529 + 0.0im
  0.0034891358210049914 - 0.0im
 -0.0027633710614652324 + 0.0im
   0.004483217941394411 - 0.0im
 -0.0012812573517353436 + 0.0im
   0.005408854034677167 - 0.0im
                        ⋮
   0.017760808273048677 - 0.0im
   0.017576118974646303 + 0.0im
   0.017813950945910262 - 0.0im
   0.017643856335506514 + 0.0im
   0.017863575855931582 - 0.0im
   0.017706926033852457 + 0.0im
   0.017909914934081783 - 0.0im
    0.01776565037329605 + 0.0im
    0.01795318489371729 - 0.0im

julia> calculatesignal_gbloch_ide(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6; grad_list=[grad_R1f(), grad_T2s()], output=:realmagnetization)
100×15 transpose(::Matrix{Float64}) with eltype Float64:
 -0.00796627   0.0   0.000637773  …   0.0   -10.8757  -335.26   0.0
  0.00125903  -0.0  -0.00700671      -0.0   125.882   -326.977  0.0
 -0.00608882   0.0   0.00185086       0.0   -30.4187  -317.56   0.0
  0.00241873  -0.0  -0.00520622      -0.0    96.1776  -309.906  0.0
 -0.00436133   0.0   0.00296471       0.0   -47.5803  -302.948  0.0
  0.003489    -0.0  -0.00354518   …  -0.0    69.5148  -298.697  0.0
 -0.00276366   0.0   0.00399588       0.0   -62.8453  -294.886  0.0
  0.00448273  -0.0  -0.00200673      -0.0    45.3179  -292.783  0.0
 -0.00128187   0.0   0.00495478       0.0   -76.6573  -290.321  0.0
  0.00540814  -0.0  -0.000578836     -0.0    23.1756  -289.245  0.0
  ⋮                               ⋱
  0.0177563   -0.0   0.0175372       -0.0  -290.779   -349.855  0.0
  0.0175716    0.0   0.0177845        0.0  -295.347   -350.002  0.0
  0.0178094   -0.0   0.0176073       -0.0  -292.44    -350.163  0.0
  0.0176393    0.0   0.0178359        0.0  -296.668   -350.3    0.0
  0.017859    -0.0   0.0176727    …  -0.0  -294.001   -350.451  0.0
  0.0177024    0.0   0.0178838        0.0  -297.914   -350.579  0.0
  0.0179053   -0.0   0.0177335       -0.0  -295.467   -350.72   0.0
  0.0177611    0.0   0.0179286        0.0  -299.09    -350.84   0.0
  0.0179486   -0.0   0.0177902       -0.0  -296.845   -350.972  0.0

calculatesignal_graham_ode(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s[; grad_list, 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.

Arguments

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

Optional:

• grad_list=[]: Vector that specifies the gradients that are calculated; the vector can either be empty [] for no gradient, or contain any subset/order of grad_list=[grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1(), 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> calculatesignal_graham_ode(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6)
100×1 Matrix{ComplexF64}:
  -0.008073451191813547 + 0.0im
    0.00126864329048296 - 0.0im
  -0.006178694437239949 + 0.0im
   0.002435865817867096 - 0.0im
  -0.004437476277201384 + 0.0im
  0.0035164646498403165 - 0.0im
 -0.0028315542520438355 + 0.0im
   0.004523902250379148 - 0.0im
 -0.0013422299483389299 + 0.0im
   0.005454562034842184 - 0.0im
                        ⋮
   0.018148222040953982 - 0.0im
    0.01795769661414948 + 0.0im
    0.01820486087666198 - 0.0im
   0.018029363102618986 + 0.0im
   0.018257820142562945 - 0.0im
    0.01809616891246802 + 0.0im
   0.018307337729977703 - 0.0im
   0.018158444501177453 + 0.0im
   0.018353636229654635 - 0.0im

julia> calculatesignal_graham_ode(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6; grad_list=[grad_R1f(), grad_T2s()], output=:realmagnetization)
100×15 transpose(::Matrix{Float64}) with eltype Float64:
 -0.0080756    0.0   0.000643162  …   0.0   -10.4986  -323.634  0.0
  0.00126867  -0.0  -0.00710692      -0.0   123.078   -316.358  0.0
 -0.00618067   0.0   0.00186482       0.0   -29.4458  -307.862  0.0
  0.00243634  -0.0  -0.0052899       -0.0    94.1692  -300.821  0.0
 -0.00443746   0.0   0.00298646       0.0   -46.1422  -294.116  0.0
  0.00351386  -0.0  -0.00361718   …  -0.0    68.262   -289.421  0.0
 -0.00282882   0.0   0.00402793       0.0   -61.0236  -285.242  0.0
  0.00451808  -0.0  -0.00206741      -0.0    44.7888  -282.867  0.0
 -0.00133608   0.0   0.00499316       0.0   -74.4075  -280.79   0.0
  0.00544896  -0.0  -0.00062663      -0.0    23.2536  -280.148  0.0
  ⋮                               ⋱
  0.0181755   -0.0   0.0179431       -0.0  -285.071   -337.48   0.0
  0.0179911    0.0   0.0181963        0.0  -289.591   -337.634  0.0
  0.0182343   -0.0   0.018019        -0.0  -286.777   -337.804  0.0
  0.0180645    0.0   0.0182531        0.0  -290.964   -337.948  0.0
  0.0182893   -0.0   0.0180898    …  -0.0  -288.383   -338.106  0.0
  0.0181329    0.0   0.0183062        0.0  -292.261   -338.241  0.0
  0.0183407   -0.0   0.0181559       -0.0  -289.895   -338.39   0.0
  0.0181968    0.0   0.0183559        0.0  -293.488   -338.516  0.0
  0.0183888   -0.0   0.0182175       -0.0  -291.318   -338.656  0.0

calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT[; grad_list=[undef], rfphase_increment=[π], m0=:periodic, output=:complexsignal, isInversionPulse = [true; falses(length(α)-1)]])

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.

Arguments

• α::Vector{<:Number}: Array of flip angles in radians
• TRF::Vector{<:Number}: Array of the RF-pulse durations in seconds
• TR::Number: Repetition time in seconds
• ω0::Number: Off-resonance frequency in rad/s
• B1::Number: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
• m0s::Number: Fractional size of the semi-solid pool; should be in range of 0 to 1
• R1f::Number: Longitudinal relaxation rate of the free pool in 1/seconds
• R2f::Number: Transversal relaxation rate of the free pool in 1/seconds
• Rx::Number: Exchange rate between the two spin pools in 1/seconds
• R1f::Number: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
• T2s::Number: 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:

• grad_list=[undef]: Vector that specifies the gradients that are calculated; the vector elements can either be undef for no gradient, or any subset/order of grad_list=[grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]; the derivative wrt. to apparent R1a = R1f = R1s can be calculated with grad_R1a()
• rfphase_increment=[π]::Vector{<:Number}: Increment of the RF phase between consecutive pulses. The default value π, together with $ω0=0$ corresponds to the on-resonance condition. When more than one value is supplied, their resulting signal is stored along the second dimension of the output array.
• 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.
• 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]
• isInversionPulse::Vector{Bool}: Indicates all inversion pulses; must have the same length as α; the default = [true; falses(length(α)-1)] indicates that the first pulse is an inversion pulse and all others are not

Examples

julia> R2slT = precompute_R2sl();


julia> calculatesignal_linearapprox(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.1, 1, 15, 30, 6.5, 10e-6, R2slT)
100×1×1 Array{ComplexF64, 3}:
[:, :, 1] =
  -0.029305987774458163 - 2.6090233074873602e-18im
   0.004329424678273618 - 9.228587093606358e-19im
  -0.022761409373843695 + 1.3684406715212817e-18im
   0.008280330224850314 - 3.913888063486213e-18im
  -0.016751305727868086 + 4.640955399485587e-18im
   0.011921649143708494 - 6.302663734240699e-18im
   -0.01119057989041819 + 7.310583447081545e-18im
   0.015305608440952356 - 8.175490153983926e-18im
 -0.0060267918180664974 + 9.464626759459644e-18im
   0.018463027499697613 - 9.608537136498027e-18im
                        ⋮
   0.061531473509660414 + 1.4962743251121995e-18im
    0.06081400768357244 + 5.6874647068595666e-18im
     0.0617257515931871 + 1.6576028430814591e-18im
   0.061064216222629225 + 5.55793788363693e-18im
     0.0619074570531536 + 1.8084438857930338e-18im
    0.06129750535994419 + 5.436837755621582e-18im
   0.062077399758002534 + 1.9493976022865502e-18im
   0.061515021919442275 + 5.323692528202896e-18im
    0.06223633770306246 + 2.0810405785134972e-18im

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.095033670965

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),
    Rx   = (   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{<:Number}: Array of flip angles in radians
• TRF::Vector{<:Number}: Array of the RF-pulse durations in seconds
• TR::Number: Repetition time in seconds
• ω0::Number: Off-resonance frequency in rad/s

Optional Keyword Arguments:

• reM0::Union{Number, Tuple{Number, Number, Number}}: Real part of M0; either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• imM0::Union{Number, Tuple{Number, Number, Number}}: Imaginary part of M0; either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• m0s::Union{Number, Tuple{Number, Number, Number}}: Fractional size of the semi-solid pool (should be in range of 0 to 1); either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• R1f::Union{Number, Tuple{Number, Number, Number}}: Longitudinal relaxation rate of the free pool in 1/s; only used in combination with fit_apparentR1=false; either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• R2f::Union{Number, Tuple{Number, Number, Number}}: Transversal relaxation rate of the free pool in 1/s; either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• Rx::Union{Number, Tuple{Number, Number, Number}}: Exchange rate between the two spin pools in 1/s; either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• R1s::Union{Number, Tuple{Number, Number, Number}}: Longitudinal relaxation rate of the semi-solid pool in 1/s; only used in combination with fit_apparentR1=false; either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• T2s::Union{Number, Tuple{Number, Number, Number}}: Transversal relaxation time of the semi-solid pool in s; either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• ω0::Union{Number, Tuple{Number, Number, Number}}: Off-resonance frequency in rad/s; either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• B1::Union{Number, Tuple{Number, Number, Number}}: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field; either fixed value as a Number or fit limits thereof as a Tuple with the elements (min, start, max)
• R1a::Union{Number, Tuple{Number, Number, Number}}: Apparent longitudinal relaxation rate in 1/s; only used in combination with fit_apparentR1=true; either fixed value as a Number 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_ω1_TRF(x; ω1_min = 0, ω1_max = 2e3π, TRF_min = 100e-6, TRF_max = 500e-6)

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{<:Number}: Control vector of length = 2Npulses with values in the range [-Inf, Inf]

Optional Keyword Arguments:

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

Examples

julia> x = 1000 * randn(2 * 100);

julia> ω1, TRF = MRIgeneralizedBloch.get_bounded_ω1_TRF(x)
([6283.185307179586, 0.0, 6283.185307179586, 0.0, 6283.185307179586, 0.0, 6283.185307179586, 0.0, 6283.185307179586, 6283.185307179586  …  6283.185307179586, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6283.185307179586, 0.0, 0.0], [0.0005, 0.0005, 0.0005, 0.0001, 0.0005, 0.0001, 0.0001, 0.0005, 0.0005, 0.0005  …  0.0001, 0.0001, 0.0005, 0.0005, 0.0001, 0.0005, 0.0001, 0.0005, 0.0001, 0.0001])

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.1471246868094442

hamiltonian_linear(ω1, B1, ω0, T, m0s, R1f, R2f, Rx, 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::Number: Rabi frequency in rad/s (rotation about the y-axis)
• B1::Number: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
• ω0::Number: Larmor (or off-resonance) frequency in rad/s (rotation about the z-axis)
• T::Number: Time in seconds; this can, e.g., be the RF-pulse duration, or the time of free precession with ω1=0
• m0s::Number: Fractional size of the semi-solid pool; should be in range of 0 to 1
• R1f::Number: Longitudinal relaxation rate of the free pool in 1/seconds
• R2f::Number: Transversal relaxation rate of the free pool in 1/seconds
• Rx::Number: Exchange rate between the two spin pools in 1/seconds
• R1s::Number: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
• R2s::Number: 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:

• dR2sdT2s::Number: 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::Number: 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_Rx(), 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> Rx = 30;

julia> R1s = 6.5;

julia> R2s = 1e5;

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

julia> (xf, yf, zf, xs, zs, _) = exp(hamiltonian_linear(ω1, B1, ω0, T, m0s, R1f, R2f, Rx, R1s, R2s)) * m0
6-element StaticArrays.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.1471246868094442

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


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

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.5, 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::Number: lower bound of the RF-pulse duration range in seconds
• TRF_max::Number: upper bound of the RF-pulse duration range in seconds
• T2s_min::Number: lower bound of the T2s range in seconds
• T2s_max::Number: upper bound of the T2s range in seconds
• ω1_max::Number: upper bound of the Rabi frequency ω1, the default is the frequency of a 500μs long π-pulse
• B1_max::Number: 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.3, greens=greens_gaussian);

F = second_order_α!(grad_ω1, grad_TRF, ω1, TRF; λ = 1)

Calculate second order penalty of variations of the flip angle α and over-write the gradients in place.

Arguments

• grad_ω1::Vector{<:Number}: Gradient of control, which will be over-written in place
• grad_TRF::Vector{<:Number}: Gradient of control, which will be over-written in place
• ω1::Vector{<:Number}: Control vector
• TRF::Vector{<:Number}: Control vector

Optional Keyword Arguments:

• λ::Number: regularization parameter

Examples

julia> ω1 = 4000π * rand(100);

julia> TRF = 500e-6 * rand(100);

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

julia> grad_TRF = similar(ω1);

julia> F = MRIgeneralizedBloch.second_order_α!(grad_ω1, grad_TRF, ω1, TRF; λ = 1e-3)
0.3394498406831758