Non-Linear Least Square Fitting

This section gives a brief overview of the interface to fit the generalized Bloch model to hybrid-state free precession data. We use the LsqFit.jl package and supply the algorithm with analytic gradients that are calculated with the calculatesignal_linearapprox function that implements the linear approximation of the generalized Bloch model for a train of rectangular RF pulses.

Basic Interface

This tutorial uses the following packages:

using MRIgeneralizedBloch
using MAT
using LinearAlgebra
using Plots

and we demonstrate the concept at the example of the RF pulse train that we published in the paper Rapid quantitative magnetization transfer imaging: utilizing the hybrid state and the generalized Bloch model:

control = matread(normpath(joinpath(pathof(MRIgeneralizedBloch), "../../docs/control_MT_v3p2_TR3p5ms_discretized.mat")))
α   = control["alpha"]
TRF = control["TRF"]

TR = 3.5e-3
t = TR .* (1:length(TRF));

As an example we can assume the following ground truth parameters

m0s = 0.15
R1f = 0.5 # 1/s
R2f = 17 # 1/s
Rx = 30 # 1/s
R1s = 3 # 1/s
T2s = 12e-6 # s
ω0 = 100 # rad/s
B1 = 0.9; # in units of B1_nominal

precompute the linear approximation

R2slT = precompute_R2sl();

and simulate the signal:

s = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT)
s = vec(s)

1142-element Vector{ComplexF64}:
    -0.3252173153487624 + 0.009389118419681007im
   -0.23144917764360645 + 0.004770811454990871im
   -0.12903774678498184 + 0.011349840151064092im
   -0.05678766513455557 + 0.0007970042626797638im
   -0.02675559139202654 + 0.00398275590019067im
  -0.016730390075552413 - 0.0022525327725481914im
  -0.012251559494301616 + 0.002743865569592166im
  -0.007219555756390814 - 0.001820502894584987im
  -0.009658488390964151 + 0.0012105896643954384im
 -0.0048500068341275114 - 0.0003889196669722614im
                        ⋮
   0.009356175012659677 + 0.00011234338454070591im
   0.010555785673980063 + 0.00020142661069044352im
   0.013338246077603656 + 0.0004095199938106825im
     0.0195565290569489 + 0.0008491700809056665im
    0.03288202710045477 + 0.0017547678348315706im
    0.06717246407976572 + 0.00474565268336824im
    0.15801530357449664 + 0.012316528152782029im
     0.2654078013422411 + 0.00913125790667977im
    0.35381344155560773 + 0.010214696894837762im

To make this example a bit more realistic, we add complex valued Gaussian noise:

s .+= 0.01 * randn(ComplexF64, size(s));

Now we can fit the model to the noisy data:

qM = fit_gBloch(s, α, TRF, TR; R2slT=R2slT)

MRIgeneralizedBloch.qMTparam(1.0000603362909752 - 0.0016297797749979179im, 0.14431223878563115, 0.28830663685586855, 15.874126170136815, 67.53915969012161, 4.070421488280747, 8.775154888316857e-6, 132.0549603612573, 0.8875832369391727, 0.3340982795273344, LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}}([1.0000603362909752, -0.0016297797749979179, 0.14431223878563115, 0.28830663685586855, 15.874126170136815, 67.53915969012161, 4.070421488280747, 8.775154888316857e-6, 132.0549603612573, 0.8875832369391727], [0.006340280067428261, 0.013034254889553593, 0.0016228534386894256, 0.003620649626918042, 0.0007116455690389413, 0.006541237752984989, -0.002345384396716478, -0.00036783917415302535, -0.006616730555333988, 0.00863286519863581  …  0.0038118764666522464, 0.008541565921625143, -0.007054058960158308, -0.009304693091517968, 0.003072296244234729, 0.004851555479870059, 0.0043876779702491184, 0.00915764145656504, -0.005324080630240749, 0.00442239526172844], [-0.32280463128340364 -0.011776153564477879 … -1.6113686618821224e-5 -0.2979406592774891; -0.22305683160930237 -0.008732995366004137 … -6.618388311559142e-5 -0.3973836958694736; … ; 0.012319812328615572 0.2649050640374011 … 0.00010886406098130218 0.0018428770580807027; 0.012845342671861396 0.35211294436639967 … 7.442192011068194e-5 -0.005123102548658863], false, Float64[]))

The last keyword argument is optional. It allows to recycle the precomputed R2sl, which improves speed. If not specified, it is re-calculated internally.

The results are stored in a struct and we can access the fitted parameters with

qM.m0s

0.14431223878563115

qM.R1f # 1/s

0.28830663685586855

qM.R2f # 1/s

15.874126170136815

qM.Rx # 1/s

67.53915969012161

qM.R1s # 1/s

4.070421488280747

1e6qM.T2s # μs

8.775154888316857

qM.ω0 # rad/s

132.0549603612573

qM.B1 # 1/B1_nominal

0.8875832369391727

We can also simulate the signal with the fitted parameters

