API

In the following, you find the documentation of all exported functions of the MRIgeneralizedBloch.jl package:

MRIgeneralizedBloch.apply_hamiltonian_gbloch!
MRIgeneralizedBloch.apply_hamiltonian_sled!
MRIgeneralizedBloch.calculatesignal_gbloch_ide
MRIgeneralizedBloch.calculatesignal_graham_ode
MRIgeneralizedBloch.calculatesignal_linearapprox
MRIgeneralizedBloch.dG_o_dT2s_x_T2s_gaussian
MRIgeneralizedBloch.dG_o_dT2s_x_T2s_lorentzian
MRIgeneralizedBloch.dG_o_dT2s_x_T2s_superlorentzian
MRIgeneralizedBloch.greens_gaussian
MRIgeneralizedBloch.greens_lorentzian
MRIgeneralizedBloch.greens_superlorentzian
MRIgeneralizedBloch.hamiltonian_linear
MRIgeneralizedBloch.interpolate_greens_function
MRIgeneralizedBloch.precompute_R2sl

MRIgeneralizedBloch.apply_hamiltonian_gbloch! — Method

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.

Arguemnts

∂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 fuction; 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, R1, R2f, T2s, Rx, 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 -R1::Number: Apparent longitudinal spin relaxation rate of both pools in 1/seconds -R2f::Number: Trasversal spin relaxation rate of the free pool in 1/seconds -T2s::Number: Trasversal spin relaxation time of the semi-solid pool in seconds -Rx::Number: Exchange rate between the two pools in 1/seconds -g::Function: Green's function of the form G(κ) = G((t-τ)/T2s) or (ω1, B1, ω0, m0s, R1, R2f, T2s, Rx, 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, R1, R2f, T2s, Rx, 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 calucualted; any subset of [grad_m0s(), grad_R1(), grad_R2f(), grad_Rx(), grad_T2s(), grad_ω0(), grad_B1()]; length of the vector must be n (cf. arguments m and ∂m∂t)
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 inteded 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.1;

julia> R1 = 1;

julia> R2f = 15;

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, R1, R2f, T2s, Rx, G)))
retcode: Success
Interpolation: automatic order switching interpolation
t: 9-element Vector{Float64}:
 0.0
 1.220281289257312e-7
 1.342309418183043e-6
 7.538223396809993e-6
 2.0264275992449307e-5
 3.806391790996879e-5
 6.131350253044042e-5
 8.985713095394202e-5
 0.0001
u: 9-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.9, 0.1, 1.0]
 [0.0017251293948764095, 0.0, 0.8999983466235149, 0.0999998162918461, 1.0]
 [0.018974855260649, 0.0, 0.8997999500976795, 0.0999777787537773, 1.0]
 [0.10631425327309292, 0.0, 0.8936981899494069, 0.0993062488692229, 1.0]
 [0.28162336805181576, 0.0, 0.8547938784075204, 0.09527129619977613, 1.0]
 [0.5064779297334128, 0.0, 0.7438963250738748, 0.08537276746991054, 1.0]
 [0.7385374607670901, 0.0, 0.5140010917423798, 0.06896091618108181, 1.0]
 [0.888016931522836, 0.0, 0.14315426117917585, 0.04808883373953286, 1.0]
 [0.899347249517465, 0.0, 0.00048588799488130663, 0.04106861705554077, 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, R1, R2f, T2s, Rx, G, dG_o_dT2s_x_T2s, grad_list)));


julia> plot(sol);

source

MRIgeneralizedBloch.apply_hamiltonian_sled! — Method

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

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

Arguemnts

∂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 9, Any}: (ω1, B1, ω0, R1, T2s, g) for a simulating an isolated semi-solid pool or (ω1, B1, ω0, m0s, R1, R2f, T2s, Rx, 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 -R1::Number: Longitudinal spin relaxation rate in 1/seconds -R2f::Number: Trasversal spin relaxation rate of the free pool in 1/seconds -T2s::Number: Trasversal spin relaxation time in seconds -Rx::Number: Exchange rate between the two pools in 1/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> R1 = 1;

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, R1, T2s, G)), Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 3-element Vector{Float64}:
 0.0
 7.475658194333419e-5
 0.0001
u: 3-element Vector{Vector{Float64}}:
 [1.0]
 [0.6313685535188782]
 [0.48951919836592006]

