Skip to content

Commit

Permalink
some minor changes fpor better functionality and simpler
Browse files Browse the repository at this point in the history 
added BM simulation_paramspace
  • Loading branch information
mzaiss committed Apr 10, 2018
1 parent e827d33 commit 6bf803d
Show file tree
Hide file tree
Showing 4 changed files with 187 additions and 12 deletions.
4 changes: 3 additions & 1 deletion doc/2_BMsimulation.m
View file
Expand Up @@ -6,6 +6,7 @@


%%

Sim.n_cest_pool=1;
Sim.dwA=0;
Sim.R1A=1;
Expand Down Expand Up @@ -41,7 +42,8 @@
Sim.flipangle=5;
Sim.readout='FID';
Sim.spoilf=0;
Sim.TTM_rep =0;



num = NUMERIC_SIM(Sim);
figure(1), plot(num.x,num.zspec,'.'); hold on;
Expand Down
166 changes: 166 additions & 0 deletions doc/3_BMsimulation_param_space.m
View file
@@ -0,0 +1,166 @@
%%%%%%%%%%% BATCH_simulation %%%%%%%%%%%
% - run it section by section
% - results are NOT saved automatically
% - numeric (BM) solution is calculated in NUMERIC_SIM
% it supports 1 water (A), 1 MT (C) and up to 5 CEST/NOE pools (B,D,E,F,G)
% - analytic solution is calculated in ANALYTIC_SIM
% it supports 1 water (A), 1 MT (C) and up to 5 CEST/NOE pools (B,D,E,F,G)
% - always BOTH solutions are calculated and plotted

% last change: 2018/04/10

%% clear all
clear all
clc



%% BASIC SIMULATION PARAMETERS
Sim=init_Sim(struct());
Sim.analytic = 1; % calculate analtical solution? 1=yes, 0=no
Sim.numeric = 0; % calculate numerical solution? 1=yes, 0=no
Sim.all_offsets = 1; % 1 = z-spectrum, 0 = simulate only +- offset, 2 = on-resonant
Sim.offset = 70; % offset range in ppm
Sim.n_cest_pool = 5; % number of CEST/NOE pools (CEST pools: B,D,E,F,G...)


%% SEQUENCE AND SCANNER PARAMETERS


Sim.FREQ = 7*gamma_; % frequency (=B0[T] * gamma)
Sim.B1 = 5; % standard B1 value in µT
Sim.Trec = 2.5; % standard recover time in s
Sim.Zi = 0.1; % initial magnetisation (should be between -1 and +1)
Sim.shape = 'block'; % cases: SPINLOCK, seq_gauss, block, AdiaSL, AdiaSinCos, AdiaInversion,
% block_trap, gauss, sech, sinc_1, sinc_2, sinc_3, sinc_4
Sim.pulsed = 0; % 0 = cw saturation, 1 = pulsed saturation

% settings for pulsed saturation
if Sim.pulsed
Sim.n = 7; % number of saturation pulses
Sim.tp = 0.1; % saturation time per pulse in s
Sim.DC = 0.70; % duty cycle

% settings for continuous wave (cw) saturatation
else
Sim.tp = 3; % saturation time in s
Sim.n = 1; % choose n=1 for cw saturation
Sim.shape = 'block'; % choose 'block' for cw saturation
Sim.DC = 1.0; % choose DC=1 for cw saturation
end;


%% CHOOSE TISSUE-TYPE AND CEST-AGENT AND LOAD PARAMETERS WITH getSim

% Pool system parameters
% water pool A
Sim.dwA=0;
Sim.R1A=1/3; % PBS
Sim.R2A=1/2; % PBS
% Sim.R1A=1/1.820; % GM3T
% Sim.R2A=1/0.099; % GM3T


% first CEST pool B (paraCEST pool)
Sim.fB=0.00009; % rel. conc 10mM/111M
Sim.kBA=1000; % exchange rate in Hz ( the fast one, kBA is calculated by this and fB)
Sim.dwB=55; % ppm relative to dwA
Sim.R1B=1; % R1B relaxation rate [Hz]
Sim.R2B=50; % R2B relaxation rate [Hz]

% second CEST pool D (amide)
Sim.fD=0.00009; % rel. conc 10mM/111M
Sim.kDA=6000; % exchange rate in Hz ( the fast one, kBA is calculated by this and fB)
Sim.dwD=-5; % ppm relative to dwA
Sim.R1D=1; % R1B relaxation rate [Hz]
Sim.R2D=50; % R2B relaxation rate [Hz]

Sim.n_cest_pool=2;
Sim.MT = 0;


%MT
% Sim.MT = 1; % 1 = with MT pool (pool C), 0 = no MT pool
% Sim.MT_lineshape = 'SuperLorentzian'; % MT lineshape - cases: SuperLorentzian, Gaussian, Lorentzian
% Sim.R1C=1;
% Sim.fC=0.05;
% Sim.R2C=109890; % 1/9.1µs
% Sim.kCA=40; Sim.kAC=Sim.kCA*Sim.fC;