s_fitted = calculatesignal_linearapprox(α, TRF, TR, qM.ω0, qM.B1, qM.m0s, qM.R1f, qM.R2f, qM.Rx, qM.R1s, qM.T2s, R2slT)
s_fitted = vec(s_fitted);

and compare it to the noisy data:

p = plot(xlabel="t (s)", ylabel="signal (normalized)", legend=:topleft)
plot!(p, t, real.(s), label="Re(s)")
plot!(p, t, imag.(s), label="Im(s)")
plot!(p, t, real.(s_fitted), label="Re(s_fitted)")
plot!(p, t, imag.(s_fitted), label="Im(s_fitted)")

Bounds and Fixed Parameters

Above example uses the default bounds

reM0 = (-Inf, 1, Inf)

imM0 = (-Inf, 0, Inf)

m0s = ( 0, 0.2, 1)

R1f = ( 0, 0.3, Inf)

R2f = ( 0, 15, Inf)

Rx = ( 0, 20, Inf)

R1s = ( 0, 3, Inf)

T2s = (8e-6,1e-5,12e-6)

ω0 = (-Inf, 0, Inf)

B1 = ( 0, 1, 1.5)

where the three entries refer to (minimum, start_value, maximum) (cf. fit_gBloch).

With keyword arguments, one can modify each of these bounds. For example:

qM = fit_gBloch(s, α, TRF, TR; R2slT=R2slT, m0s  = (0.1, 0.3, 0.5))

MRIgeneralizedBloch.qMTparam(1.000005006112668 - 0.0016172678478107248im, 0.14452050519863077, 0.28506533311857174, 15.85294308229628, 68.05410167568446, 4.082317259602439, 8.724773271705966e-6, 131.9817372448438, 0.8878780721404812, 0.3340983788961961, LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}}([1.000005006112668, -0.0016172678478107248, 0.14452050519863077, 0.28506533311857174, 15.85294308229628, 68.05410167568446, 4.082317259602439, 8.724773271705966e-6, 131.9817372448438, 0.8878780721404812], [0.006295251421042092, 0.013058163823236074, 0.0016370699688502133, 0.003639821991260707, 0.0007280629147803294, 0.0065364920242235725, -0.0023248695618289066, -0.00037953304253952706, -0.006600478931770599, 0.008627735882463516  …  0.003811847595387638, 0.0085415100829839, -0.007054183202997388, -0.009304976780414007, 0.0030716715279403475, 0.004850233741324036, 0.004384157263911317, 0.009147957760500213, -0.005333163743035605, 0.004417290901067331], [-0.32286736458018805 -0.011770878629227361 … -1.6109339703788192e-5 -0.2973651059073065; -0.22304518050368022 -0.008748615580181696 … -6.624598694055412e-5 -0.3966848662503803; … ; 0.012308076844468191 0.2648929357793602 … 0.00010879971756912456 0.001836806875388898; 0.012836517920540663 0.35209714091353966 … 7.442930610805302e-5 -0.005121713909451141], false, Float64[]))

starts the fit at m0s = 0.3 and uses a lower bound of 0.1 and an upper bound of 0.5. Alternatively, one also fix parameters to specified values:

qM = fit_gBloch(s, α, TRF, TR; R2slT=R2slT, ω0 = 0, B1 = 1)

MRIgeneralizedBloch.qMTparam(0.9612349967889942 + 0.004885749537742629im, 0.17459881134324323, 0.4430301625033134, 14.404162000539145, 37.630720871643916, 2.915192041087354, 8.0e-6, 0, 1, 0.3392428806228742, LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}}([0.9612349967889942, 0.004885749537742629, 0.17459881134324323, 0.4430301625033134, 14.404162000539145, 37.630720871643916, 2.915192041087354, 8.0e-6], [-0.00444278457875652, -0.0042960927359487755, -0.0034149417557093886, -0.006647279349637797, 0.0030327256214909767, 0.0011317602811418172, 0.00018196204079172062, -0.004260484414360185, -0.004080351312701211, 0.004879820208237875  …  0.003775299963735004, 0.008460675830904655, -0.007243883128235125, -0.009749289764075719, 0.0020913252600905136, 0.0027707058106057215, -0.0013717422375055372, -0.005979676299153148, -0.015865263868428185, -0.0060666282538192635], [-0.3470410270102214 2.456415114176918e-17 … -0.03765392337456717 2489.464089029737; -0.2500807875149091 1.705429078120691e-17 … -0.02718482189043301 1629.6009735723615; … ; 1.3940495556581356e-17 0.27582983876490774 … 0.00015229649026882817 -9.640063240681108; 1.8863540119732893e-17 0.3649854805289113 … 0.00020128278280038078 -13.307677251501708], false, Float64[]))

In this case, the derivatives wrt. ω0 and B1 are not calculated and the result is accordingly

qM.ω0

0

qM.B1

1

Linear Compression

As originally suggested by McGivney et al. for MR Fingerprinting, the manifold of signal evolution or fingerprints is low rank and it is often beneficial to reconstruct images directly in this domain. We can calculate a low rank basis with