julia> using Plots

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

source

MRIgeneralizedBloch.calculatesignal_gbloch_ide — Method

calculatesignal_gbloch_ide(α, TRF, TR, ω0, B1, m0s, R1, R2f, Rx, T2s[; grad_list, Ncyc=2, output=:complexsignal])

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

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

Arguemnts

α::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
R1::Number: Apparent longitudinal relaxation rate of the free and semi-solid 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
T2s::Number: Transversal relaxationt time of the semi-solid pool in seconds

Optional:

grad_list=[]: Vector to indicate which gradients should be calculated; the vector can either be empty [] for no gradient, or contain any subset/order of grad_list=[grad_m0s(), grad_R1(), grad_R2f(), grad_Rx(), grad_T2s(), grad_ω0(), grad_B1()]
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 defaul 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.1, 1, 15, 30, 10e-6)
100×1 Matrix{ComplexF64}:
 -0.024657762441422027 + 0.0im
 0.0037348678313655435 - 0.0im
 -0.019057736703007047 + 0.0im
  0.007146413346758778 - 0.0im
 -0.013913423956595785 + 0.0im
  0.010291046549792265 - 0.0im
 -0.009153866378612775 + 0.0im
  0.013213045210360654 - 0.0im
 -0.004734258510785772 + 0.0im
   0.01593906991792929 - 0.0im
                       ⋮
   0.05321851165156517 - 0.0im
   0.05261662009092025 + 0.0im
  0.053387874462524944 - 0.0im
  0.052832959843114265 + 0.0im
   0.05354631440847341 - 0.0im
  0.053034722397620235 + 0.0im
   0.05369453263373485 - 0.0im
   0.05322289238188484 + 0.0im
   0.05383318553569216 - 0.0im

julia> calculatesignal_gbloch_ide(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.1, 1, 15, 30, 10e-6; grad_list=[grad_R1(), grad_T2s()], output=:realmagnetization)
100×15 transpose(::Matrix{Float64}) with eltype Float64:
 -0.0246575    0.0   0.00191834  …   0.0    -8.09518  -146.358   0.0
  0.00373479  -0.0  -0.0217727      -0.0    83.521    -135.807   0.0
 -0.0190576    0.0   0.00549726      0.0   -21.4616   -119.477   0.0
  0.00714631  -0.0  -0.0164074      -0.0    64.4139   -106.11    0.0
 -0.0139133    0.0   0.0087715       0.0   -30.9968    -92.3956  0.0
  0.0102909   -0.0  -0.0114604   …  -0.0    49.4596    -83.7202  0.0
 -0.00915379   0.0   0.0118022       0.0   -37.5344    -75.2901  0.0
  0.0132129   -0.0  -0.00687472     -0.0    37.423     -71.2706  0.0
 -0.00473424   0.0   0.0146242       0.0   -42.288     -66.9254  0.0
  0.0159389   -0.0  -0.00261256     -0.0    27.1804    -66.0774  0.0
  ⋮                              ⋱
  0.0532167   -0.0   0.0525262      -0.0  -206.85     -155.403   0.0
  0.0526148    0.0   0.0533279       0.0  -210.517    -155.62    0.0
  0.053386    -0.0   0.0527502      -0.0  -209.167    -155.899   0.0
  0.0528312    0.0   0.0534917       0.0  -212.584    -156.103   0.0
  0.0535445   -0.0   0.052959    …  -0.0  -211.372    -156.364   0.0
  0.0530329    0.0   0.0536449       0.0  -214.557    -156.555   0.0
  0.0536927   -0.0   0.0531539      -0.0  -213.471    -156.799   0.0
  0.053221     0.0   0.0537883       0.0  -216.439    -156.979   0.0
  0.0538313   -0.0   0.0533355      -0.0  -215.468    -157.207   0.0

source

MRIgeneralizedBloch.calculatesignal_graham_ode — Method

calculatesignal_graham_ode(α, TRF, TR, ω0, B1, m0s, R1, R2f, Rx, T2s[; grad_list, Ncyc=2, output=:complexsignal])

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

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

Arguemnts

α::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
R1::Number: Apparent longitudinal relaxation rate of the free and semi-solid 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
T2s::Number: Transversal relaxationt time of the semi-solid pool in seconds

Optional:

