GitHub

Global fit

First, we load the required packages and include some Helper functions

include("helper_functions.jl")

load all T₁-mapping methods:

include("T1_mapping_methods.jl")

and initialize the plot:

p = plot([0, 2], [0, 2], xlabel="simulated T1 (s)", ylabel="literature T1 (s)", legend=:topleft, label=:none, xlim=(0.5, 1.1), ylim=(0.5, 1.1))
marker_list = [(seq_type_i == :IR) ? (:circle) : ((seq_type_i == :LL) ? (:cross) : ((seq_type_i == :vFA) ? (:diamond) : ((seq_type_i == :SR) ? (:dtriangle) : (:x)))) for seq_type_i in seq_type]

Mono-exponential model

We simulate the mono-exponential model as an MT model with a vanishing semi-solid spin pool. In this case, the underlying MT model is irrelevant, and we choose Graham's model for speed purposes:

MT_model = Graham()

The following parameters are hard-coded, except for R1f, which serves as an initialization for the global fit.

m0s = 0
R1f = 1 / 1.084  # 1/s
R1s = R1f        # 1/s
R2f = 1 / 0.0769 # 1/s
T2s = 0
Rx = 0           # 1/s
ω0 = 0
B1 = 1

We define a model for the global fit. It takes the global set of parameters p as an input and returns a vector of T₁ estimates that correspond to the T₁ mapping methods described by the vector T1_functions

function model(iseq, p)
    R1f = p[1]
    T1 = similar(iseq, Float64)
    Threads.@threads for i in eachindex(iseq)
        T1[i] = T1_functions[iseq[i]](m0s, R1f, R2f, Rx, R1s, T2s)
    end
    return T1
end

Perform the fit and calculate the simulated T₁:

fit_mono = curve_fit(model, 1:length(T1_literature), T1_literature, [R1f], x_tol=1e-3)
T1_simulated = model(1:length(T1_literature), fit_mono.param)

For this model, the single fitted global parameter is:

T1f = 1 / fit_mono.param[1] # s
0.8894530882606787

The following plot visualizes the quality of the fit and replicates Fig. 1 in the manuscript:

scatter!(p, T1_simulated, T1_literature, label="mono-exponential model", markershape=marker_list, hover=seq_name)

Akaike (AIC) and Bayesian (BIC) information criteria

The information criteria depend on the number of measurements:

n = length(T1_literature)
25

the number of parameters:

k = length(fit_mono.param)
1

and the squared sum of the residuals:

RSS = norm(fit_mono.resid)^2
0.36165228007243183

As the AIC and BIC are only informative in the difference between two models, we use the mono-exponential model as the reference:

AIC_mono = n * log(RSS / n) + 2k
-103.8986976484162
BIC_mono = n * log(RSS / n) + k * log(n)
-102.679821823548

In this case, ΔAIC is by definition zero:

ΔAIC = n * log(RSS / n) + 2k - AIC_mono
0.0

as is ΔBIC:

ΔBIC = n * log(RSS / n) + k * log(n) - BIC_mono
0.0

Graham's MT model (unconstrained R₁ˢ)

MT_model = Graham()

The following parameters are hard-coded, except for m0s, R1f, and R1s, which serve as an initialization for the global fit.

m0s = 0.25
R1f = 1 / 1.84   # 1/s
R1s = 1 / 0.34   # 1/s
R2f = 1 / 0.0769 # 1/s
T2s = 12.5e-6    # s
Rx = 13.6        # 1/s
ω0 = 0
B1 = 1

We define a model for the global fit. It takes the global set of parameters p as an input and returns a vector of T₁ estimates that correspond to the T₁ mapping methods described by the vector T1_functions

function model(iseq, p)
    m0s, R1f, R1s = p
    T1 = similar(iseq, Float64)
    Threads.@threads for i in eachindex(iseq)
        T1[i] = T1_functions[iseq[i]](m0s, R1f, R2f, Rx, R1s, T2s)
    end
    return T1
end

