Inversion Recovery Experiments

The following code replicates the NMR data analysis in Figs. 4-6 and complements the paper with additional analyses that are not shown in the paper in the interest of brevity.

For this analysis we need the following packages:

using MRIgeneralizedBloch
using DifferentialEquations
using LinearAlgebra
using LsqFit
using Statistics
import Pingouin
using Printf
using Formatting
using Plots

WARNING: Method definition anderson(Array{var"#s34", N} where N where var"#s34"<:Number) in module Pingouin at /home/runner/.julia/packages/Pingouin/L0EE5/src/distributions.jl:131 overwritten at /home/runner/.julia/packages/Pingouin/L0EE5/src/distributions.jl:140.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
┌ Warning: Replacing docs for `Pingouin.convert_effsize :: Tuple{Float64, String, String}` in module `Pingouin`
└ @ Base.Docs docs/Docs.jl:243
┌ Warning: Replacing docs for `Pingouin.compute_effsize_from_t :: Tuple{Float64}` in module `Pingouin`
└ @ Base.Docs docs/Docs.jl:243
┌ Warning: Replacing docs for `Pingouin.compute_esci :: Tuple{}` in module `Pingouin`
└ @ Base.Docs docs/Docs.jl:243

The raw data is stored in a separate github repository and the following functions return the URL to the individual files:

MnCl2_data(TRF_scale) = string("https://github.com/JakobAsslaender/MRIgeneralizedBloch_NMRData/blob/main/20210419_1mM_MnCl2/ja_IR_v2%20(", TRF_scale, ")/1/data.2d?raw=true")
BSA_data(TRF_scale)   = string("https://github.com/JakobAsslaender/MRIgeneralizedBloch_NMRData/blob/main/20210416_15%25BSA_2ndBatch/ja_IR_v2%20(", TRF_scale, ")/1/data.2d?raw=true");

which can be loaded with functions implemented in this file:

include(string(pathof(MRIgeneralizedBloch), "/../../docs/src/load_NMR_data.jl"));

MnCl$_2$ Sample

$T_2^{*,f}$ Estimation

We estimate $T_2^{*,f}$ by fitting a mono-exponential decay curve to the FID of the acquisition with $T_\text{RF} = 22.8$μs and $T_\text{i} = 5$s.

M = load_Data(MnCl2_data(1))
M = M[:,1]; # select Tᵢ = 5s

The data was measured at the following time points in units of seconds:

T_dwell = 100e-6 # s
TE = T_dwell * ((1:length(M)) .+ 7) # s

0.0008:0.0001:1.6378

Note that the signal is an FID, so the phrase echo time is a bit misleading.

The function curve_fit from the LsqFit.jl package is only implemented for real-valued models. To accommodate this, we need to split the data into its real and imaginary part:

TEreal = [TE;TE]
Mreal = [real(M);imag(M)];

Here, we are using a simple mono-exponential model with a complex-valued scaling factor p[1] + 1im p[2], the decay time $T_2^{*,f} =$ p[3], and the Larmor frequency p[4]:

FID_model(t, p) = @. [p[1] * exp(- t[1:end ÷ 2] / p[3]) * cos(p[4] * t[1:end ÷ 2]); p[2] * exp(- t[end ÷ 2 + 1:end] / p[3]) * sin(p[4] * t[end ÷ 2 + 1:end])];

Fitting this model to the NMR data estimates $T_2^{*,f}$:

fit = curve_fit(FID_model, TEreal, Mreal, [1, 1, 0.1, 0])
T₂star_MnCl2 = fit.param[3] # s

0.07497933747410264

seconds and its uncertainty (also in units of seconds)

stderror(fit)[3] # s

1.7377694953694763e-5

Visually, the plot and the data show good agreement:

Mfitted = FID_model(TEreal, fit.param)
Mfitted = Mfitted[1:end÷2] + 1im * Mfitted[end÷2+1:end]
p = plot(xlabel="TE [s]", ylabel="|FID(TE)| [a.u.]")
plot!(p, TE, abs.(M), label="data")
plot!(p, TE, abs.(Mfitted), label=@sprintf("fit with T₂* = %2.3f ms", 1e3 * T₂star_MnCl2))

