From 6bf803d78d88d17327bc83d4155c121423a40be6 Mon Sep 17 00:00:00 2001 From: mzaiss Date: Tue, 10 Apr 2018 02:45:18 +0200 Subject: [PATCH] some minor changes fpor better functionality and simpler added BM simulation_paramspace --- doc/2_BMsimulation.m | 4 +- doc/3_BMsimulation_param_space.m | 166 +++++++++++++++++++++++++++++++ general functions/PLOT_SPACE.m | 24 +++-- general functions/init_Sim.m | 5 +- 4 files changed, 187 insertions(+), 12 deletions(-) create mode 100644 doc/3_BMsimulation_param_space.m diff --git a/doc/2_BMsimulation.m b/doc/2_BMsimulation.m index 72ee7ee..c6baf37 100644 --- a/doc/2_BMsimulation.m +++ b/doc/2_BMsimulation.m @@ -6,6 +6,7 @@ %% + Sim.n_cest_pool=1; Sim.dwA=0; Sim.R1A=1; @@ -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; diff --git a/doc/3_BMsimulation_param_space.m b/doc/3_BMsimulation_param_space.m new file mode 100644 index 0000000..efecc10 --- /dev/null +++ b/doc/3_BMsimulation_param_space.m @@ -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); + diff --git a/general functions/PLOT_SPACE.m b/general functions/PLOT_SPACE.m index e917c66..ac38798 100644 --- a/general functions/PLOT_SPACE.m +++ b/general functions/PLOT_SPACE.m @@ -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; @@ -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); @@ -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); @@ -195,6 +195,7 @@ function PLOT_SPACE(Sim,Space,NUMERIC_SPACE,ANALYTIC_SPACE) + % % % % % % % % % % % % AREX % % % % % % % % % % % % @@ -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); @@ -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 @@ -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 @@ -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)); @@ -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)); diff --git a/general functions/init_Sim.m b/general functions/init_Sim.m index dddc716..b67ad67 100644 --- a/general functions/init_Sim.m +++ b/general functions/init_Sim.m @@ -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 @@ -33,4 +35,5 @@ Sim.linestomeasure=1; Sim.flipangle= 90; Sim.readout='FID'; - Sim.TTM_rep = 0; \ No newline at end of file + Sim.TTM_rep = 0; + Sim.DC=0.5; \ No newline at end of file