Continuous Wave Saturation Experiments

The following code analyzes data from a steady-state experiment similar to the original work of Henkelman et al. In this experiment, the magnetization of the coupled spin system is saturated with off-resonant continuous waves of the exponentially spaced frequencies:

Δ = exp.(range(log(0.01e3), log(100e3), length=20)) * 2π # rad/s

20-element Vector{Float64}:
     62.83185307179588
    102.02490149810555
    165.66566187064444
    269.00404822785634
    436.8025162599944
    709.2697655220819
   1151.6957470645957
   1870.0967647060686
   3036.619626561797
   4930.79231537499
   8006.5058674764105
  13000.777989786218
  21110.360891171575
  34278.51297096368
  55660.65201622967
  90380.47203203167
 146757.70816250978
 238301.75281093005
 386948.8431222893
 628318.5307179587

and the amplitudes:

ω1_dB = -60:5:-5 # dB
ω1 = @. 10^(ω1_dB / 20) * π / 2 / 11.4e-6 # rad/s

12-element Vector{Float64}:
   137.78915147323656
   245.02761099159085
   435.7275555173729
   774.845340363136
  1377.8915147323655
  2450.276109915909
  4357.275555173728
  7748.453403631359
 13778.915147323654
 24502.76109915909
 43572.75555173729
 77484.53403631358

The waves were applied for 7 seconds to ensure a steady state. Thereafter, the magnetization was excited with a π/2-pulse and an FID was acquired. The repetition times was 30s to ensure full recovery to thermal equilibrium.

We fit the data with Henkelman's closed form solution to this steady-state problem while assuming a Lorentzian lineshape for the free spin pool, and different lineshapes for the semi-solid spin pool:

g_Lorentzian(Δ, T2) = T2 / π / (1 + (T2 * Δ)^2)
g_Gaussian(Δ, T2) = T2 / sqrt(2π) * exp(-(T2 * Δ)^2 / 2)
g_superLorentzian(Δ, T2) = T2 * sqrt(2 / π) * quadgk(ct -> exp(- 2 *  (T2 * Δ / (3 * ct^2 - 1))^2) / abs(3 * ct^2 - 1), 0, sqrt(1 / 3), 1)[1];

For this data analysis we need the following packages:

using MRIgeneralizedBloch
using QuadGK
using LsqFit
using LinearAlgebra
using Statistics
using Printf
using Plots

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

MnCl2_data(ω1_dB) = string("https://github.com/JakobAsslaender/MRIgeneralizedBloch_NMRData/blob/main/20210419_1mM_MnCl2/ja_PreSat_v2%20(", ω1_dB, ")/1/data.2d?raw=true")
BSA_data(ω1_dB)   = string("https://github.com/JakobAsslaender/MRIgeneralizedBloch_NMRData/blob/main/20210416_15%25BSA_2ndBatch/ja_PreSat_v2%20(", ω1_dB, ")/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"));

We store the off-resonance frequencies and wave amplitudes in a single matrix for convenience:

