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 with the exception of 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 of parameter is:

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

0.8894539621298889

The following plot visualizes the quality of the fit an 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.361652277711102

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

BIC_mono = n * log(RSS / n) + k * log(n)

-102.67982198678006

In this case, ΔAIC is per 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 with the exception of 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.19436834193025312

T1f = 1 / fit_Graham.param[2] # s

2.0217349657401127

T1s = 1 / fit_Graham.param[3] # s

0.24569338339800717

The following plot visualizes the quality of the fit an 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.14656319226751763

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

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

-18.580662734862543

and similarly for the BIC:

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

-16.142911085126144

Sled's MT model (unconstrained R₁ˢ)

MT_model = Sled()

The following parameters are hard-coded with the exception of 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.09937562692628771

T1f = 1 / fit_Sled.param[2] # s

1.70480391856932

T1s = 1 / fit_Sled.param[3] # s

0.11763471088913305

The following plot visualizes the quality of the fit an 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.13996537913499815

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

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

-19.732202281241385

and similarly for the BIC:

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

-17.294450631504986

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

MT_model = gBloch()

The following parameters are hard-coded with the exception of 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.12881834367305509

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

0.9742293629275287

The following plot visualizes the quality of the fit an 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.17054658200500056

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

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

-16.791868093373466

and similarly for the BIC:

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

-15.572992268505274

Generalized Bloch model (unconstrained R₁ˢ)

MT_model = gBloch()

The following parameters are hard-coded with the exception of 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.20903499148698299

T1f = 1 / fit_uncon.param[2] # s

2.0632291957979128

T1s = 1 / fit_uncon.param[3] # s

0.26253379951628936

The following plot visualizes the quality of the fit an 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.11479629852721095

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

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

-24.688098770327016

and similarly for the BIC:

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

-22.250347120590618

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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