Contents
% implement W-F process analogous to Y. Friedmann % simulate finite niche model (W-F process) clc; tic; % TIC, pair 1 %%%%%%%% get filename %%%%%%%%%%%%%%%%% p = mfilename('fullpath'); [pathstr, name, ext] = fileparts(p); rpath = [pathstr '/Results/']; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % run description time_end = 1e5; % number of simulated generations R = 1e6; % carrying capacity of ALL niches r = 10*1e-5; % per cell rate of gene transfer m_set = 0.02; % per cell rate of migration s_v = [0.025 0.030 0.034 0.038 0.04 0.045 0.05 0.1]; % vector of selection pressures % n = 1; % number of ecological niches in patch z = 1; % affinity to own niche nol = 1 - z; % 1 - the niche overlap n_firstC = 1000; % initial number of resistant cells Dm = 1; % diversity density of the incoming migrants % early initialise BIG matrixes and vectors R_groups = R; % use less space TOG = zeros(R_groups,2 + n); % BIG matrix lo_Mpops = false(R_groups,1); % BIG logical vector lo_Cpops = false(R_groups,1); % BIG logical vector lo_existing_cells = false(R_groups,1); % BIG logical vector F1 = zeros(R_groups,1); % BIG vector F2 = zeros(R_groups,1); % BIG vector % even create enough memory for linear index vectors li_Mpops = zeros(R_groups,1); % BIG vector (changing size) li_Cpops = zeros(R_groups,1); % BIG vector (changing size) li_selected_M_pops = zeros(R_groups,1); li_free_positions = zeros(R_groups,1); % BIG vector (changing size) first_free_positions = zeros(R_groups,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Situation % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% num_of_sim_repeats = 50; % number of simulation repeats for rep = 1 : num_of_sim_repeats % through repeats of the simulation for m_c = 1 : 2 % through the 2 cases of migration switch m_c case 1 m = 0; case 2 m = m_set; end for s_c = 1 : numel(s_v)
s = s_v(s_c);
rng('shuffle') % initiate random number generator
% initiate TOG (Table of Genotypes)
% use non-sparse matrix
new_bg = 2;% start at genotype 2 because 1 is the resistant one
next_niche = 1; % first niche to start with
count2Dminv = 0;
% generate possible niche associations
obtain_recource_mat = ones(n,n)*((nol)/(n-1));
obtain_recource_mat(eye(size(obtain_recource_mat))~=0) = (1-nol);
if n == 1 % if only a single niche
obtain_recource_mat = 1;
end
% seed TOG
% 5/s initially resistant guys, according to Kaplan1989
for through_TOG = 1 : (R - n_firstC) % through TOG
TOG(through_TOG,:) = [1, 0, obtain_recource_mat(next_niche,:)];
count2Dminv = count2Dminv + 1;
if count2Dminv == 1/Dm
new_bg = new_bg + 1;
count2Dminv = 0;
end
next_niche = next_niche + 1;
if next_niche > n
next_niche = 1;
end
end % through TOG
% finally add resistant guy
pos_first_Cpop = (R - n_firstC + 1);
niche_first_Cpop = 1;
TOG(pos_first_Cpop,:) = [0, n_firstC, obtain_recource_mat(niche_first_Cpop,:)];
% useful things
n_of_saves = 2; % how often do you want to save within simulated time
scope = time_end/n_of_saves; % scope to save TOS
counter_plus = scope*1;
% reset stuff
n_HGT = 0;
t = 0; counter = 0; counterp = 0;
sum_niche_com = zeros(1,n);
% saving stuff
TOG_save = zeros(R_groups,2+n); % BIG 2D matrix to save at end
t_vec = zeros(1,n_of_saves);
% get M (non-adapted) populations and C (adapted) populations,
% and numbers
% (important spead up: create logicals first (into pre-existing vector), then create linear
% indeces since this is faster when eccessing matix entries)
lo_Mpops = (TOG(:,1) ~= 0); % BIG vect (not changing sizes)
li_Mpops = find(lo_Mpops); % BIG vect (changing size)
lo_Cpops = (TOG(:,2) ~= 0); % BIG vect (not changing sizes)
li_Cpops = find(lo_Cpops); % BIG vect (changing size)
M = sum(TOG(li_Mpops,1));
C = sum(TOG(li_Cpops,2));
li_selected_M_pops = [];
first_free_positions = [];
mig_n = R*m*3;
%new_guy_M = [ones(mig_n,1),zeros(mig_n,1),ones(mig_n,1)]; %
%only for n = 1
n_mig = 0;
% MAIN LOOP
loopStart = tic; % TIC, pair loopStart (time prediction)
keepLooping = true;
while keepLooping % THROUGH TIME
% what you want to save
if counter == 0 || (counter >= counter_old + counter_plus)
counterp = counterp + 1;
t_vec(counterp) = t; % save the time
TOG_save(:,:,counterp) = TOG; % save TOG
counter_old = counter;
end
% update time and counter
t = t + 1; % move through steps of generations
counter = counter + 1; % count generations
% get summed (total) competetive weights in each of the n niches
% for each niche sum over existing pops:(popssize*popsfitness*popsniche_association)
sum_niche_com = sum(bsxfun(@times,TOG([li_Mpops; first_free_positions],1),TOG([li_Mpops; first_free_positions],3:end)),1) +...
sum(bsxfun(@times,(TOG([li_Cpops; li_selected_M_pops],2)*(1+s)),TOG([li_Cpops; li_selected_M_pops],3:end)),1);
% how much resources are obtained (per cell in every
% existing genotype) = growth rate (per cell in every existing genotype)
Rovern = R/n;
F1([li_Mpops; first_free_positions]) = Rovern*sum(bsxfun(@ldivide,sum_niche_com,TOG([li_Mpops; first_free_positions],3:end)),2);
F2([li_Cpops; li_selected_M_pops]) = Rovern*sum(bsxfun(@ldivide,sum_niche_com,TOG([li_Cpops; li_selected_M_pops],3:end)*(1+s)),2);
% get mean new numbers+poiss in one step
% account for loss through m!!
TOG([li_Mpops; first_free_positions],1) = ...
poissrnd(F1([li_Mpops; first_free_positions]).*(1 - m).*TOG([li_Mpops; first_free_positions],1)); % only use poiss where neccess
TOG([li_Cpops; li_selected_M_pops;pos_first_Cpop],2) = ...
[poissrnd(F2([li_Cpops; li_selected_M_pops]).*(1 - m).*TOG([li_Cpops; li_selected_M_pops],2));...
round(F2(pos_first_Cpop).*(1-m).*TOG(pos_first_Cpop,2))]; % only use poiss where neccess
% is it unlikely for many genotypes, but it might N > R by
% chance
% get existing M and C populations
lo_Mpops = (TOG(:,1) ~= 0); % BIG vect (not changing sizes)
li_Mpops = find(lo_Mpops); % BIG vect (changing size)
lo_Cpops = (TOG(:,2) ~= 0); % BIG vect (not changing sizes)
li_Cpops = find(lo_Cpops); % BIG vect (changing size)
M = sum(TOG(li_Mpops,1)); C = sum(TOG(li_Cpops,2)); % = sum(TOG(li_Cpops,2)); % get cell numbers M and C
% several HGT event per generation expected (without new migrators)
% HGT happens the same now as in exact simulation, weighted by cell groups in M,
% where genotypes with more cells have higher change to HG
num_o_HGT = poissrnd(M*C/(M+C)*r); % how many HGT events?
if num_o_HGT > 0
num_Mpops = numel(li_Mpops);
if num_o_HGT > num_Mpops
num_o_HGT = num_Mpops;
disp('Corrected HGT event number'); % only happens for small R < 1e4
end
if numel(li_Mpops) == 1 % if only single M population
li_selected_M_pops = ones(num_o_HGT,1)*li_Mpops;
else li_selected_M_pops = randsample(li_Mpops,num_o_HGT,true,TOG(li_Mpops,1));
end
TOG(li_selected_M_pops,1) = TOG(li_selected_M_pops,1) - 1; % remove a cell from r_M's populations
M = M - num_o_HGT;
TOG(li_selected_M_pops,2) = TOG(li_selected_M_pops,2) + 1; % add a cell to r_M's C populations
C = C + num_o_HGT;
n_HGT = n_HGT + num_o_HGT;
else li_selected_M_pops = [];
end
% migration (I know how many cells are missing due to migration)
% find free spaces (before HGT, are also free after HGT!)
lo_existing_cells = (lo_Mpops|lo_Cpops); % BIG vector (not changing size)
li_free_positions = find(lo_existing_cells == 0); % BIG vector (lin indeces, changing size)
% % LOOK into: number of used space
% used_spaces = sum(lo_existing_cells);
% disp([num2str(disp_num_of_cells) ': ' num2str(used_spaces) ' spaces used'])
% % LOOK
n_mig = R - M - C;
if n_mig > 0
niche_vec = [next_niche:(n_mig + next_niche -1)] - floor(([next_niche:(n_mig+next_niche-1)]-1)/n)*n;
next_niche = niche_vec(end) + 1;
new_guys = [ones(n_mig,1),zeros(n_mig,1),obtain_recource_mat(niche_vec',:)];
first_free_positions = li_free_positions(1:n_mig);
TOG(first_free_positions,:) = new_guys;
% update numbers
M = M + n_mig;
else first_free_positions = [];
end
% 6 Decide if time is over or go back and sample new time-step
if t >= time_end
keepLooping = false;
end
% Time prediction unit
if counter == 500
looptimeXtimes = toc(loopStart);
time_in_sec = looptimeXtimes/500 * (time_end - 500); % expected time in seconds
time_in_hms = secs2hms(time_in_sec); % expected time in hours,minutes,secs
display1 = [datestr(now) ' (' num2str(rep) ') ' ' S = ' num2str(s) '; Time left: ' time_in_hms]; % print expected time
disp(display1);
end
end % THROUGH TIME
Saving
savename = [rpath name '_rep[' num2str(rep) ']_S[' num2str(s)... ']_m[' num2str(m) ']_r[' num2str(r) ']_t[' num2str(time_end) ']_result.mat']; save(savename,'TOG_save','-v7.3',... 't_vec','counterp','n_HGT','R','R_groups');
end % through selection end % through migration end % through reps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
From here can be run without simulation if previously run %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% Evaluate Results %%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mys_set = s_v; rep_set = 1 : num_of_sim_repeats; % reset all outputs % you need the same amount of reps for all; for to_m = 2 % through m if to_m == 1 m = 0; elseif to_m == 2 m = 0.02; end for to_c = 1 : numel(mys_set) % through selection s = mys_set(to_c); for to_repeat = 1 : numel(rep_set) % through repeats
rep = rep_set(to_repeat);
% load
load_name = [rpath name '_rep[' num2str(rep) ']_S[' num2str(s)...
']_m[' num2str(m) ']_r[' num2str(r) ']_t[' num2str(time_end) ']_result.mat'];
load(load_name);
% reset
zxc = zeros(1,counterp); C_save = zxc; M_save = zxc;
HGP_diversity2 = zxc;
for i = 1 : counterp % go to end of TOS
% M and C fractions
num_M = sum(TOG_save(:,1,i));
num_C = sum(TOG_save(:,2,i));
M_save(i) = num_M/(num_M + num_C);
C_save(i) = num_C/(num_M + num_C);
% True Diversity q = 2
% get M (non-adapted) populations and C (adapted) populations
lo_Mpops = (TOG_save(:,1,i) ~= 0); % BIG vect (not changing sizes)
li_Mpops = find(lo_Mpops); % BIG vect (changing size)
lo_Cpops = (TOG_save(:,2,i) ~= 0); % BIG vect (not changing sizes)
li_Cpops = find(lo_Cpops); % BIG vect (changing size)
all_BGs = TOG_save(:,1,i) + TOG_save(:,2,i);
BG_simps = sum((all_BGs(lo_Mpops|lo_Cpops)/(num_M+num_C)).^2);
FL_simps = sum(([TOG_save(li_Mpops,1,i);num_C]/(num_M+num_C)).^2);
HGP_diversity2(i) = FL_simps/BG_simps;
% Check how we are doing with space in TOG
num_used_slots = sum(all_BGs == 1);
my_disptext = [num2str(num_used_slots) ' out of ' num2str(R_groups) ' used'];
%disp(my_disptext);
end % from TOS-save
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% OUTPUT %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Real HGT flow
realHGTflow(to_c,to_repeat,to_m) = n_HGT;
Final DR
finalDR(to_c,to_repeat,to_m) = HGP_diversity2(counterp);
disp(['For R = ' num2str(R) ', s = ' num2str(s) ' u get DR = ' num2str(finalDR(to_c,to_repeat,1)) '!']);
Plot one simulation
if 1 == 0 magicfigure_single_line(HGP_diversity2,'DR over time'); %savefig([name '02'],gcf,'EPS'); end
end % through repeats end % through selection end % through m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Plotting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% realHGTflow_mean = mean(realHGTflow,2); realHGTflow_sd = std(realHGTflow,0,2); % finalDR_mean = mean(finalDR,2); finalDR_sd = std(finalDR,0,2);
Create figure
close all XMatrix1 = mys_set; YMatrix1 = finalDR_mean(:,:,1)'; YMatrix2 = finalDR_mean(:,:,2)'; EMatrix1 = finalDR_sd(:,:,1)'; EMatrix2 = finalDR_sd(:,:,2)'; myylabel = 'final DR'; axisfixed = 0; % Defaults values width = 3; % Width in inches height = 3; % Height in inches alw = 2; % AxesLineWidth fsz = 20; % Fontsize lw = 2.5; % LineWidth msz = 8; % MarkerSize lbs = 28; % Label Fontsize % Create figure figure; pos = get(gcf, 'Position'); set(gcf, 'Position', [pos(1) pos(2) width*115, height*115]); %<- Set size set(gca, 'FontSize', fsz, 'LineWidth', alw, 'FontName','Calibri'); %<- Set properties % Create plot % create pretty filled bit hold on fill( [XMatrix1 fliplr(XMatrix1)], [(YMatrix1 - EMatrix1) fliplr(YMatrix1 + EMatrix1)],... [0 0.6392 0.2039],'EdgeColor','w','FaceAlpha', 0.3); % This is the seethroughness of the pretty filled bit (1 = not see through) plot1 = plot(XMatrix1, YMatrix1, 'LineWidth', lw); % Plots main line set(plot1(1),'DisplayName','Simulation without migration',... 'Color',[0 124 175]/255); % fill( [XMatrix1 fliplr(XMatrix1)], [(YMatrix2 - EMatrix2) fliplr(YMatrix2 + EMatrix2)],... [1.0000 0.1608 0.1373],'EdgeColor','w','FaceAlpha', 0.3); % This is the seethroughness of the pretty filled bit (1 = not see through) plot2 = plot(XMatrix1, YMatrix2, 'LineWidth', lw); % Plots main line set(plot2(1),'DisplayName','Simulation with migration',... 'Color',[1.0000 0.1608 0.1373]); set(gca,'Layer','top') % put axis in front again % Create xlabel xlabel('S','FontSize',lbs); % Create ylabel ylabel(myylabel,'FontSize',lbs); xlim([0.023 0.102]); if axisfixed == 1 % x and ylim xlim([0 0.10]); % Set Tick Marks %set(gca,'XTick',0:0.05:0.3); %set(gca,'YTick',[0 1 2 3 4]*1e-3); firstYlim = ylim; ylim([-40 1400]) end if axisfixed == 2 % x and ylim xlim([0 0.10]); % Set Tick Marks %set(gca,'XTick',0:0.05:0.3); %set(gca,'YTick',[0 1 2 3 4]*1e-3); firstYlim = ylim; ylim([-1.4 firstYlim(2)]) end if 0 == 1 % Create legend legend1 = legend([plot1,plot2],'Simulation without migration','Simulation with migration'); set(legend1,'EdgeColor',[1 1 1],'Location','Best','box','off','YColor',[1 1 1],... 'XColor',[1 1 1],... 'Color',[1 1 1]); end % Here we preserve the size of the image when we save it. set(gcf,'InvertHardcopy','on'); set(gcf,'PaperUnits', 'inches'); papersize = get(gcf, 'PaperSize'); left = (papersize(1)- width)/2; bottom = (papersize(2)- height)/2; myfiguresize = [left, bottom, width, height]; set(gcf,'PaperPosition', myfiguresize); % make this so that saved graph looks the same as like on the screen set(gcf, 'PaperPositionMode', 'auto'); savename = [name '02.eps']; print('-depsc2',savename) %%%%%%%%%%%% show elapsed time %%%%%%%%%%%%%%%%%%%%%%%% time_whole = toc; % TOC, pair 1 time_whole_hms = secs2hms(time_whole); display2 = ['Elapsed time is ' time_whole_hms]; % print expected time disp(display2);