Perform the fit and calculate the simulated T₁:

fit_Graham = curve_fit(model, 1:length(T1_literature), T1_literature, [m0s, R1f, R1s], x_tol=1e-3)
T1_simulated = model(1:length(T1_literature), fit_Graham.param)

For this model, the fitted global set of parameters is:

m0s = fit_Graham.param[1]
0.19446378469746337
T1f = 1 / fit_Graham.param[2] # s
2.026087604670749
T1s = 1 / fit_Graham.param[3] # s
0.24535283941180797

The following plot visualizes the quality of the fit and replicates Fig. 1 in the manuscript:

scatter!(p, T1_simulated, T1_literature, label="Graham's model (unconstrained R₁ˢ)", markershape=marker_list, hover=seq_name)

Note that clicking on a legend entry removes the fit from the plot. A double click on an entry selects only this particular fit.

Akaike (AIC) and Bayesian (BIC) information criteria

The information criteria depend on the number of measurements n, the number of parameters k, and the squared sum of the residuals RSS:

n = length(T1_literature)
25
k = length(fit_Graham.param)
3
RSS = norm(fit_Graham.resid)^2
0.14656258643959397

With this information, we can calculate the AIC difference to the mono-exponential model:

ΔAIC = n * log(RSS / n) + 2k - AIC_mono
-18.580766237337997

and similarly for the BIC:

ΔBIC = n * log(RSS / n) + k * log(n) - BIC_mono
-16.1430145876016

Sled's MT model (unconstrained R₁ˢ)

MT_model = Sled()

The following parameters are hard-coded, except for m0s, R1f, and R1s, which serve as an initialization for the global fit.

m0s = 0.25
R1f = 1 / 1.84   # 1/s
R1s = 1 / 0.34   # 1/s
R2f = 1 / 0.0769 # 1/s
T2s = 12.5e-6    # s
Rx = 13.6        # 1/s
ω0 = 0
B1 = 1

We define a model for the global fit. It takes the global set of parameters p as an input and returns a vector of T₁ estimates that correspond to the T₁ mapping methods described by the vector T1_functions

function model(iseq, p)
    m0s, R1f, R1s = p
    T1 = similar(iseq, Float64)
    Threads.@threads for i in eachindex(iseq)
        T1[i] = T1_functions[iseq[i]](m0s, R1f, R2f, Rx, R1s, T2s)
    end
    return T1
end

Perform the fit and calculate the simulated T₁:

fit_Sled = curve_fit(model, 1:length(T1_literature), T1_literature, [m0s, R1f, R1s], x_tol=1e-3)
T1_simulated = model(1:length(T1_literature), fit_Sled.param)

For this model, the fitted global set of parameters is:

m0s = fit_Sled.param[1]
0.09938220860166001
T1f = 1 / fit_Sled.param[2] # s
1.7045103081206696
T1s = 1 / fit_Sled.param[3] # s
0.11768341516509102

The following plot visualizes the quality of the fit and replicates Fig. 1 in the manuscript:

scatter!(p, T1_simulated, T1_literature, label="Sled's model (unconstrained R₁ˢ)", markershape=marker_list, hover=seq_name)

Note that clicking on a legend entry removes the fit from the plot. A double click on an entry selects only this particular fit.

Akaike (AIC) and Bayesian (BIC) information criteria

The information criteria depend on the number of measurements n, the number of parameters k, and the squared sum of the residuals RSS:

n = length(T1_literature)
25
k = length(fit_Sled.param)
3
RSS = norm(fit_Sled.resid)^2
0.13996520972762952

With this information, we can calculate the AIC difference to the mono-exponential model:

ΔAIC = n * log(RSS / n) + 2k - AIC_mono
-19.73223270329035

and similarly for the BIC:

ΔBIC = n * log(RSS / n) + k * log(n) - BIC_mono
-17.29448105355395

Generalized Bloch model (R₁ˢ = R₁ᶠ)

MT_model = gBloch()

The following parameters are hard-coded, except for m0s and R1f, which serve as an initialization for the global fit.