grad_list=[]: Vector to indicate which gradients should be calculated; the vector can either be empty [] for no gradient, or contain any subset/order of grad_list=[grad_m0s(), grad_R1(), grad_R2f(), grad_Rx(), grad_T2s(), grad_ω0(), grad_B1()]
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 defaul 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.1, 1, 15, 30, 10e-6)
100×1 Matrix{ComplexF64}:
  -0.02507016283364451 + 0.0im
 0.0037430687099323156 - 0.0im
  -0.01943221197276058 + 0.0im
  0.007158922245383863 - 0.0im
 -0.014255325151363943 + 0.0im
  0.010307593338620509 - 0.0im
 -0.009486903618758917 + 0.0im
  0.013252701136887458 - 0.0im
 -0.005050073780485204 + 0.0im
  0.015978096974037494 - 0.0im
                       ⋮
   0.05452107922170709 - 0.0im
   0.05378824747234005 + 0.0im
  0.054667746773539264 - 0.0im
   0.05399168955011588 + 0.0im
   0.05480524072601236 - 0.0im
   0.05418157100953124 + 0.0im
   0.05493412868432558 - 0.0im
   0.05435879806126502 + 0.0im
  0.055054944205797034 - 0.0im

julia> calculatesignal_graham_ode(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.1, 1, 15, 30, 10e-6; grad_list=[grad_R1(), grad_T2s()], output=:realmagnetization)
100×15 transpose(::Matrix{Float64}) with eltype Float64:
 -0.0249614    0.0   0.00192295  …   0.0    -7.79895  -141.002   0.0
  0.00374263  -0.0  -0.0220613      -0.0    81.4297   -131.309   0.0
 -0.0193318    0.0   0.00550778      0.0   -20.7551   -116.146   0.0
  0.00715882  -0.0  -0.0166685      -0.0    62.8368   -103.541   0.0
 -0.014162     0.0   0.00878636      0.0   -30.0836    -90.4788  0.0
  0.010306    -0.0  -0.0117047   …  -0.0    48.2494    -81.8859  0.0
 -0.00938721   0.0   0.0118269       0.0   -36.5357    -73.5855  0.0
  0.0132393   -0.0  -0.00710492     -0.0    36.5226    -69.4018  0.0
 -0.0049538    0.0   0.0146514       0.0   -41.1928    -65.1257  0.0
  0.0159661   -0.0  -0.00282248     -0.0    26.5349    -64.1602  0.0
  ⋮                              ⋱
  0.0541053   -0.0   0.053379       -0.0  -202.073    -149.963   0.0
  0.0535064    0.0   0.0541954       0.0  -205.711    -150.191   0.0
  0.054291    -0.0   0.0536194      -0.0  -204.431    -150.476   0.0
  0.0537389    0.0   0.054375        0.0  -207.825    -150.691   0.0
  0.054465    -0.0   0.0538441   …  -0.0  -206.681    -150.958   0.0
  0.0539561    0.0   0.0545433       0.0  -209.847    -151.16    0.0
  0.0546281   -0.0   0.0540539      -0.0  -208.826    -151.411   0.0
  0.054159     0.0   0.054701        0.0  -211.78     -151.601   0.0
  0.0547809   -0.0   0.05425        -0.0  -210.87     -151.835   0.0

source

MRIgeneralizedBloch.calculatesignal_linearapprox — Method

calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1, R2f, Rx, T2s, R2slT[; grad_list=[undef], rfphase_increment=[π], m0=:antiperiodic, output=:complexsignal])

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

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

Arguemnts

α::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
R1::Number: Apparent longitudinal relaxation rate of the free and semi-solid 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
T2s::Number: Transversal relaxationt 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 to indicate which gradients should be calculated; the vector elements can either be undef for no gradient, or any subset/order of grad_list=[grad_m0s(), grad_R1(), grad_R2f(), grad_Rx(), grad_T2s(), grad_ω0(), grad_B1()]
rfphase_increment=[π]::Vector{<:Number}: Increment of the RF phase between consequtive 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=:antiperiodic: With the default keyword :antiperiodic, the signal and their derivatives are calcualted 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 initalization is followed an inversion pulse of duration TRF[1], (set α[1]=π) and a α[2]/2 - TR/2 prep pulse.
output=:complexsignal: The defaul 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> R2slT = precompute_R2sl(4e-4, 6e-4, 5e-6, 15e-6, 0, π, 0.9, 1.1);


