From fce24ba502490d0f6052a9e4fd9490e9e946e568 Mon Sep 17 00:00:00 2001
From: Mathieu Boudreau <emb6150@gmail.com>
Date: Sat, 5 Oct 2024 23:10:54 -0300
Subject: [PATCH] Add Appendix B code/data

---
 .DS_Store                                     | Bin 10244 -> 10244 bytes
 06-MT/03-MTsat/binder/postBuild               |   2 +-
 .../content/MTsat_simulation_scripts.m        |  69 ++++++++++
 .../content/betterMTsatsimulationscripts.m    |  71 ++++++++++
 06-MT/03-MTsat/content/blog_simulation_1.m    | 128 ++++++++++++++++++
 .../content/blog_simulation_1_wmlesion.m      | 128 ++++++++++++++++++
 06-MT/03-MTsat/content/blog_simulation_2.m    |  27 ++++
 06-MT/03-MTsat/content/blog_simulation_3.m    |  20 +++
 .../content/blog_simulation_4_final.m         |  19 +++
 9 files changed, 463 insertions(+), 1 deletion(-)
 create mode 100644 06-MT/03-MTsat/content/MTsat_simulation_scripts.m
 create mode 100644 06-MT/03-MTsat/content/betterMTsatsimulationscripts.m
 create mode 100644 06-MT/03-MTsat/content/blog_simulation_1.m
 create mode 100644 06-MT/03-MTsat/content/blog_simulation_1_wmlesion.m
 create mode 100644 06-MT/03-MTsat/content/blog_simulation_2.m
 create mode 100644 06-MT/03-MTsat/content/blog_simulation_3.m
 create mode 100644 06-MT/03-MTsat/content/blog_simulation_4_final.m

diff --git a/.DS_Store b/.DS_Store
index d11e5c69beb320010550958ecb9974f9d20e8346..89750ac2ff7462623384ce058aaf73bd21a86e74 100644
GIT binary patch
delta 179
zcmZn(XbG6$&nUMsU^hRb++-dBso2<Q3=9k`40;Ud44Diix%n<GNjdpRKrxOVj@;24
zk;fcS<x}v=7i1U)C+FuD0QE31x;&WNEbx+%VRNovGKX49adJ*letr%kP>PkofWeGG
hm%*1IBqtq7hfMF}`;rMvasitirPi}<W>@&j4gk|DEr|dC

delta 149
zcmZn(XbG6$&nUYwU^hRb>|`DRsaP)?1_lNe20ey!hD?T%+<X_8q@4UDpcuz)u5x?k
z$YYME@+o-b3o;CYlk;;6fO;4hT>>UI3%q1x+?*?z%rW`7Bp)Nw<ZqG*OtKD}Bc#@|
NPAm}L%&zd49RR<xCa?ei

diff --git a/06-MT/03-MTsat/binder/postBuild b/06-MT/03-MTsat/binder/postBuild
index 7fd3d49..92c771d 100644
--- a/06-MT/03-MTsat/binder/postBuild
+++ b/06-MT/03-MTsat/binder/postBuild
@@ -3,5 +3,5 @@
 
 git clone https://www.github.com/qMRLab/qMRLab.git ../qMRLab
 cd ../qMRLab
-git checkout 5c64c747a7c08102201bb1404c868d6caf51d08d
+git checkout be314a4b44ca655e409a151bd8ee69a27a0d1ce0 
 cd binder