m0s = 0.139
R1f = 1 / 1.084  # 1/s
R1s = 1          # 1/s
R2f = 1 / 0.0769 # 1/s
T2s = 12.5e-6    # s
Rx = 23          # 1/s
ω0 = 0
B1 = 1

We define a model for the global fit. It takes the global set of parameters p as an input and returns a vector of T₁ estimates that correspond to the T₁ mapping methods described by the vector T1_functions

function model(iseq, p)
    m0s, R1f = p
    R1s = R1f
    T1 = similar(iseq, Float64)
    Threads.@threads for i in eachindex(iseq)
        T1[i] = T1_functions[iseq[i]](m0s, R1f, R2f, Rx, R1s, T2s)
    end
    return T1
end

Perform the fit and calculate the simulated T₁:

fit_constr = curve_fit(model, 1:length(T1_literature), T1_literature, [m0s, R1f], x_tol=1e-3)
T1_simulated = model(1:length(T1_literature), fit_constr.param)

For this model, the fitted global set of parameters is:

m0s = fit_constr.param[1]
0.12881793551863777
T1f = 1 / fit_constr.param[2] # s
0.9742290010278117

The following plot visualizes the quality of the fit and replicates Fig. 1 in the manuscript:

scatter!(p, T1_simulated, T1_literature, label="generalized Bloch model (R₁ˢ = R₁ᶠ constraint)", markershape=marker_list, hover=seq_name)

Note that clicking on a legend entry removes the fit from the plot. A double click on an entry selects only this particular fit.

Akaike (AIC) and Bayesian (BIC) information criteria

The information criteria depend on the number of measurements n, the number of parameters k, and the squared sum of the residuals RSS:

n = length(T1_literature)
25
k = length(fit_constr.param)
2
RSS = norm(fit_constr.resid)^2
0.17054658474454984

With this information, we can calculate the AIC difference to the mono-exponential model:

ΔAIC = n * log(RSS / n) + 2k - AIC_mono
-16.791867855021792

and similarly for the BIC:

ΔBIC = n * log(RSS / n) + k * log(n) - BIC_mono
-15.5729920301536

Generalized Bloch model (unconstrained R₁ˢ)

MT_model = gBloch()

The following parameters are hard-coded, except for m0s, R1f, and R1s, which serve as an initialization for the global fit.

m0s = 0.25
R1f = 1 / 1.84   # 1/s
R1s = 1 / 0.34   # 1/s
R2f = 1 / 0.0769 # 1/s
T2s = 12.5e-6    # s
Rx = 13.6        # 1/s
ω0 = 0
B1 = 1

We define a model for the global fit. It takes the global set of parameters p as an input and returns a vector of T₁ estimates that correspond to the T₁ mapping methods described by the vector T1_functions

function model(iseq, p)
    m0s, R1f, R1s = p
    T1 = similar(iseq, Float64)
    Threads.@threads for i in eachindex(iseq)
        T1[i] = T1_functions[iseq[i]](m0s, R1f, R2f, Rx, R1s, T2s)
    end
    return T1
end

Perform the fit and calculate the simulated T₁:

fit_uncon = curve_fit(model, 1:length(T1_literature), T1_literature, [m0s, R1f, R1s], x_tol=1e-3)
T1_simulated = model(1:length(T1_literature), fit_uncon.param)

For this model, the fitted global set of parameters is:

m0s = fit_uncon.param[1]
0.20903106717408043
T1f = 1 / fit_uncon.param[2] # s
2.063020626243208
T1s = 1 / fit_uncon.param[3] # s
0.262550660513381

The following plot visualizes the quality of the fit and replicates Fig. 1 in the manuscript:

scatter!(p, T1_simulated, T1_literature, label="generalized Bloch model (unconstrained R₁ˢ)", markershape=marker_list, hover=seq_name)

Note that clicking on a legend entry removes the fit from the plot. A double click on an entry selects only this particular fit.

Akaike (AIC) and Bayesian (BIC) information criteria