The relative residual norm of the fit, i.e. $||\text{residual}||_2/||M||_2$ is

norm(fit.resid) / norm(M)

0.013503650656413443

Despite its small $\ell_2$-norm, the Shapiro-Wilk test indicates that the residual is not Gaussian or normal distributed at a significance level of α=0.05

Pingouin.normality(fit.resid, α=0.05)

1×3 DataFrame

Row	W	pval	normal
	Float64	Float64	Bool
1	0.483365	6.38071e-131	false

We note that mono-exponential $T_2^*$ decays assume a Lorentzian distributed magnetic field, which is in general an assumption rather than a well-founded theory.

Mono-Exponential IR Model

We performed several experiments in which we inverted the thermal equilibrium magnetization with rectangular π-pulses with the following pulse durations (in seconds):

Tʳᶠmin = 22.8e-6 # s - shortest Tʳᶠ possible on the NMR
TRF_scale = [1;2;5:5:40] # scaling factor
Tʳᶠ = TRF_scale * Tʳᶠmin # s

10-element Vector{Float64}:
 2.28e-5
 4.56e-5
 0.00011399999999999999
 0.00022799999999999999
 0.00034199999999999996
 0.00045599999999999997
 0.00057
 0.0006839999999999999
 0.000798
 0.0009119999999999999

and acquired inversion recovery data at exponentially spaced inversion times (in seconds):

Tᵢ = exp.(range(log(3e-3), log(5), length=20)) # s
Tᵢ .+= 12 * Tʳᶠmin + (13 * 15.065 - 5) * 1e-6 # s - correction factors

20-element Vector{Float64}:
 0.003464444999999999
 0.0048973889644219615
 0.007014775730568361
 0.010143528025680238
 0.014766722559942762
 0.02159817666227864
 0.03169266107267046
 0.046608755653003126
 0.06864949280054741
 0.10121794370375319
 0.1493426496577329
 0.22045402459049154
 0.32553160469379794
 0.4807992795321372
 0.7102302468802864
 1.049248454199443
 1.550198026187646
 2.2904251533959354
 3.3842202786824003
 5.0004644449999995

We calculate the Rabi frequencies of the RF pulses and a finer grid of $T_\text{i}$ to plot the IR model:

ω₁ = π ./ Tʳᶠ # rad/s
Tᵢplot = exp.(range(log(Tᵢ[1]), log(Tᵢ[end]), length=500)); # s

After loading and normalizing the data

M = zeros(Float64, length(Tᵢ), length(TRF_scale))
for i = 1:length(TRF_scale)
    M[:,i] = load_first_datapoint(MnCl2_data(TRF_scale[i]))
end
M ./= maximum(M);

we analyze each inversion recovery curve that corresponds to a different $T_\text{RF}$ separately. This allows us to fit a simple mono-exponential model

standard_IR_model(t, p) = @. p[1] - p[3] * exp(- t * p[2]);

where p[1] is the thermal equilibrium magnetization, p[2] $= T_1$, and p[1] - p[3] is the magnetization right after the inversion pulse or, equivalently, Minv = p[1] / p[3] - 1 is the inversion efficiency, which is 1 for an ideal π-pulse and smaller otherwise. The parameters are initialized with

p0 = [1.0, 1.0, 2.0];

and we can loop over $T_\text{RF}$ to perform the fits:

R₁ = similar(M[1,:])
Minv = similar(R₁)
residual = similar(R₁)
p = plot(xlabel="Tᵢ [s]", ylabel="zᶠ(Tʳᶠ, Tᵢ) [a.u.]")
for i = 1:length(TRF_scale)
    Mi = @view M[:,i]

    fit = curve_fit(standard_IR_model, Tᵢ, Mi, p0)

    R₁[i] = fit.param[2]
    Minv[i] = fit.param[3] / fit.param[1] - 1

    residual[i] = norm(fit.resid) / norm(Mi)

    scatter!(p, Tᵢ, Mi, label=@sprintf("Tʳᶠ = %1.2es - data", Tʳᶠ[i]), color=i)
    plot!(p, Tᵢplot, standard_IR_model(Tᵢplot, fit.param), label=@sprintf("fit with R₁ = %.3f/s; MInv = %.3f", R₁[i], Minv[i]), color=i)