diff --git a/06-MT/03-MTsat/content/MTsat_simulation_scripts.m b/06-MT/03-MTsat/content/MTsat_simulation_scripts.m
new file mode 100644
index 0000000..e87ccb8
--- /dev/null
+++ b/06-MT/03-MTsat/content/MTsat_simulation_scripts.m
@@ -0,0 +1,69 @@
+
+T1_true = 1
+
+
+
+%% Get PDw/T1w ratio from analytical
+
+PDw_Model = vfa_t1;
+
+params.EXC_FA = 7;
+params.T1 = T1_true; % Could improve by caclulating T1meas from qMT values
+params.TR = 0.025; % ms
+
+PDw_anal = vfa_t1.analytical_solution(params);
+
+T1w_Model = vfa_t1;
+
+paramsT1w.EXC_FA = 20;
+paramsT1w.T1 = T1_true; % ms
+paramsT1w.TR = 0.018; % ms
+
+T1w_anal = vfa_t1.analytical_solution(paramsT1w);
+
+PDwT1w_ratio = PDw_anal/T1w_anal
+%%
+
+Model = mt_sat;
+FlipAngle = 7;
+TR  = 0.025;
+Model.Prot.MTw.Mat = [ FlipAngle TR ];
+FlipAngle = 20;
+TR = 0.018;
+Model.Prot.T1w.Mat = [ FlipAngle TR];
+FlipAngle = 7;
+TR = 0.025;
+Model.Prot.PDw.Mat = [ FlipAngle TR];
+
+data = struct();
+%data.MTw=MT_norm*PDw;
+data.MTw=MTw
+data.T1w=PDw/PDwT1w_ratio;
+data.PDw=PDw;
+FitResults = FitData(data,Model,0);
+MTsats = FitResults.MTSAT
+MTRs = FitResults.MTR
+T1s = FitResults.T1
+
+%%
+
+
+Model = mt_sat;
+FlipAngle = 7;
+TR  = 0.025;
+Model.Prot.MTw.Mat = [ FlipAngle TR ];
+FlipAngle = 20;
+TR = 0.018;
+Model.Prot.T1w.Mat = [ FlipAngle TR];
+FlipAngle = 7;
+TR = 0.025;
+Model.Prot.PDw.Mat = [ FlipAngle TR];
+
+data = struct();
+data.MTw=MTw;
+data.T1w=T1w;
+data.PDw=PDw;
+FitResults = FitData(data,Model,0);
+MTsats = FitResults.MTSAT
+MTRs = FitResults.MTR
+T1s = FitResults.T1
diff --git a/06-MT/03-MTsat/content/betterMTsatsimulationscripts.m b/06-MT/03-MTsat/content/betterMTsatsimulationscripts.m
new file mode 100644
index 0000000..b77591a
--- /dev/null
+++ b/06-MT/03-MTsat/content/betterMTsatsimulationscripts.m
@@ -0,0 +1,71 @@
+
+T1_true = 1
+
+
+
+%% Get PDw/T1w ratio from analytical
+
+PDw_Model = vfa_t1;
+
+params.EXC_FA = 7;
+params.T1 = T1_true; % Could improve by caclulating T1meas from qMT values
+params.TR = 0.025; % ms
+
+PDw_anal = vfa_t1.analytical_solution(params);
+
+T1w_Model = vfa_t1;
+
+paramsT1w.EXC_FA = 20;
+paramsT1w.T1 = T1_true; % ms
+paramsT1w.TR = 0.018; % ms
+
+T1w_anal = vfa_t1.analytical_solution(paramsT1w);
+
+PDwT1w_ratio = (PDw_anal/sind(params.EXC_FA))/(T1w_anal/sind(paramsT1w.EXC_FA))
+%%
+
+Model = mt_sat;
+FlipAngle = 7;
+TR  = 0.025;
+Model.Prot.MTw.Mat = [ FlipAngle TR ];
+FlipAngle = 20;
+TR = 0.018;
+Model.Prot.T1w.Mat = [ FlipAngle TR];
+FlipAngle = 7;
+TR = 0.025;
+Model.Prot.PDw.Mat = [ FlipAngle TR];
+
+data = struct();
+%data.MTw=MT_norm*PDw;
+data.MTw=MTw*cosd(params.EXC_FA)
+data.T1w=PDw/PDwT1w_ratio*cosd(paramsT1w.EXC_FA);
+data.PDw=PDw*cosd(params.EXC_FA);
+FitResults = FitData(data,Model,0);
+MTsats = FitResults.MTSAT
+MTRs = FitResults.MTR
+T1s = FitResults.T1
+
+%%
+
+%%
+
+
+Model = mt_sat;
+FlipAngle = 7;
+TR  = 0.025;
+Model.Prot.MTw.Mat = [ FlipAngle TR ];
+FlipAngle = 20;
+TR = 0.018;
+Model.Prot.T1w.Mat = [ FlipAngle TR];
+FlipAngle = 7;
+TR = 0.025;
+Model.Prot.PDw.Mat = [ FlipAngle TR];
+
+data = struct();
+data.MTw=MTw*cosd(params.EXC_FA);
+data.T1w=T1w*cosd(paramsT1w.EXC_FA);
+data.PDw=PDw*cosd(params.EXC_FA);
+FitResults = FitData(data,Model,0);
+MTsats = FitResults.MTSAT
+MTRs = FitResults.MTR
+T1s = FitResults.T1
\ No newline at end of file
diff --git a/06-MT/03-MTsat/content/blog_simulation_1.m b/06-MT/03-MTsat/content/blog_simulation_1.m
new file mode 100644
index 0000000..371bf27
--- /dev/null
+++ b/06-MT/03-MTsat/content/blog_simulation_1.m
@@ -0,0 +1,128 @@
+clear all, close all, clc
+
+%% Load protocol
+
+fname = 'configs/mtsat-protocols.json';
+fid = fopen(fname);
+raw = fread(fid,inf);
+str = char(raw');
+fclose(fid);
+val = loadjson(str);
+
+protocols = val.karakuzu2022.siemens1
+
+%% Load tissues
+
+fname = 'configs/tissues.json';
+fid = fopen(fname);
+raw = fread(fid,inf);
+str = char(raw');
+fclose(fid);
+val = loadjson(str);
+
+tissue = val.sled2001.healthywhitematter;
+
+%% T1 range
+
+T1_true = 1
+
+%% qMT-SPGR experiment
+
+MTsats = zeros(1,length(T1_true))
+MTRs = zeros(1,length(T1_true))
+T1s = zeros(1,length(T1_true))
+
+
+protocol = protocols.pdw
+
+fa = protocol.fa
+tr = protocol.tr/1000
+te = protocol.te/1000
+offset = protocol.offset
+mt_shape = protocol.mtshape
+mt_duration = protocol.mtduration/1000
+mt_angle = protocol.mtangle
+
+Model = qmt_spgr;
+Model.Prot.MTdata.Mat = [mt_angle, offset];
+Model.Prot.TimingTable.Mat(5) = tr ;
+Model.Prot.TimingTable.Mat(1) = mt_duration;
+Model.Prot.TimingTable.Mat(4) = Model.Prot.TimingTable.Mat(5) - (Model.Prot.TimingTable.Mat(1) + Model.Prot.TimingTable.Mat(2) + Model.Prot.TimingTable.Mat(3)) ;
+Model.options.Readpulsealpha = fa;
+Model.options.MT_Pulse_Shape = mt_shape
+
+params = tissue{1}
+x = struct;
+x.F = params.F.mean;
+x.kr = params.kf.mean / x.F;
+x.R1f = 1/T1_true;
+x.R1r = 1;
+x.T2f = params.T2f.mean/1000;
+x.T2r = params.T2r.mean/(10^6);
+
+Opt.SNR = 1000;
+Opt.Method = 'Bloch sim';
+Opt.ResetMz = false;
+
+[FitResult, Smodel] = Model.Sim_Single_Voxel_Curve(x,Opt); % NOTE: this uses a modified version of the qmt_spgr.m file where the additional output is included. Not all version of qMRLab has this; if yours doesn't, go to the file and add the additional function output accordingly.
+
+%% Cleanup
+
+%Smodel is the normalized MT-SPGR value, that is, the signal with the MT
+%pulse on divided by the signal from the same sequence with the MT pulse
+%off. Since in the MTsat experiment, the PD-weighted pulse sequence is the
+%latter case above, we then define:
+
+Signal_MT = Smodel;
+Signal_PDw = 1;
+
+%% Find scaling value for T1w signal
+
+
+% Get PDw/T1w ratio from analytical
+
+PDw_Model = vfa_t1;
+
+params.EXC_FA = protocols.pdw.fa;
+params.T1 = T1_true; % Could improve by caclulating T1meas from qMT values
+params.TR = protocols.pdw.tr/1000; % ms
+
+PDw_anal = vfa_t1.analytical_solution(params);
+
+T1w_Model = vfa_t1;
+
+paramsT1w.EXC_FA = protocols.t1w.fa;
+paramsT1w.T1 = T1_true; % ms
+paramsT1w.TR = protocols.t1w.tr/1000; % ms
+
+T1w_anal = vfa_t1.analytical_solution(paramsT1w);
+
+T1wPDw_ratio = T1w_anal/PDw_anal
+
+%% Cleanup
+
+% Since Signal_PDw = 1, then it's clear that 
+
+Signal_T1w = T1wPDw_ratio
+
+%% Calculate MTsat from signals
+
+Model = mt_sat;
+FlipAngle =  protocols.pdw.fa;
+TR =  protocols.pdw.tr/1000;
+Model.Prot.MTw.Mat = [ FlipAngle TR ];
+FlipAngle = protocols.t1w.fa;
+TR =  protocols.t1w.tr/1000;
+Model.Prot.T1w.Mat = [ FlipAngle TR];
+FlipAngle =  protocols.pdw.fa;
+TR =  protocols.pdw.tr/1000;
+Model.Prot.PDw.Mat = [ FlipAngle TR];
+
+data = struct();
+data.MTw=Signal_MT;
+data.T1w=Signal_T1w;
+data.PDw=Signal_PDw;
+FitResults = FitData(data,Model,0);
+MTsats = FitResults.MTSAT
+MTRs = FitResults.MTR
+T1s = FitResults.T1
\ No newline at end of file
diff --git a/06-MT/03-MTsat/content/blog_simulation_1_wmlesion.m b/06-MT/03-MTsat/content/blog_simulation_1_wmlesion.m
new file mode 100644
index 0000000..b7c7e37
--- /dev/null
+++ b/06-MT/03-MTsat/content/blog_simulation_1_wmlesion.m
@@ -0,0 +1,128 @@
+clear all, close all, clc
+
+%% Load protocol
+
+fname = 'configs/mtsat-protocols.json';
+fid = fopen(fname);
+raw = fread(fid,inf);
+str = char(raw');
+fclose(fid);
+val = loadjson(str);
+
+protocols = val.karakuzu2022.siemens1
+
+%% Load tissues
+
+fname = 'configs/tissues.json';
+fid = fopen(fname);
+raw = fread(fid,inf);
+str = char(raw');
+fclose(fid);
+val = loadjson(str);
+
+tissue = val.sled2001.latewmlesion;
+
+%% T1 range
+
+T1_true = 1
+
+%% qMT-SPGR experiment
+
+MTsats = zeros(1,length(T1_true))
+MTRs = zeros(1,length(T1_true))
+T1s = zeros(1,length(T1_true))
+
+
+protocol = protocols.pdw
+
+fa = protocol.fa
+tr = protocol.tr/1000
+te = protocol.te/1000
+offset = protocol.offset
+mt_shape = protocol.mtshape
+mt_duration = protocol.mtduration/1000
+mt_angle = protocol.mtangle
+
+Model = qmt_spgr;
+Model.Prot.MTdata.Mat = [mt_angle, offset];
+Model.Prot.TimingTable.Mat(5) = tr ;
+Model.Prot.TimingTable.Mat(1) = mt_duration;
+Model.Prot.TimingTable.Mat(4) = Model.Prot.TimingTable.Mat(5) - (Model.Prot.TimingTable.Mat(1) + Model.Prot.TimingTable.Mat(2) + Model.Prot.TimingTable.Mat(3)) ;
+Model.options.Readpulsealpha = fa;
+Model.options.MT_Pulse_Shape = mt_shape
+
+params = tissue{1}
+x = struct;
+x.F = params.F.mean;
+x.kr = params.kf.mean / x.F;
+x.R1f = 1/T1_true;
+x.R1r = 1;
+x.T2f = params.T2f.mean/1000;
+x.T2r = params.T2r.mean/(10^6);
+
+Opt.SNR = 1000;
+Opt.Method = 'Bloch sim';
+Opt.ResetMz = false;
+
+[FitResult, Smodel] = Model.Sim_Single_Voxel_Curve(x,Opt); % NOTE: this uses a modified version of the qmt_spgr.m file where the additional output is included. Not all version of qMRLab has this; if yours doesn't, go to the file and add the additional function output accordingly.
+
+%% Cleanup
+
+%Smodel is the normalized MT-SPGR value, that is, the signal with the MT
+%pulse on divided by the signal from the same sequence with the MT pulse
+%off. Since in the MTsat experiment, the PD-weighted pulse sequence is the
+%latter case above, we then define:
+
+Signal_MT = Smodel;
+Signal_PDw = 1;
+
+%% Find scaling value for T1w signal
+
+
+% Get PDw/T1w ratio from analytical
+
+PDw_Model = vfa_t1;
+
+params.EXC_FA = protocols.pdw.fa;
+params.T1 = T1_true; % Could improve by caclulating T1meas from qMT values
+params.TR = protocols.pdw.tr/1000; % ms
+
+PDw_anal = vfa_t1.analytical_solution(params);
+
+T1w_Model = vfa_t1;
+
+paramsT1w.EXC_FA = protocols.t1w.fa;
+paramsT1w.T1 = T1_true; % ms
+paramsT1w.TR = protocols.t1w.tr/1000; % ms
+
+T1w_anal = vfa_t1.analytical_solution(paramsT1w);
+
+T1wPDw_ratio = T1w_anal/PDw_anal
+
+%% Cleanup
+
+% Since Signal_PDw = 1, then it's clear that 
+
+Signal_T1w = T1wPDw_ratio
+
+%% Calculate MTsat from signals
+
+Model = mt_sat;
+FlipAngle =  protocols.pdw.fa;
+TR =  protocols.pdw.tr/1000;
+Model.Prot.MTw.Mat = [ FlipAngle TR ];
+FlipAngle = protocols.t1w.fa;
+TR =  protocols.t1w.tr/1000;
+Model.Prot.T1w.Mat = [ FlipAngle TR];
+FlipAngle =  protocols.pdw.fa;
+TR =  protocols.pdw.tr/1000;
+Model.Prot.PDw.Mat = [ FlipAngle TR];
+
+data = struct();
+data.MTw=Signal_MT;
+data.T1w=Signal_T1w;
+data.PDw=Signal_PDw;
+FitResults = FitData(data,Model,0);
+MTsats = FitResults.MTSAT
+MTRs = FitResults.MTR
+T1s = FitResults.T1
\ No newline at end of file
diff --git a/06-MT/03-MTsat/content/blog_simulation_2.m b/06-MT/03-MTsat/content/blog_simulation_2.m
new file mode 100644
index 0000000..1af96a2
--- /dev/null
+++ b/06-MT/03-MTsat/content/blog_simulation_2.m
@@ -0,0 +1,27 @@
+clear all, close all, clc
+
+%% Load Mz before and after for each TR
+
+load('sim2.mat')
+
+%% Plot Mz before and after for each TR
+
+figure()
+plot(Mz_before, 'r')
+hold on
+plot(Mz_after, 'b')
+legend('Mz_{before}', 'Mz_{after}')
+
+figure()
+plot(Mz_before(end-10+1:end), 'r')
+hold on
+plot(Mz_after(end-10+1:end), 'b')
+legend('Mz_{before}', 'Mz_{after}')
+
+figure()
+plot(1-Mz_after./Mz_before)
+legend('1-Mz_{after}/Mz_{before}')
+
+figure()
+plot((1-Mz_after./Mz_before)*100)
+legend('1-Mz_{after}/Mz_{before}')
diff --git a/06-MT/03-MTsat/content/blog_simulation_3.m b/06-MT/03-MTsat/content/blog_simulation_3.m
new file mode 100644
index 0000000..94dce16
--- /dev/null
+++ b/06-MT/03-MTsat/content/blog_simulation_3.m
@@ -0,0 +1,20 @@
+clear all, close all, clc
+
+%% Load Mz before and after for each TR
+
+load('sim2.mat')
+
+%% Plot Mz before and after for each TR
+
+figure()
+
+figure()
+plot((1-Mz_after./Mz_before)*100)
+hold on
+plot((1-(Mz_after-delta_Mz_T1relax)./Mz_before)*100)
+
+legend('MTsat before T1 correction', 'MTsat after T1 correction')
+
+figure()
+plot((1-(Mz_after-delta_Mz_T1relax)./Mz_before)*100)
+legend('1-Mz_{after}/Mz_{before}')
diff --git a/06-MT/03-MTsat/content/blog_simulation_4_final.m b/06-MT/03-MTsat/content/blog_simulation_4_final.m
new file mode 100644
index 0000000..df44313
--- /dev/null
+++ b/06-MT/03-MTsat/content/blog_simulation_4_final.m
@@ -0,0 +1,19 @@
+clear all, close all, clc
+
+%% Load Mz before and after for each TR
+
+load('sim2.mat')
+
+%% Plot Mz before and after for each TR
+
+figure()
+
+figure()
+plot((1-(Mz_after-delta_Mz_T1relax)./Mz_before)*100)
+
+hold on
+plot((1-(M0_remainingTR_free-delta_Mz_T1relax_remaining)./Mz_before)*100)
+plot((1-(Mz_after-delta_Mz_T1relax)./Mz_before)*100+(1-(M0_remainingTR_free-delta_Mz_T1relax_remaining)./Mz_before)*100)
+
+legend('MTsat contribution from MT pulse event', 'MTsat contribution from cross-relaxation event', 'Total MTsat for TR')
+