The information criteria depend on the number of measurements n, the number of parameters k, and the squared sum of the residuals RSS:

n = length(T1_literature)
25
k = length(fit_uncon.param)
3
RSS = norm(fit_uncon.resid)^2
0.11479629307483916

With this information, we can calculate the AIC difference to the mono-exponential model:

ΔAIC = n * log(RSS / n) + 2k - AIC_mono
-24.6881001209606

and similarly for the BIC:

ΔBIC = n * log(RSS / n) + k * log(n) - BIC_mono
-22.2503484712242

Data analysis

The span of T₁ estimates in the literature is

extrema(T1_literature)
(0.64, 1.0855)

and the coefficient of variation is

variation(T1_literature)
0.1408807728048427

Median absolute deviation

In the paper, the median absolute deviation from the mean value is used, as it is more robust to outliers compared to the mean absolute deviation or the standard deviation. Note, however, that the median absolute deviation of the mono-exponential fit is dominated by an outlier, artificially inflating the corresponding reduction.

A mono-exponential model explains the following fraction of the T₁ variability in the literature:

1 - mad(fit_mono.resid; center=mean(fit_mono.resid)) / mad(T1_literature; center=mean(T1_literature))
0.07118935312824093

Graham's spectral model without constraints on R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - mad(fit_Graham.resid; center=mean(fit_Graham.resid)) / mad(T1_literature; center=mean(T1_literature))
0.6071335860151669

Sled's model without constraints on R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - mad(fit_Sled.resid; center=mean(fit_Sled.resid)) / mad(T1_literature; center=mean(T1_literature))
0.5567822162017504

The generalized Bloch model constrained by R₁ˢ = R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - mad(fit_constr.resid; center=mean(fit_constr.resid)) / mad(T1_literature; center=mean(T1_literature))
0.4684203732531683

The generalized Bloch model without constraints on R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - mad(fit_uncon.resid; center=mean(fit_uncon.resid)) / mad(T1_literature; center=mean(T1_literature))
0.6171961683507089

Mean absolute deviation

For comparison, the mean absolute deviation is analyzed: A mono-exponential model explains the following fraction of the T₁ variability in the literature:

1 - mean(abs.(fit_mono.resid .- mean(fit_mono.resid))) / mean(abs.(T1_literature .- mean(T1_literature)))
0.023426403041203026

Graham's spectral model without constraints on R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - mean(abs.(fit_Graham.resid .- mean(fit_Graham.resid))) / mean(abs.(T1_literature .- mean(T1_literature)))
0.4834831199140681

Sled's model without constraints on R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - mean(abs.(fit_Sled.resid .- mean(fit_Sled.resid))) / mean(abs.(T1_literature .- mean(T1_literature)))
0.48265447199825884

The generalized Bloch model constrained by R₁ˢ = R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - mean(abs.(fit_constr.resid .- mean(fit_constr.resid))) / mean(abs.(T1_literature .- mean(T1_literature)))
0.42795873321794164

The generalized Bloch model without constraints on R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - mean(abs.(fit_uncon.resid .- mean(fit_uncon.resid))) / mean(abs.(T1_literature .- mean(T1_literature)))
0.5261587981865726

Standard deviation

For comparison, the standard deviation is analyzed: A mono-exponential model explains the following fraction of the T₁ variability in the literature:

1 - std(fit_mono.resid) / std(T1_literature)
0.01826093355907199

Graham's spectral model without constraints on R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - std(fit_Graham.resid) / std(T1_literature)
0.37502858305844344

Sled's spectral model without constraints on R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - std(fit_Sled.resid) / std(T1_literature)
0.38924942271187124

The generalized Bloch model constrained by R₁ˢ = R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - std(fit_constr.resid) / std(T1_literature)
0.32587111172659133

The generalized Bloch model without constraints on R₁ᶠ explains the following fraction of the T₁ variability in the literature:

1 - std(fit_uncon.resid) / std(T1_literature)
0.4468849451044564