function [p, q, Matrix_ls, Matrix_ms] = stochseq_DEBIPM_Manta(E_Ybad, E_Ygood); % [p, q, Matrix_ls, Matrix_ms] = stochseq_DEBIPM_Manta(0.5, 0.9); % check value of Ystdev steps = 6; % p: chance of moving from good to bad environment % q; % chance of moving from bad to good environment p=linspace(0,1,steps); q = linspace(0,1,steps); transient=200; % initial transient of length tlimit=1000; % Number of time steps trun=transient+tlimit; %time to run is desired length of sample path plus transients MatrixSize = 200; % number of columns and rows % DEB parameters Lmin = 130; % mm used to define limits of matrix approximation Lb = 380; % disc width at birth in cm Lm = 550; % maximum disc width in mm Lp = 7.29; % disc width at puberty in mm E_Ystdev=0.1; % sigma(Y) kappa = 0.80; % fraction energy allocation to respiration (as opposed to reproduction) mu = 0.05; % mortality rate (per year) rb = 0.18; % von Bertalanffy growth rate (per year) Rm = 1; % maximum reproduction of offspring produced per year at L_m for k = 1:length(p), disp(['lambda' num2str(k/length(p))]) % rows/ yaxis, p for l = 1:length(q) % columns / xaxis, q % habitat transition matrix H; 1=good; 2=bad; p379 Caswell 2001 H = [1-p(k) q(l); p(k) 1-q(l)]; % make Markov chain of environments Env=habitat_matrix(H,trun); %calculate vector of unperturbed stochastic growh rate [growth, loglambda_s, MeanSize] = stochseq(Env,E_Ygood,E_Ybad, MatrixSize,Lmin,Lp,rb,Rm,Lm,E_Ystdev,kappa,Lb,mu,transient,tlimit); Matrix_ls(k,l) = loglambda_s; Matrix_ms(k,l) = mean(MeanSize); end end % imagesc(q,p,Matrix_ls); colorbar; set(gca,'ydir','normal'); % q: x-axis; p: y-axis figure contour(q,p,Matrix_ls,'k','Showtext','on'); set(gca,'ydir','normal'); axis square; hold on; contour(q,p,Matrix_ls,[0,0],'r','LineWidth',2); xlabel('q: probability of moving to good environment') ylabel('p: probability of moving to bad environment') title('log lambda_s') figure contour(q,p,Matrix_ms,'k','Showtext','on'); set(gca,'ydir','normal'); axis square; xlabel('q: probability of moving to good environment') ylabel('p: probability of moving to bad environment') title('Mean body size') %%%%%%%%%%%%%%%%%%%% END OF FILE %%%%%%%%%%%%%%%%%%%%%%%%%%% function [growth, loglam_S, MeanSize]= stochseq(Env,E_Ygood, E_Ybad, MatrixSize,Lmin,Lp,rb, Rm,Lm,E_Ystdev,kappa,Lb,mu,transient,tlimit) % specify intial population structure vector and reproductive value vector vec1=15/MatrixSize*ones(1,MatrixSize); % 15 inds spread over all size bins vec1=vec1'; popvecAbs = vec1; trun=transient+tlimit; %time to run is desired length of sample path plus transients %set up a new cumulative number for all trun for i = 1:trun i2=Env(i); if i2==1, E_Y=E_Ygood; else E_Y=E_Ybad; end Linf = E_Y*Lm; % ultimate length [S, R, G, D, meshpts] = BigMatrixNotShrink(MatrixSize,Lmin,Lm,Lp,Rm,Lm,E_Y,E_Ystdev,rb,kappa,Lb,mu); mat1 = G*S + D*R; vec1 = mat1 * vec1; %numberE1 = sum(vec1.*(meshpts<1.1*Lmin)'); %numberJ1 = sum(vec1.*(meshpts>=1.1*Lmin & meshpts=Lp)'); growth1 = sum(vec1); vec1= vec1/growth1; % scale vec1 if( i > transient) %only store the ones after the transients i1 = i - transient; growth(i1) = growth1; [Mu, Var] = MeanVar(meshpts,vec1'); MeanSize(i1) = Mu; end end loglam_S=mean(log(growth)); % lam_S=exp(a); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%