sv = Array{ComplexF64}(undef, length(s), 50)
for i=1:size(sv,2)
    sv[:,i] = calculatesignal_linearapprox(α, TRF, TR, 500randn(), 0.8 + 0.4rand(), rand(), rand(), 20rand(), 30rand(), 3rand(), 8e-6+5e-6rand(), R2slT)
end
u, _, _ = svd(sv)
u = u[:,1:9];

where the rank 9 was chosen heuristically. The noisy signal can be compressed with

sc = u' * s

9-element Vector{ComplexF64}:
    2.1123828531977833 - 0.06033398193103418im
   0.27111782249324695 - 0.04414735125280569im
     0.375882769568456 - 0.01966546873826104im
 -0.030214071636922815 + 0.013886721680368091im
   0.12889999202146554 - 0.0031313488125832706im
 -0.010844539022797424 + 0.008567109651773im
  -0.07786380064406072 + 0.036639743935972434im
 -0.001666039620501986 - 0.035034109388607786im
 -0.017826804963322528 - 0.01131690634226613im

and fitted by calling fit_gBloch with the keyword argument u=u:

qM = fit_gBloch(sc, α, TRF, TR; R2slT=R2slT, u=u)

MRIgeneralizedBloch.qMTparam(1.0058434307845368 + 9.53310502979697e-5im, 0.14439481746152108, 0.3402991696778543, 16.646456711637455, 56.18796507762062, 3.692847079639061, 8.0e-6, 94.05429202084534, 0.8907566533528819, 0.030018548733129622, LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}}([1.0058434307845368, 9.53310502979697e-5, 0.14439481746152108, 0.3402991696778543, 16.646456711637455, 56.18796507762062, 3.692847079639061, 8.0e-6, 94.05429202084534, 0.8907566533528819], [0.0003533275832450222, -0.004909175631886531, -0.0018467844114915533, 0.007824092652957526, -0.00195956653871876, 0.0006444480395836217, 0.0002335804940261227, 0.0021122437993329538, 0.003122876809344477, -0.0024385839791879957, 0.022341799452417183, 0.0022065607188874077, -0.00028411352897537816, -0.011899398086542869, -0.009487499892802169, -0.000808384492127931, 0.0009726003092792659, 0.007104757913569293], [2.100456341147843 0.06260696515183412 … 0.00011919907070132382 0.16370702305874807; 0.26466005508705365 0.021703956553540675 … -0.00018260668220554212 -0.4026467743081026; … ; -0.03386367104549263 0.000440402458221174 … 2.4229649251361333e-5 0.2732250222138224; -0.0041862924829812465 -0.014618902691612881 … -1.5487830947028385e-6 0.06064221119558118], false, Float64[]))

Apparent $R_1$

Above fits tread R1f and R1s of the free and the semi-solid as independent parameters. As we discussed in our paper, many publications in the literature assume an apparent R1a = R1f = R1s. The corresponding model can be fitted by specifying fit_apparentR1=true:

R1a = 1 # 1/s
s = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1a, R2f, Rx, R1a, T2s, R2slT)
qM = fit_gBloch(vec(s), α, TRF, TR; fit_apparentR1=true, R1a = (0, 0.7, Inf), R2slT=R2slT)

MRIgeneralizedBloch.qMTparam(0.9999999999999625 - 3.5756561556125403e-17im, 0.15000000000000366, 1.0000000000000313, 16.999999999999964, 29.999999999998877, 1.0000000000000313, 1.1999999999998677e-5, 100.0000000000009, 0.9000000000000017, 3.98018078758594e-14, LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vector{Float64}}([0.9999999999999625, -3.5756561556125403e-17, 0.15000000000000366, 1.0000000000000313, 16.999999999999964, 29.999999999998877, 1.1999999999998677e-5, 100.0000000000009, 0.9000000000000017], [-1.2601031329495527e-14, -9.43689570931383e-15, -5.051514762044462e-15, -2.248201624865942e-15, -1.1483869410966463e-15, -4.753142324176451e-16, -6.38378239159465e-16, -9.8879238130678e-17, -5.065392549852277e-16, -7.719519468096792e-17  …  5.55653613398821e-18, 3.550762114890027e-18, 7.968885967768458e-18, 2.1846673775582914e-17, 3.2742905609062234e-17, 9.237402509576498e-17, 1.9255430583342559e-16, 5.585809592645319e-16, 3.8163916471489756e-16, 4.0766001685454967e-16], [-0.3660243243186932 -0.010586026795917107 … -9.513949151233961e-6 -0.33851676101566414; -0.2591710987881927 -0.005589974863073777 … -5.161513042198635e-5 -0.4144904053839076; … ; 0.01027201277001976 0.29861787950551993 … 0.00011260594964832189 0.0024409814452963593; 0.011516848569537719 0.39820858167216944 … 9.988694484535086e-5 -0.003324282121720797], true, Float64[]))

Here, we specified the limits of R1a explicitly, which is optional.

This page was generated using Literate.jl.