end

Here, the data measured with different $T_\text{RF}$ are indicated by markers in different colors, and the corresponding fits are the line plots in the same color. The fitted parameters are denoted in the legend. In the paper, we highlight the estimated inversion efficiency and the relaxation rate of the dataset acquired with $T_\text{RF}=22.8$μs:

Minv[1]

0.9303588578422617

R₁[1] # 1/s

1.4771257292672177

and of the dataset acquired with $T_\text{RF}=912$μs:

Minv[end]

0.9092660344729475

R₁[end] # 1/s

1.479057638066805

The mean value of all R₁ estimates is

mean(R₁) # 1/s

1.478741342058

1/s and its standard deviation in units of 1/s is

std(R₁) # 1/s

0.001590791201701839

The relative residual norm of the fits is on average

mean(residual)

0.0013767952719624893

Global IR Fit

As an alternative to individual fits to the inversion recovery curves with different $T_\text{RF}$, we can also perform a global fit that accounts for the $T_2^{*,f}$ decay during the inversion pulse. The model first simulates the $T_2^{*,f}$ decay during the inversion pulse, followed by $T_1$ recovery:

function Bloch_IR_model(p, Tʳᶠ, Tᵢ, T2)
    (m0, m0_inv, R₁) = p
    R2 = 1 / T2

    M = zeros(Float64, length(Tᵢ), length(Tʳᶠ))
    for i = 1:length(Tʳᶠ)
        # simulate inversion pulse
        ω₁ = π / Tʳᶠ[i]
        H = [-R2 -ω₁ 0 ;
              ω₁ -R₁ R₁;
               0   0 0 ]

        m_inv = m0_inv * (exp(H * Tʳᶠ[i]) * [0,1,1])[2]

        # simulate T1 recovery
        H = [-R₁ R₁*m0;
               0     0]

        for j = 1:length(Tᵢ)
            M[j,i] = m0 * (exp(H .* (Tᵢ[j] - Tʳᶠ[i] / 2)) * [m_inv,1])[1]
        end
    end
    return vec(M)
end;

We use the previously estimated $T_2^{*,f}$ value for the fit:

fit = curve_fit((x, p) -> Bloch_IR_model(p, Tʳᶠ, Tᵢ, T₂star_MnCl2), 1:length(M), vec(M), [ 1, .8, 1])

p = plot(xlabel="Tᵢ [s]", ylabel="zᶠ(Tʳᶠ, Tᵢ) [a.u.]")
for i=1:length(Tʳᶠ)
    scatter!(p, Tᵢ, M[:,i], label=@sprintf("Tʳᶠ = %1.2es", Tʳᶠ[i]), color=i)
    plot!(p, Tᵢplot, Bloch_IR_model(fit.param, Tʳᶠ[i], Tᵢplot, T₂star_MnCl2), label=@sprintf("Tʳᶠ = %1.2es", Tʳᶠ[i]), color=i)
end

With this global fit, we get a very similar relaxation rate in units of 1/s

R₁_MnCl2 = fit.param[3] # 1/s

1.4787326527483415

with an uncertainty (also in units of 1/s) of

stderror(fit)[3] # 1/s

0.0025253461723147427

Note that the relative residual norm is somewhat increased compared to individual fits to each inversion recovery curve:

norm(fit.resid) / norm(M)

0.0069842810791774365

Bovine Serum Albumin Sample

$T_2^{*,f}$ Estimation

We repeat the $T_2^{*,f}$ estimation for the bovine serum albumin (BSA) sample by fitting a mono-exponential decay curve to the FID of the acquisition with $T_\text{RF} = 22.8$μs and $T_\text{i} = 5$s.

M = load_Data(BSA_data(1));
M = M[:,1] # select Tᵢ = 5s
Mreal = [real(M);imag(M)]

fit = curve_fit(FID_model, TEreal, Mreal, [1.0, 1.0, .1, 0.0]);

The estimated $T_2^{*,f}$ of the BSA sample is

