Skip to content

Commit ce92721

Browse files
mathieuboudreaumathieuboudreau
committed
Add MTsat
1 parent 4fe9f3a commit ce92721

11 files changed

+1126
-0
lines changed

‎.DS_Store‎

0 Bytes
Binary file not shown.

‎06-MT/03-MTsat/binder/postBuild‎

Whitespace-only changes.

‎06-MT/03-MTsat/binder/requirements.txt‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
jupyter-book==0.13.0
2+
lxml_html_clean
3+
nibabel
4+
pandas
5+
plotly
6+
scipy
7+
tabulate

‎06-MT/03-MTsat/binder/runtime.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
python-3.8

‎06-MT/03-MTsat/content/mtsat_fig_1_t1range_mtsat_3t_alternatesim.m‎

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
% NOTE TO SELF: Possibilty of simulation validation of MTsat through the
2+
% bloch simulator (i.e. know the true MTsat that happened during the
3+
% simulation, prior to fitting)
4+
clear all, close all, clc
5+
6+
%% Load protocol
7+
8+
fname = 'configs/mtsat-protocols.json';
9+
fid = fopen(fname);
10+
raw = fread(fid,inf);
11+
str = char(raw');
12+
fclose(fid);
13+
val = loadjson(str);
14+
15+
protocols = val.karakuzu2022.siemens1
16+
17+
%% Load tissues
18+
19+
fname = 'configs/tissues.json';
20+
fid = fopen(fname);
21+
raw = fread(fid,inf);
22+
str = char(raw');
23+
fclose(fid);
24+
val = loadjson(str);
25+
26+
tissue = val.sled2001.healthywhitematter;
27+
28+
%% T1 range
29+
30+
T1_true = 1
31+
T1_min = T1_true*0.7
32+
T1_max = T1_true*1.3
33+
34+
T1_range = linspace(T1_min, T1_max, 21)
35+
36+
37+
%%
38+
39+
MTsats = zeros(1,length(T1_range))
40+
MTRs = zeros(1,length(T1_range))
41+
T1s = zeros(1,length(T1_range))
42+
43+
for ii=1:length(T1_range)
44+
protocol = protocols.pdw
45+
46+
fa = protocol.fa
47+
tr = protocol.tr/1000
48+
te = protocol.te/1000
49+
offset = protocol.offset
50+
mt_shape = protocol.mtshape
51+
mt_duration = protocol.mtduration/1000
52+
mt_angle = protocol.mtangle
53+
54+
Model = qmt_spgr;
55+
Model.Prot.MTdata.Mat = [mt_angle, offset];
56+
Model.Prot.TimingTable.Mat(5) = tr ;
57+
Model.Prot.TimingTable.Mat(1) = mt_duration;
58+
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)) ;
59+
Model.options.Readpulsealpha = fa;
60+
Model.options.MT_Pulse_Shape = mt_shape
61+
62+
params = tissue{1}
63+
x = struct;
64+
x.F = params.F.mean;
65+
x.kr = params.kf.mean / x.F;
66+
x.R1f = 1/T1_range(ii);
67+
x.R1r = 1;
68+
x.T2f = params.T2f.mean/1000;
69+
x.T2r = params.T2r.mean/(10^6);
70+
71+
Opt.SNR = 1000;
72+
Opt.Method = 'Bloch sim';
73+
Opt.ResetMz = false;
74+
75+
[FitResult, MT_norm, PDw] = Model.Sim_Single_Voxel_Curve(x,Opt);
76+
77+
78+
79+
%% Get PDw/T1w ratio from analytical
80+
81+
PDw_Model = vfa_t1;
82+
83+
params.EXC_FA = 6;
84+
params.T1 = T1_range(ii); % Could improve by caclulating T1meas from qMT values
85+
params.TR = 0.032; % ms
86+
87+
PDw_anal = vfa_t1.analytical_solution(params);
88+
89+
T1w_Model = vfa_t1;
90+
91+
paramsT1w.EXC_FA = 20;
92+
paramsT1w.T1 = T1_range(ii); % ms
93+
paramsT1w.TR = 0.018; % ms
94+
95+
T1w_anal = vfa_t1.analytical_solution(paramsT1w);
96+
97+
PDwT1w_ratio = PDw_anal/T1w_anal
98+
%%
99+
100+
Model = mt_sat;
101+
FlipAngle = 6;
102+
TR = 0.032;
103+
Model.Prot.MTw.Mat = [ FlipAngle TR ];
104+
FlipAngle = 20;
105+
TR = 0.018;
106+
Model.Prot.T1w.Mat = [ FlipAngle TR];
107+
FlipAngle = 6;
108+
TR = 0.032;
109+
Model.Prot.PDw.Mat = [ FlipAngle TR];
110+
111+
data = struct();
112+
data.MTw=MT_norm*PDw;
113+
data.T1w=PDw/PDwT1w_ratio;
114+
data.PDw=PDw;
115+
FitResults = FitData(data,Model,0);
116+
MTsats(1,ii) = FitResults.MTSAT
117+
MTRs(1,ii) = FitResults.MTR
118+
T1s(1,ii) = FitResults.T1
119+
120+
end
121+
122+
%%
123+
close all
124+
125+
figure(1)
126+
plot(squeeze(T1_range), squeeze(MTRs),'LineWidth', 5)
127+
ylabel('MTR')
128+
legend('Location','northoutside')
129+
structHandler.figure = figure(1);
130+
structHandler.xlabel = xlabel('T1 (s)');
131+
structHandler.ylabel = ylabel('MTR');
132+
structHandler.legend = legend('Karakuzu2022 (Siemens 1)');
133+
figureProperties_plot(structHandler)
134+
135+
figure(2)
136+
plot(squeeze(T1_range), squeeze(MTsats),'LineWidth', 5)
137+
legend('Location','northoutside')
138+
structHandler.figure = figure(2);
139+
structHandler.xlabel = xlabel('T1 (s)');
140+
structHandler.ylabel = ylabel('MTsat');
141+
structHandler.legend = legend('Karakuzu2022 (Siemens 1)');
142+
figureProperties_plot(structHandler)
143+
144+
figure(3)
145+
plot(squeeze(T1_range), squeeze(T1s),'LineWidth', 5)
146+
legend('Location','northoutside')
147+
structHandler.figure = figure(3);
148+
structHandler.xlabel = xlabel('T1 (s)');
149+
structHandler.ylabel = ylabel('T1');
150+
structHandler.legend = legend('Karakuzu2022 (Siemens 1)');
151+
figureProperties_plot(structHandler)
152+
153+
save("fig1_mtsat.mat", "T1_range", "MTRs", "MTsats", "T1s")