julia> calculatesignal_linearapprox(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.1, 1, 15, 30, 10e-6, R2slT)
100×1×1 Array{ComplexF64, 3}:
[:, :, 1] =
  -0.02534278046134154 - 6.748418787451518e-18im
  0.003747524854957436 + 3.4325653058017187e-18im
   -0.0196833907083505 - 2.8004553842504324e-18im
  0.007162011566093652 + 3.672271808540125e-19im
 -0.014489517679344808 + 4.801739729077009e-19im
  0.010305804983291571 - 2.1233287976034503e-18im
 -0.009686634863202483 + 3.1882740494855564e-18im
   0.01322536234967787 - 4.1197910372278614e-18im
 -0.005228086823125368 + 5.405192374297665e-18im
  0.015948493046190188 - 5.693388849090533e-18im
                       ⋮
   0.05320472714326833 + 1.3125707519115202e-18im
   0.05258500575198769 + 4.900122476070182e-18im
  0.053374187025433915 + 1.4500192394684467e-18im
   0.05280280912207337 + 4.790222663609598e-18im
   0.05353272287002102 + 1.5786388776856578e-18im
   0.05300592940528556 + 4.687377647481522e-18im
   0.05368103548081201 + 1.6989211160502906e-18im
   0.05319535891146354 + 4.591205423715975e-18im
   0.05381978096985631 + 1.811340158316642e-18im

source

MRIgeneralizedBloch.dG_o_dT2s_x_T2s_gaussian — Method

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

source

MRIgeneralizedBloch.dG_o_dT2s_x_T2s_lorentzian — Method

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

source

MRIgeneralizedBloch.dG_o_dT2s_x_T2s_superlorentzian — Method

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

source

MRIgeneralizedBloch.greens_gaussian — Method

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

source

MRIgeneralizedBloch.greens_lorentzian — Method

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

source

MRIgeneralizedBloch.greens_superlorentzian — Method

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

source

MRIgeneralizedBloch.hamiltonian_linear — Method

hamiltonian_linear(ω1, B1, ω0, T, m0s, R1, R2f, Rx, 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

Arguemnts

ω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
R1::Number: Apparent longitudinal relaxation rate of the free and semi-solid 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
R2s::Number: Transversal relaxationt rate of the semi-solid pool in 1/seconds; this number can be calcualated 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 calcualated 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 calcualated with the third function returned by precompute_R2sl
grad_type::grad_param: grad_m0s(), grad_R1(), 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> R1 = 1;

julia> R2f = 15;

julia> Rx = 30;

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, R1, R2f, Rx, R2s)) * m0
6-element StaticArrays.SVector{6, Float64} with indices SOneTo(6):
  0.0010646925712316103
  0.0
 -0.8957848933541458
  0.005125086137871529
  0.08119617921987109
  1.0

source

MRIgeneralizedBloch.interpolate_greens_function — Method

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

source

MRIgeneralizedBloch.precompute_R2sl — Method

precompute_R2sl(TRF_min, TRF_max, T2s_min, T2s_max, α_min, α_max, B1_min, B1_max[; greens=greens_superlorentzian])

Pre-compute and interpolate the linearized R2sl(TRF, α, B1, T2s) and its derivatives dR2sldB1(TRF, α, B1, T2s) and R2sldT2s(TRF, α, B1, T2s) 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, calulates the linearized R2sl that minimizes the error of zs at the end of the RF-pulse, and interpolates between the different samples.

Arguemnts

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
α_min::Number: lower bound of the flip angle range in radians
α_max::Number: upper bound of the flip angle range in radians
B1_min::Number: lower bound of the B1 range, normalized so that B1 = 1 corresponds to a perfectly calibrated RF field
B1_max::Number: upper bound of the B1 range, normalized so that B1 = 1 corresponds to a perfectly calibrated RF field

Optional:

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> (R2sl, dR2sldB1, R2sldT2s) = precompute_R2sl(100e-6, 1e-3, 5e-6, 15e-6, 0, π, 0.7, 1.3);


julia> (R2sl, dR2sldB1, R2sldT2s) = precompute_R2sl(100e-6, 1e-3, 5e-6, 15e-6, 0, π, 0.7, 1.3; greens=greens_gaussian);

source