T₂star_BSA = fit.param[3] # s

0.05240554901160865

seconds with an uncertainty of

stderror(fit)[3] # s

7.0295104330382624e-6

seconds.

Visually, the plot and the data align well for the BSA sample, too:

Mfitted = FID_model(TEreal, fit.param)
Mfitted = Mfitted[1:end÷2] + 1im * Mfitted[end÷2+1:end]
p = plot(xlabel="TE [s]", ylabel="|FID(TE)| [a.u.]")
plot!(p, TE, abs.(M), label="data")
plot!(p, TE, abs.(Mfitted), label=@sprintf("fit with T₂* = %2.3f ms", 1e3 * T₂star_BSA))

The relative residual norm ($||\text{residual}||_2/||M||_2$) is

norm(fit.resid) / norm(M)

0.006257402534109852

Despite the small residual, the Shapiro-Wilk test indicates that the residual is not normal distributed for this sample either:

Pingouin.normality(fit.resid, α=0.05)

1×3 DataFrame

Row	W	pval	normal
	Float64	Float64	Bool
1	0.534965	7.08003e-128	false

Mono-Exponential IR Model

We also fit a mono-exponential model to each inversion recovery curve of the BSA data:

M = zeros(Float64, length(Tᵢ), length(TRF_scale))
for i = 1:length(TRF_scale)
    M[:,i] = load_first_datapoint(BSA_data(TRF_scale[i]))
end
M ./= maximum(M)


p = plot(xlabel="Tᵢ [s]", ylabel="zᶠ(Tʳᶠ, Tᵢ) [a.u.]")
for i = 1:length(TRF_scale)
    Mi = @view M[:,i]

    fit = curve_fit(standard_IR_model, Tᵢ, Mi, p0)

    R₁[i] = fit.param[2]
    Minv[i] = fit.param[3] / fit.param[1] - 1
    residual[i] = norm(fit.resid) / norm(Mi)

    scatter!(p, Tᵢ, Mi, label=@sprintf("Tʳᶠ = %1.2es - data", Tʳᶠ[i]), color=i)
    plot!(p, Tᵢplot, standard_IR_model(Tᵢplot, fit.param), label=@sprintf("fit with R₁ = %.3f/s; MInv = %.3f", R₁[i], Minv[i]), color=i)
end

Zooming into the early phase of the recovery curve reveals the poor fit quality, in particular for long $T_\text{RF}$. This is also reflected by a substantially larger relative residual norm compared to the MnCl$_2$ sample:

mean(residual)

0.01924546690167816

In the paper, we highlight the estimated inversion efficiency and the relaxation rate of the dataset acquired with $T_\text{RF}=22.8$μs

Minv[1]

0.9324506391277547

R₁[1] # 1/s

1.1962243323294965

and of the dataset acquired with $T_\text{RF}=912$μs

Minv[end]

0.8626963380362451

R₁[end] # 1/s

1.2614605480338448

The mean value of all R₁ estimates is

mean(R₁) # 1/s

1.2417851203291097

and its standard deviation is substantially larger compared to the same fit of the MnCl$_2$ sample:

std(R₁) # 1/s

0.023621750246010924

Global IR Fit - Generalized Bloch Model

In order to repeat the global fit that includes all $T_\text{RF}$ values, we have to account for the spin dynamics in the semi-solid pool during the RF-pulse. First, we do this with the proposed generalized Bloch model:

function gBloch_IR_model(p, G, Tʳᶠ, TI, R2f)
    (m0, m0f_inv, m0s, R₁, T₂ˢ, Rx) = p
    m0f = 1 - m0s
    ω₁ = π ./ Tʳᶠ

    m0vec = [0, 0, m0f, m0s, 1]
    m_fun(p, t; idxs=nothing) = typeof(idxs) <: Number ? 0.0 : zeros(5)


    H = [-R₁-m0s*Rx     m0f*Rx R₁*m0f;
             m0s*Rx -R₁-m0f*Rx R₁*m0s;
              0          0         0 ]

    M = zeros(Float64, length(TI), length(Tʳᶠ))
    for i = 1:length(Tʳᶠ)
        param = (ω₁[i], 1, 0, m0s, R₁, R2f, Rx, R₁, T₂ˢ, G)
        prob = DDEProblem(apply_hamiltonian_gbloch!, m0vec, m_fun, (0.0, Tʳᶠ[i]), param)
        m = solve(prob).u[end]

        for j = 1:length(TI)
            M[j,i] = m0 * (exp(H .* (TI[j] - Tʳᶠ[i] / 2)) * [m0f_inv * m[3],m[4],1])[1]
        end
    end
    return vec(M)