‎06-MT/03-MTsat/content/mtsat_fig_2_b1range_mtsat_3t_alternatesim.m‎

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
% NOTE TO SELF: Possibilty of simulation validation of MTsat through the
2+
% bloch simulator (i.e. know the true MTsat that happened during the
3+
% simulation, prior to fitting)
4+
clear all, close all, clc
5+
6+
%% Load protocol
7+
8+
fname = 'configs/mtsat-protocols.json';
9+
fid = fopen(fname);
10+
raw = fread(fid,inf);
11+
str = char(raw');
12+
fclose(fid);
13+
val = loadjson(str);
14+
15+
protocols = val.karakuzu2022.siemens1
16+
17+
%% Load tissues
18+
19+
fname = 'configs/tissues.json';
20+
fid = fopen(fname);
21+
raw = fread(fid,inf);
22+
str = char(raw');
23+
fclose(fid);
24+
val = loadjson(str);
25+
26+
tissue = val.sled2001.healthywhitematter;
27+
28+
%% B1 range
29+
30+
B1_true = 1
31+
B1_min = B1_true*0.7
32+
B1_max = B1_true*1.3
33+
34+
B1_range = linspace(B1_min, B1_max, 21)
35+
36+
37+
%%
38+
39+
MTsats = zeros(1,length(B1_range))
40+
MTRs = zeros(1,length(B1_range))
41+
T1s = zeros(1,length(B1_range))
42+
43+
for ii=1:length(B1_range)
44+
protocol = protocols.pdw
45+
46+
fa = protocol.fa*B1_range(ii)
47+
tr = protocol.tr/1000
48+
te = protocol.te/1000
49+
offset = protocol.offset
50+
mt_shape = protocol.mtshape
51+
mt_duration = protocol.mtduration/1000
52+
mt_angle = protocol.mtangle*B1_range(ii)
53+
54+
Model = qmt_spgr;
55+
Model.Prot.MTdata.Mat = [mt_angle, offset];
56+
Model.Prot.TimingTable.Mat(5) = tr ;
57+
Model.Prot.TimingTable.Mat(1) = mt_duration;
58+
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)) ;
59+
Model.options.Readpulsealpha = fa;
60+
Model.options.MT_Pulse_Shape = mt_shape
61+
62+
params = tissue{1}
63+
x = struct;
64+
x.F = params.F.mean;
65+
x.kr = params.kf.mean / x.F;
66+
x.R1f = params.R1f.mean;
67+
x.R1r = 1;
68+
x.T2f = params.T2f.mean/1000;
69+
x.T2r = params.T2r.mean/(10^6);
70+
71+
Opt.SNR = 1000;
72+
Opt.Method = 'Bloch sim';
73+
Opt.ResetMz = false;
74+
75+
[FitResult, MT_norm, PDw] = Model.Sim_Single_Voxel_Curve(x,Opt);
76+
77+
78+
79+
%% Get PDw/T1w ratio from analytical
80+
81+
PDw_Model = vfa_t1;
82+
83+
params.EXC_FA = 6*B1_range(ii);
84+
params.T1 = 1/params.R1f.mean; % Could improve by caclulating T1meas from qMT values
85+
params.TR = 0.032; % ms
86+
87+
PDw_anal = vfa_t1.analytical_solution(params);
88+
89+
T1w_Model = vfa_t1;
90+
91+
paramsT1w.EXC_FA = 20*B1_range(ii);
92+
paramsT1w.T1 = 1/params.R1f.mean; % ms
93+
paramsT1w.TR = 0.018; % ms
94+
95+
T1w_anal = vfa_t1.analytical_solution(paramsT1w);
96+
97+
PDwT1w_ratio = PDw_anal/T1w_anal
98+
%%
99+
100+
Model = mt_sat;
101+
FlipAngle = 6; % These are the nominal flip angles used for fitting, so not corrected in these sims
102+
TR = 0.032;
103+
Model.Prot.MTw.Mat = [ FlipAngle TR ];
104+
FlipAngle = 20;
105+
TR = 0.018;
106+
Model.Prot.T1w.Mat = [ FlipAngle TR];
107+
FlipAngle = 6;
108+
TR = 0.032;
109+
Model.Prot.PDw.Mat = [ FlipAngle TR];
110+
111+
data = struct();
112+
data.MTw=MT_norm*PDw;
113+
data.T1w=PDw/PDwT1w_ratio;
114+
data.PDw=PDw;
115+
116+
FitResults = FitData(data,Model,0);
117+
MTsats(1,ii) = FitResults.MTSAT
118+
MTRs(1,ii) = FitResults.MTR
119+
T1s(1,ii) = FitResults.T1
120+
121+
end
122+
123+
%%
124+
close all
125+
126+
figure(1)
127+
plot(squeeze(B1_range), squeeze(MTRs), 'LineWidth', 5)
128+
legend('Location','northoutside')
129+
structHandler.figure = figure(1);
130+
structHandler.xlabel = xlabel('B1 (n.u.)');
131+
structHandler.ylabel = ylabel('MTR');
132+
structHandler.legend = legend('Karakuzu2022 (Siemens 1)');
133+
figureProperties_plot(structHandler)
134+
135+
figure(2)
136+
plot(squeeze(B1_range), squeeze(MTsats), 'LineWidth', 5)
137+
legend('Location','northoutside')
138+
structHandler.figure = figure(2);
139+
structHandler.xlabel = xlabel('B1 (n.u.)');
140+
structHandler.ylabel = ylabel('MTsat');
141+
structHandler.legend = legend('Karakuzu2022 (Siemens 1)');
142+
figureProperties_plot(structHandler)
143+
144+
figure(3)
145+
plot(squeeze(B1_range), squeeze(T1s), 'LineWidth', 5)
146+
legend('Location','northoutside')
147+
structHandler.figure = figure(3);
148+
structHandler.xlabel = xlabel('B1 (n.u.)');
149+
structHandler.ylabel = ylabel('T1');
150+
structHandler.legend = legend('Karakuzu2022 (Siemens 1)');
151+
figureProperties_plot(structHandler)
152+
153+
save("fig2a_mtsat.mat", "B1_range", "MTRs", "MTsats", "T1s")

0 commit comments

Comments
 (0)