SimStart = Sim;

%% CHOOSE PARAMETER'S VALUES TO BE SIMULATED
% !!! dont delete parameters - just comment them out !!!
clear Space
% settings for all_offsets = 1 (complete Z-spectrum)
if Sim.all_offsets == 1

Space.tp = [0.5 1 2 3 4];
Space.B1 = [5 10 15 25 50 100];
Space.Trec = [0 1 2 5 10 20];


else % settings for all_offsets = 0 / 2

Space.tp = linspace(0.05,20,50);
Space.B1 = logspace(-1,2,50);
Space.Trec = linspace(0,20,50);

end;

names = fieldnames(Space);

for jj = 1:numel(names)
Sim = SimStart; % set back to standard Parameter
field = names{jj};

for ii = 1:numel(Space.(field))

% do NOT change
Sim.(field) = Space.(field)(ii);
Sim.kAB = Sim.kBA*Sim.fB;
Sim.kAC = Sim.kCA*Sim.fC;
Sim.kAD = Sim.kDA*Sim.fD;
Sim.kAE = Sim.kEA*Sim.fE;
Sim.kAF = Sim.kFA*Sim.fF;
Sim.kAG = Sim.kGA*Sim.fG;
Sim.td = calc_td(Sim.tp,Sim.DC);
Sim.Zi= 1- (1-SimStart.Zi)*exp(-Sim.R1A*Sim.Trec);

% adjust your spectral resolution here
if Sim.all_offsets == 1
Sim.xZspec = [Sim.dwA-Sim.offset:0.5:Sim.dwA+Sim.offset];
Sim.xxZspec = [Sim.dwA-Sim.offset:0.5:Sim.dwA+Sim.offset];

elseif Sim.all_offsets == 0 % !!! CHANGED FROM +- dwB to +- Offset !!!
Sim.xZspec = [Sim.offset -Sim.offset];
Sim.xxZspec = [Sim.offset -Sim.offset];

elseif Sim.all_offsets == 2
Sim.xZspec = 0;
Sim.xxZspec = 0;
end;

% do NOT change
if Sim.numeric && Sim.analytic
NUMERIC_SPACE.(field){ii} = NUMERIC_SIM(Sim);
ANALYTIC_SPACE.(field){ii} = ANALYTIC_SIM(Sim);
elseif Sim.numeric && ~Sim.analytic
NUMERIC_SPACE.(field){ii} = NUMERIC_SIM(Sim);
ANALYTIC_SPACE = NUMERIC_SPACE;
elseif ~Sim.numeric && Sim.analytic
ANALYTIC_SPACE.(field){ii} = ANALYTIC_SIM(Sim);
NUMERIC_SPACE = ANALYTIC_SPACE;
end
end
end

clear ii jj names field


%% PLOT SIMULATED ASYM- AND Z-SPECTRA
% - numerical solutions are displayed as diamonds
% - analytical solutions are displayed as solid lines
Sim.modelfield = 'zspec'; %standard value = 'zspec'
PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE);

24 changes: 14 additions & 10 deletions general functions/PLOT_SPACE.m
View file
Expand Up @@ -6,12 +6,12 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE)
assert( ~isempty( fieldnames(NUMERIC_SPACE) ) || ~isempty( fieldnames(ANALYTIC_SPACE) ), 'No data to plot given.' );

% adjust param_xscale for all_offsets = 0 plots
% param_xscale={'lin','lin','lin','lin','log','log','lin','log','log'};
param_xscale={'lin','log','log','log'};
param_xscale={'lin','lin','lin','lin','lin','lin','lin','lin','lin'};


% set FontSize and LineWidth for legend/label/title...
myLineWidth = 1;
myMarkerSize = 2;
myMarkerSize = 6;
fontsize_legend = 12;
fontsize_label = 14;
fontsize_ticks = 14;
Expand Down Expand Up @@ -117,7 +117,7 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE)
if Sim.numeric && Sim.analytic
% numerical solution
% plot(xsim,Zsim,'Marker','d','Markersize', 5,'MarkerFaceColor','black','color',cc(ii,:),'LineStyle','none');
plot(xsim,Zsim,'Marker','o','Markersize', myMarkerSize,'MarkerFaceColor',cc(ii,:),'color',cc(ii,:),'LineStyle','none');
plot(xsim,Zsim,'Marker','.','Markersize', myMarkerSize,'MarkerFaceColor',cc(ii,:),'color',cc(ii,:),'LineStyle','none');
hold on;
% analytical solution
h = plot(xmod,Zmod,'LineWidth',myLineWidth,'color',cc(ii,:),'DisplayName',str);
Expand All @@ -127,7 +127,7 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE)
% plot numerical solution only
elseif Sim.numeric && ~Sim.analytic
% h = plot(xsim,abs(Zsim),'Marker','d','Markersize', 5,'MarkerFaceColor','black','color',cc(ii,:),'LineStyle','-','DisplayName',str);
h = plot(xsim,abs(Zsim),'Marker','o','Markersize', myMarkerSize,'MarkerFaceColor',cc(ii,:),'color',cc(ii,:),'LineStyle','-','LineWidth',myLineWidth,'DisplayName',str);
h = plot(xsim,abs(Zsim),'Marker','.','Markersize', myMarkerSize,'MarkerFaceColor',cc(ii,:),'color',cc(ii,:),'LineStyle','-','LineWidth',myLineWidth,'DisplayName',str);
hold on;
hvec(ii)= h; % handle vector for the legend
% title({sprintf('%s' ,PLabel.(field));'\fontsize{8} numerical (\diamondsuit) solutions'},'FontSize',12);
Expand Down Expand Up @@ -195,6 +195,7 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE)




