%% Stochastic Simulations G=0.1:0.1:2.5; D=0.1:0.1:2.5; L=3; S=1; for m=1:length(G) for n=1:length(D) tic try parfor i=1:200 [N1,N2]=Stochastic_FKPP_withDeath(G(m),D(n),L,S,200,100); NN1D{m,n,i}=N1; NN2D{m,n,i}=N2; end catch end toc end end for m=1:length(G) for n=1:length(D) tic try parfor i=1:200 [N1,N2]=Stochastic_FKPP(G(m),D(n),L,S,200,100); NN1{m,n,i}=N1; NN2{m,n,i}=N2; end catch end toc end end %% Analysis: computation of fixation probability, based on total biomass trends. for m=1:length(G) for n=1:length(D) fix_prob=[]; parfor i=1:200 try if length(NN1{m,n,i})<100 fix_prob(i)=0; else p1=polyfit(1:100,NN1{m,n,i}(end-99:end),1); p2=polyfit(1:100,NN2{m,n,i}(end-99:end),1); fix_prob(i)=1-(p1(1)>p2(1)); end catch fix_prob(i)=NaN; end end score(m,n)=nanmean(fix_prob); fix_probD=[]; parfor i=1:200 try if length(NN1D{m,n,i})<100 fix_probD(i)=0; else p1=polyfit(1:100,NN1D{m,n,i}(end-99:end),1); p2=polyfit(1:100,NN2D{m,n,i}(end-99:end),1); fix_probD(i)=1-(p1(1)>p2(1)); end catch fix_probD(i)=NaN; end end scoreD(m,n)=nanmean(fix_probD); end end %% Plotting hold off subplot(1,2,1) S=scoreD; [DD,GG]=meshgrid(D,G); scatter(DD(:),GG(:),100,S(:),'filled') colormap jet axis xy axis square caxis([0 1]) xlabel('Diffusion') ylabel('Growth') colorbar axis([0 2.35 0 2.35]) set(gca,'xtick',[0 1 2 3],'ytick',[0 1 2 3],'fontsize',20,'box','on') title('With Death') subplot(1,2,2) S=score; [DD,GG]=meshgrid(D,G); scatter(DD(:),GG(:),100,S(:),'filled') colormap jet axis xy axis square caxis([0 1]) xlabel('Diffusion') ylabel('Growth') colorbar axis([0 2.35 0 2.35]) set(gca,'xtick',[0 1 2 3],'ytick',[0 1 2 3],'fontsize',20,'box','on') title('Without Death') function [N1,N2,Edge_ratio]=Stochastic_FKPP(G,D,LL,S,Ltotal,K) %% parameters L=Ltotal; % size of space dx=1; X=0:dx:L; % space g1 = 0.1; % growth rate of species 1 m1 = 0.3; % diffusion coefficient of species 1 lambda = sqrt(m1/g1); % lengthscale n=length(X); %% initialization c1=zeros(1,n); c1(X0 & X>1 & XK*0.05,1,'first'); t=t+1; end %% implantation c2=zeros(1,n); c2(edge+round(LL*lambda))=S; g2 = g1*G; m2 = m1*D; N1=0; % total biomass for species 1 N2=0; % total biomass for species 2 Edge_ratio=0; t=1; edge=0; while edge0 % Growth and Death growth1=poissrnd(max(0,g1*c1.*(K-c1-c2)/K)); growth2=poissrnd(max(0,g2*c2.*(K-c1-c2)/K)); c1 = c1 + growth1; c2 = c2 + growth2; % migration c1 NZ=find(c1>0 & X>1 & X0 & X>1 & XK/20,1,'first'); Ratio=c2(edge)/(c1(edge)+c2(edge)); % just in case c1(isnan(c1))=0; c2(isnan(c2))=0; N1(t)=sum(c1); % biomass species 1 N2(t)=sum(c2); % biomass species 2 Edge_ratio(t)=Ratio; t=t+1; end end function [N1,N2,Edge_ratio]=Stochastic_FKPP_withDeath(G,D,LL,S,Ltotal,K) %% parameters L=Ltotal; % size of space dx=1; X=0:dx:L; % space g1 = 0.1; % growth rate of species 1 m1 = 0.3; % diffusion coefficient of species 1 lambda = sqrt(m1/g1); % lengthscale n=length(X); %% initialization c1=zeros(1,n); c1(X0 & X>1 & XK*0.05,1,'first'); t=t+1; end %% implantation c2=zeros(1,n); c2(edge+round(LL*lambda))=S; g2 = g1*G; m2 = m1*D; N1=0; % total biomass for species 1 N2=0; % total biomass for species 2 Edge_ratio=0; t=1; edge=0; while edge0 % Growth growth1 = poissrnd(max(0,g1*c1)); growth2 = poissrnd(max(0,g2*c2)); c1 = c1 + growth1; c2 = c2 + growth2; % Death death1 = binornd(c1, max(0,g1.*(c1+c2)/K)); death2 = binornd(c2, max(0,g2.*(c1+c2)/K)); c1 = c1 - death1; c2 = c2 - death2; % migration c1 NZ=find(c1>0 & X>1 & X0 & X>1 & XK/20,1,'first'); Ratio=c2(edge)/(c1(edge)+c2(edge)); % just in case c1(isnan(c1))=0; c2(isnan(c2))=0; N1(t)=sum(c1); % biomass species 1 N2(t)=sum(c2); % biomass species 2 Edge_ratio(t)=Ratio; t=t+1; end end