end;

Here, we use assume a super-Lorentzian lineshape, whose Green's function is interpolated to speed up the fitting routine:

T₂ˢ_min = 5e-6 # s
G_superLorentzian = interpolate_greens_function(greens_superlorentzian, 0, maximum(Tʳᶠ)/T₂ˢ_min);

The fit is initialized with p0 = [m0, m0f_inv, m0_s, R₁, T2_s, Rx] and we set some reasonable bounds to the fitted parameters:

p0   = [  1, 0.932,  0.1,   1, 10e-6, 50]
pmin = [  0, 0.100,   .0, 0.3,  1e-9, 10]
pmax = [Inf,   Inf,  1.0, Inf, 20e-6,1e3]

fit_gBloch = curve_fit((x, p) -> gBloch_IR_model(p, G_superLorentzian, Tʳᶠ, Tᵢ, 1/T₂star_BSA), [], vec(M), p0, lower=pmin, upper=pmax);

Visually, the plot and the data align well:

p = plot(xlabel="Tᵢ [s]", ylabel="zᶠ(Tʳᶠ, Tᵢ) [a.u.]")
for i=1:length(Tʳᶠ)
    scatter!(p, Tᵢ, M[:,i], label=@sprintf("Tʳᶠ = %1.2es", Tʳᶠ[i]), color=i)
    plot!(p, Tᵢplot, gBloch_IR_model(fit_gBloch.param, G_superLorentzian, Tʳᶠ[i], Tᵢplot, 1/T₂star_BSA), label=@sprintf("Tʳᶠ = %1.2es", Tʳᶠ[i]), color=i)
end

which becomes particularly apparent when zooming into the beginning of the inversion recovery curves. Further, the relative residual norm is much smaller compared to the mono-exponential fit:

norm(fit_gBloch.resid) / norm(M)

0.004136056037877003

The estimated parameters are

m0 = fit_gBloch.param[1]

1.095292088036571

Minv = fit_gBloch.param[2]

0.9298337864526318

m0s = fit_gBloch.param[3]

0.08544861731233476

R₁ = fit_gBloch.param[4] # 1/s

1.1977731971880268

T₂ˢ = 1e6fit_gBloch.param[5] # μs

12.95430561288214

Rx = fit_gBloch.param[6] # 1/s

71.25330815815983

with the uncertainties (in the same order)

stderror(fit_gBloch)

6-element Vector{Float64}:
 0.001385902405274573
 0.0012241787869572095
 0.0009954674893831012
 0.0014821532962963702
 6.07295159602924e-7
 1.9919636474365263

Global IR Fit - Graham's Spectral Model

For comparison, we repeat the same fit with Graham's spectral model:

function Graham_IR_model(p, Tʳᶠ, TI, R2f)
    (m0, m0f_inv, m0s, R₁, T₂ˢ, Rx) = p
    m0f = 1 - m0s
    ω₁ = π ./ Tʳᶠ

    m0vec = [0, 0, m0f, m0s, 1]

    H = [-R₁-m0s*Rx     m0f*Rx R₁*m0f;
             m0s*Rx -R₁-m0f*Rx R₁*m0s;
              0          0         0 ]

    M = zeros(Float64, length(TI), length(Tʳᶠ))
    for i = 1:length(Tʳᶠ)
        param = (ω₁[i], 1, 0, Tʳᶠ[i], m0s, R₁, R2f, Rx, R₁, T₂ˢ)
        prob = ODEProblem(apply_hamiltonian_graham_superlorentzian!, m0vec, (0.0, Tʳᶠ[i]), param)
        m = solve(prob).u[end]

        for j = 1:length(TI)
            M[j,i] = m0 * (exp(H .* (TI[j] - Tʳᶠ[i] / 2)) * [m0f_inv * m[3],m[4],1])[1]
        end
    end
    return vec(M)