% % % % % % % % %
% % % AREX % % %
% % % % % % % % %
Expand Down Expand Up @@ -222,7 +223,7 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE)
% % plot numerical solution only
% elseif Sim.numeric && ~Sim.analytic
% % h2 = plot(x_arex,y_arex,'Marker','d','Markersize', 5,'MarkerFaceColor','black','color',cc(ii,:),'LineStyle','none','DisplayName',str);
% h2 = plot(x_arex,y_arex,'Marker','o','Markersize', 5,'color',cc(ii,:),'LineStyle','-','LineWidth',1.5,'DisplayName',str);
% h2 = plot(x_arex,y_arex,'Marker','.','Markersize', 5,'color',cc(ii,:),'LineStyle','-','LineWidth',1.5,'DisplayName',str);
% hold on;
% hvec2(ii)= h2; % handle vector for the legend
% % title({sprintf('%s' ,PLabel.(field));'\fontsize{8} numerical (\diamondsuit) solutions'},'FontSize',12);
Expand Down Expand Up @@ -252,7 +253,7 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE)
% plot analytical solution only
h2 = plot(xsim,R1rho,'LineWidth',myLineWidth,'color',cc(ii,:),'DisplayName',str);
hold on;
hvec2(ii)= h2; % handle vector for the legend
hvec3(ii)= h2; % handle vector for the legend
% title({sprintf('%s' ,PLabel.(field));'\fontsize{8} analytical (---) solutions'},'FontSize',12);

% general plot settings
Expand Down Expand Up @@ -331,10 +332,13 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE)
% plot Z(offset) over varyvals
figure(333),
subplot(plot_dim(1), plot_dim(2), jj),


if Sim.numeric
h = plot(Onres.(field).x,Onres.(field).Zss(:,1),'.','color','blue','DisplayName','numerical solution');
hold on

if Sim.numeric

plot(Onres.(field).x,Onres.(field).Zss(:,2),'.','color','cyan','DisplayName','numerical solution');
end

Expand All @@ -353,7 +357,7 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE)
end

xlabel(PLabel.(field));
ylabel('Z ( \Delta\omega_B )');
ylabel(['Z ( \Delta\omega=' sprintf('%.2f )',Sim.offset) ' ppm']);
set(gca,'XScale',param_xscale{jj});
set(gca,'YScale','lin');
str = sprintf('Z (%s)', PLabel.(field));
Expand Down Expand Up @@ -382,7 +386,7 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE)
end

xlabel(PLabel.(field));
ylabel('MTR_{asym} ( \Delta\omega_B )');
ylabel(['MTR_{asym} ( \Delta\omega=' sprintf('%.2f )',Sim.offset) ' ppm']);
set(gca,'XScale',param_xscale{jj});
set(gca,'YScale','lin');
str = sprintf('MTR_{asym} (%s)', PLabel.(field));
Expand Down
5 changes: 4 additions & 1 deletion general functions/init_Sim.m
View file
Expand Up @@ -20,6 +20,8 @@

% simulation initializations
Sim.Rex_sol = 'Hyper'; % cases: 'Hyper', 'Lorentz' , 'minilorentz'
Sim.MT_sol_type = 'Rex_MT'; % Rex_MT solution type - cases: 'Rex_MT'
Sim.MT_lineshape= 'Lorentzian';
Sim.spoilf = 0; % spoilingfactor (0 = full spoiling, 1= no spoiling)
%XXXXXXXXXXXXXXX DONT CHANGE ! XXXXXXXXXXXXXXXX
Sim.B1cwpe_quad = -1; % this is the B1 power equivalent flag: -1 means that P.B1 is the average B1 over a single pulse
Expand All @@ -33,4 +35,5 @@
Sim.linestomeasure=1;
Sim.flipangle= 90;
Sim.readout='FID';
Sim.TTM_rep = 0;
Sim.TTM_rep = 0;
Sim.DC=0.5;

0 comments on commit 6bf803d

Please sign in to comment.