x = zeros(Float64, length(ω1) * length(Δ), 2)
x[:,1] = repeat(Δ, length(ω1))
x[:,2] = vec(repeat(ω1, 1, length(Δ))');

MnCl$_2$ Sample

We load the first data point of each FID:

M = zeros(Float64, length(Δ), length(ω1))
for i = 1:length(ω1_dB)
    M[:,i] = load_first_datapoint(MnCl2_data(ω1_dB[i]); set_phase=:abs)
end
M ./= maximum(M);

In contrast to the inversion-recovery experiment, the phase of the signal was not stable. Therefore, we took the absolute value of the signal by setting the flag set_phase=:abs.

The MnCl$_2$-data can be described with a single compartment model:

function single_compartment_model(x, p)
    (m0, R1, T2) = p

    Δ  = @view x[:,1]
    ω1 = @view x[:,2]

    Rrf = @. π * ω1^2 * g_Lorentzian(Δ, T2)
    m = @. m0 * R1 / (R1 + Rrf)
    return m
end;

(cf. Eqs. (14) and (15) in the paper).

As this model is merely a function of the relaxation times T₁ and T₂, we forgo a fitting routine and use the estimates from the Inversion Recovery Experiments instead:

R1 = 1.479 # 1/s
T2 = 0.075 # s

Visually, this model describes the data well:

p = plot(xlabel="Δ [rad/s]", ylabel="M / max(M)", xaxis=:log, legend=:none)
[scatter!(p, Δ, M[:,i], color=i) for i=1:length(ω1)]
[plot!(p, Δ, reshape(single_compartment_model(x, [1,R1,T2]), length(Δ), length(ω1))[:,i], color=i) for i=1:length(ω1)]

Bovine Serum Albumin Sample

We acquired the same data for the BSA sample, which we load:

for i = 1:length(ω1_dB)
    M[:,i] = load_first_datapoint(BSA_data(ω1_dB[i]); set_phase=:abs)
end
M ./= maximum(M);

We model the steady-state magnetization as described by Henkelman et al.:

function Henkelman_model(x, p; lineshape=:superLorentzian)
    (m0, m0s, R1f, R1s, T2f, T2s, Rx) = p

    m0s /= 1 - m0s # switch from m0s + m0f = 1 to m0f = 1 normalization

    Δ  = @view x[:,1]
    ω1 = @view x[:,2]

    Rrf_f = @. π * ω1^2 * g_Lorentzian(Δ, T2f)

    if lineshape == :Lorentzian
        Rrf_s = @. π * ω1^2 * g_Lorentzian(Δ, T2s)
    elseif lineshape == :Gaussian
        Rrf_s = @. π * ω1^2 * g_Gaussian(Δ, T2s)
    elseif lineshape == :superLorentzian
        Rrf_s = @. π * ω1^2 * g_superLorentzian(Δ, T2s)
    end

    m = @. m0 * (R1s * Rx * m0s + Rrf_s * R1f + R1f * R1s + R1f * Rx) / ((R1f + Rrf_f + Rx * m0s) * (R1s + Rrf_s + Rx) - Rx^2 * m0s)
    return m
end;

Here, we use a fitting routine to demonstrate the best possible fit with each of the three lineshapes. We define an initialization for the fitting routine p0 = [m0, m0s, R1f, R1s, T2f, T2s, Rx] and set some reasonable bounds:

p0   = [  1,0.01,   1,   5,0.052,  1e-5, 40]
pmin = [  0,   0,   0,   0,0.052,  1e-6,  1]
pmax = [Inf,   1, Inf, Inf,0.052, 10e-3,100];

Note that we fixed T₂ᶠ = 52ms to the value estimated with the Inversion Recovery Experiments as T₂ᶠ is poorly defined by this saturation experiment.

Super-Lorentzian Lineshape

Fitting the model with a super-Lorentzian lineshape to the data achieves good concordance:

fit = curve_fit((x, p) -> Henkelman_model(x, p; lineshape=:superLorentzian), x, vec(M), p0, lower=pmin, upper=pmax)
fit_std = stderror(fit)

p = plot(xlabel="Δ [rad/s]", ylabel="M / max(M)", xaxis=:log, legend=:none)
[scatter!(p, Δ, M[:,i], color=i) for i=1:length(ω1)]
[plot!(p, Δ, reshape(Henkelman_model(x, fit.param), length(Δ), length(ω1))[:,i], color=i) for i=1:length(ω1)]

Lorentzian Lineshape

The Lorentzian lineshape, on the other hand, does not fit the data well:

fit = curve_fit((x, p) -> Henkelman_model(x, p; lineshape=:Lorentzian), x, vec(M), p0, lower=pmin, upper=pmax)
fit_std = stderror(fit)

p = plot(xlabel="Δ [rad/s]", ylabel="M / max(M)", xaxis=:log, legend=:none)
[scatter!(p, Δ, M[:,i], color=i) for i=1:length(ω1)]
[plot!(p, Δ, reshape(Henkelman_model(x, fit.param), length(Δ), length(ω1))[:,i], color=i) for i=1:length(ω1)]

Gaussian Lineshape

And the Gaussian lineshape does not not fit the data well either:

fit = curve_fit((x, p) -> Henkelman_model(x, p; lineshape=:Lorentzian), x, vec(M), p0, lower=pmin, upper=pmax)
fit_std = stderror(fit)

p = plot(xlabel="Δ [rad/s]", ylabel="M / max(M)", xaxis=:log, legend=:none)
[scatter!(p, Δ, M[:,i], color=i) for i=1:length(ω1)]
[plot!(p, Δ, reshape(Henkelman_model(x, fit.param), length(Δ), length(ω1))[:,i], color=i) for i=1:length(ω1)]

This page was generated using Literate.jl.