end

fit_Graham = curve_fit((x, p) -> Graham_IR_model(p, Tʳᶠ, Tᵢ, 1/T₂star_BSA), [], vec(M), p0, lower=pmin, upper=pmax);

Visually, the plot and the data align substantially worse:

p = plot(xlabel="Tᵢ [s]", ylabel="zᶠ(Tʳᶠ, Tᵢ) [a.u.]")
for i=1:length(Tʳᶠ)
    scatter!(p, Tᵢ, M[:,i], label=@sprintf("Tʳᶠ = %1.2es", Tʳᶠ[i]), color=i)
    plot!(p, Tᵢplot, Graham_IR_model(fit_Graham.param, Tʳᶠ[i], Tᵢplot, 1/T₂star_BSA), label=@sprintf("Tʳᶠ = %1.2es", Tʳᶠ[i]), color=i)
end

which becomes particularly apparent when zooming into the beginning of the inversion recovery curves. Further, the relative residual norm is much larger compared to the generalized Bloch fit:

norm(fit_Graham.resid) / norm(M)

0.0183281513944232

The estimated parameters are

m0 = fit_Graham.param[1]

1.083303509574129

Minv = fit_Graham.param[2]

0.9405873264118447

m0s = fit_Graham.param[3]

0.06992529329580319

R₁ = fit_Graham.param[4] # 1/s

1.1790784853335596

T₂ˢ = 1e6fit_Graham.param[5] # μs

2.895146002636496

Rx = fit_Graham.param[6] # 1/s

73.86783909974669

with the uncertainties (in the same order)

stderror(fit_Graham)

6-element Vector{Float64}:
 0.008271785461146819
 0.006169372316640898
 0.006411522420095765
 0.006696597594263232
 5.553347041214335e-6
 8.551126824965007

Global IR Fit - Sled's Model

We also performed the fit with Sled's model:

function Sled_IR_model(p, G, Tʳᶠ, TI, R2f)
    (m0, m0f_inv, m0s, R₁, T₂ˢ, Rx) = p
    m0f = 1 - m0s
    ω₁ = π ./ Tʳᶠ

    m0vec = [0, 0, m0f, m0s, 1]

    H = [-R₁-m0s*Rx     m0f*Rx R₁*m0f;
             m0s*Rx -R₁-m0f*Rx R₁*m0s;
              0          0         0 ]

    M = zeros(Float64, length(TI), length(Tʳᶠ))
    for i = 1:length(Tʳᶠ)
        param = (ω₁[i], 1, 0, m0s, R₁, R2f, Rx, R₁, T₂ˢ, G)
        prob = ODEProblem(apply_hamiltonian_sled!, m0vec, (0.0, Tʳᶠ[i]), param)
        m = solve(prob).u[end]

        for j = 1:length(TI)
            M[j,i] = m0 * (exp(H .* (TI[j] - Tʳᶠ[i] / 2)) * [m0f_inv * m[3],m[4],1])[1]
        end
    end
    return vec(M)
end

fit_Sled = curve_fit((x, p) -> Sled_IR_model(p, G_superLorentzian, Tʳᶠ, Tᵢ, 1/T₂star_BSA), [], vec(M), p0, lower=pmin, upper=pmax);

Visually, the plot and the data do not align well either:

p = plot(xlabel="Tᵢ [s]", ylabel="zᶠ(Tʳᶠ, Tᵢ) [a.u.]")
for i=1:length(Tʳᶠ)
    scatter!(p, Tᵢ, M[:,i], label=@sprintf("Tʳᶠ = %1.2es", Tʳᶠ[i]), color=i)
    plot!(p, Tᵢplot, Sled_IR_model(fit_Sled.param, G_superLorentzian, Tʳᶠ[i], Tᵢplot, 1/T₂star_BSA), label=@sprintf("Tʳᶠ = %1.2es", Tʳᶠ[i]), color=i)
end

which becomes particularly apparent when zooming into the beginning of the inversion recovery curves. Further, the relative residual norm is also large compared to the generalized Bloch fit:

norm(fit_Sled.resid) / norm(M)

0.01963558188079065

The estimated parameters are

m0 = fit_Sled.param[1]

1.0833069200503345

Minv = fit_Sled.param[2]

0.9368004034508981

m0s = fit_Sled.param[3]

0.07153499085682131

R₁ = fit_Sled.param[4] # 1/s

1.1781823629093933

T₂ˢ = 1e6fit_Sled.param[5] # μs

5.026439766352263

Rx = fit_Sled.param[6] # 1/s

60.019636791183366

with the uncertainties (in the same order)

stderror(fit_Sled)

6-element Vector{Float64}:
 0.005856522519818164
 0.006385157722960878
 0.0038387727466202308
 0.007379838048503139
 1.149119979136867e-6
 8.072347778316988

Analysis of the Residuals

In order to visualize how well the three models align with the data at different $T_\text{RF}$, we calculate the $\ell_2$-norm of the residuals after subtracting the modeled from the measured signal and normalize it by the $\ell_2$-norm of the signal:

resid_gBlo = similar(Tʳᶠ)
resid_Sled = similar(Tʳᶠ)
resid_Grah = similar(Tʳᶠ)
for i=1:length(Tʳᶠ)
    resid_gBlo[i] = norm(gBloch_IR_model(fit_gBloch.param, G_superLorentzian, Tʳᶠ[i], Tᵢ, 1/T₂star_BSA) .- M[:,i]) / norm(M[:,i])
    resid_Grah[i] = norm(Graham_IR_model(fit_Graham.param, Tʳᶠ[i], Tᵢ, 1/T₂star_BSA)                    .- M[:,i]) / norm(M[:,i])
    resid_Sled[i] = norm(Sled_IR_model(  fit_Sled.param,   G_superLorentzian, Tʳᶠ[i], Tᵢ, 1/T₂star_BSA) .- M[:,i]) / norm(M[:,i])
end

p = plot(xlabel="Tʳᶠ [s]", ylabel="relative residual")
scatter!(p, Tʳᶠ, resid_gBlo, label="generalized Bloch model")
scatter!(p, Tʳᶠ, resid_Grah, label="Graham's spectral model")
scatter!(p, Tʳᶠ, resid_Sled, label="Sled's model")

This analysis examines the residuals from the actual fits, i.e. it uses the biophysical parameters of respective fit to model the signal. The disadvantage of this approach is that residuals at long $T_\text{RF}$ are negatively affected by the poor fits of Graham's and Sled's models at short $T_\text{RF}$. This problem is overcome by subtracting the measured signal from signal that is simulated with the biophysical parameters that were estimated by fitting the generalized Bloch model:

for i=1:length(Tʳᶠ)
    resid_gBlo[i] = norm(gBloch_IR_model(fit_gBloch.param, G_superLorentzian, Tʳᶠ[i], Tᵢ, 1/T₂star_BSA) .- M[:,i]) / norm(M[:,i])
    resid_Grah[i] = norm(Graham_IR_model(fit_gBloch.param, Tʳᶠ[i], Tᵢ, 1/T₂star_BSA)                    .- M[:,i]) / norm(M[:,i])
    resid_Sled[i] = norm(Sled_IR_model(  fit_gBloch.param, G_superLorentzian, Tʳᶠ[i], Tᵢ, 1/T₂star_BSA) .- M[:,i]) / norm(M[:,i])
end

p = plot(xlabel="Tʳᶠ [s]", ylabel="relative residual")
scatter!(p, Tʳᶠ, resid_gBlo, label="generalized Bloch model")
scatter!(p, Tʳᶠ, resid_Grah, label="Graham's spectral model")
scatter!(p, Tʳᶠ, resid_Sled, label="Sled's model")

One can observe reduced residuals for Graham's and Sled's models for long $T_\text{RF}$ as a trade off for larger residuals at short $T_\text{RF}$. Yet, the residuals at long $T_\text{RF}$ are still substantially larger compared to ones of the generalized Bloch model.

This page was generated using Literate.jl.