function [p, q, Matrix_loglambda_sens1, Matrix_loglambda_sens2, Matrix_loglambda_sens1_val, Matrix_loglambda_sens2_val] = stochseq_DEBIPM_sens_Orchestia(E_Ybad, E_Ygood); % pars perturbed are: Lb, Lp, Lm, Rm, rb, kappa, mu 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=500; % 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 = 3; % mm used to define limits of matrix approximation Lb = 3.79; % length at birth in mm Lm = 15.61; % maximum body length in mm Lp = 7.29; % length at puberty in mm E_Ystdev=0.5; % sigma(Y) kappa = 0.80; % fraction energy allocation to respiration (as opposed to reproduction) mu = 0.27; % mortality rate (per month) rb = 0.13; % von Bertalanffy growth rate (per month) Rm = 32; % maximum reproduction of offspring produced per month at L_m (per month) fracperturb = 0.01; % perturbation of DEB pars count=ones(1,7); % Lb, Lp, Lm, kappa, mu, rb, Rm for k = 1:length(p), disp(['count' 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 loglambda_s_orig = stochseq(Env,E_Ygood,E_Ybad, MatrixSize,Lmin,Rm,Lm,Lp,E_Ystdev,kappa,Lb,mu,rb,transient,tlimit); for i=1:length(count); if i==6, LogLambda_s_perturb(i,1)=loglambda_s_orig; % skip kappa else perturbation = ones(1,7); perturbation(i) = perturbation(i) + fracperturb; loglambda_s_pert = stochseqperturb(Env,E_Ygood,E_Ybad, MatrixSize,Lmin,Rm,Lm,Lp,E_Ystdev,kappa,Lb,mu,rb,transient,tlimit, perturbation); LogLambda_s_perturb(i,1) = loglambda_s_pert; end end %calculate highest elasticityvalue LogLambda_s_Sens = (LogLambda_s_perturb - loglambda_s_orig)./fracperturb; LogLambda_s_Sens = abs(LogLambda_s_Sens); % use absolute elasticity values LogLambda_s_Sens2 = LogLambda_s_Sens; % also explore second highest elasticity value [a b] = max(LogLambda_s_Sens); if b==1, Matrix_loglambda_sens1(k,l) = 1; Matrix_loglambda_sens1_val(k,l) = LogLambda_s_Sens(b); LogLambda_s_Sens2(b)=-100; elseif b==2; Matrix_loglambda_sens1(k,l)=2; Matrix_loglambda_sens1_val(k,l) = LogLambda_s_Sens(b); LogLambda_s_Sens2(b)=-100; elseif b==3, Matrix_loglambda_sens1(k,l)=3; Matrix_loglambda_sens1_val(k,l) = LogLambda_s_Sens(b); LogLambda_s_Sens2(b)=-100; elseif b==4, Matrix_loglambda_sens1(k,l)=4; Matrix_loglambda_sens1_val(k,l) = LogLambda_s_Sens(b); LogLambda_s_Sens2(b)=-100; elseif b==5, Matrix_loglambda_sens1(k,l)=5; Matrix_loglambda_sens1_val(k,l) = LogLambda_s_Sens(b); LogLambda_s_Sens2(b)=-100; elseif b==6, Matrix_loglambda_sens1(k,l)=6; Matrix_loglambda_sens1_val(k,l) = LogLambda_s_Sens(b); LogLambda_s_Sens2(b)=-100; else Matrix_loglambda_sens1(k,l)=7; Matrix_loglambda_sens1_val(k,l) = LogLambda_s_Sens(b); LogLambda_s_Sens2(b)=-100; end; % rows are p (yaxis) and colums are q (xaxis) %calculate second highest elasticity value [a b] = max(LogLambda_s_Sens2); if b==1, Matrix_loglambda_sens2(k,l) = 1; Matrix_loglambda_sens2_val(k,l) = LogLambda_s_Sens2(b); elseif b==2; Matrix_loglambda_sens2(k,l)=2; Matrix_loglambda_sens2_val(k,l) = LogLambda_s_Sens2(b); elseif b==3, Matrix_loglambda_sens2(k,l)=3; Matrix_loglambda_sens2_val(k,l) = LogLambda_s_Sens2(b); elseif b==4, Matrix_loglambda_sens2(k,l)=4; Matrix_loglambda_sens2_val(k,l) = LogLambda_s_Sens2(b); elseif b==5, Matrix_loglambda_sens2(k,l)=5; Matrix_loglambda_sens2_val(k,l) = LogLambda_s_Sens2(b); elseif b==6, Matrix_loglambda_sens2(k,l)=6; Matrix_loglambda_sens2_val(k,l) = LogLambda_s_Sens2(b); else Matrix_loglambda_sens2(k,l)=7; Matrix_loglambda_sens2_val(k,l) = LogLambda_s_Sens2(b); end; % rows are p (yaxis) and colums are q (xaxis) end end imagesc(q,p,Matrix_loglambda_sens1); colorbar; set(gca,'ydir','normal'); % q: x-axis; p: y-axis xlabel('q: probability of moving to good environment') ylabel('p: probability of moving to bad environment') title('highest elasticity of loglambda_s') figure imagesc(q,p,Matrix_loglambda_sens2); colorbar; set(gca,'ydir','normal'); % q: x-axis; p: y-axis xlabel('q: probability of moving to good environment') ylabel('p: probability of moving to bad environment') title('second highest elasticity of loglambda_s') %%%%%%%%%%%%%%%%%%%% END OF FILE %%%%%%%%%%%%%%%%%%%%%%%%%%% function loglam_S= stochseq(Env,E_Ygood, E_Ybad, MatrixSize,Lmin,Rm,Lm,Lp,E_Ystdev,kappa,Lb,mu,rb,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 [S, R, G, D, meshpts] = BigMatrixShrink(MatrixSize,Lmin,Lm,Lp,Rm,Lm,E_Y,E_Ystdev,rb,kappa,Lb,mu); mat1 = G*S + D*R; vec1 = mat1 * vec1; growth1 = sum(vec1); vec1= vec1/growth1; if( i > transient) %only store the ones after the transients i1 = i - transient; growth(i1) = growth1; end end loglam_S=mean(log(growth)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function loglam_S= stochseqperturb(Env,E_Ygood, E_Ybad, MatrixSize,Lmin,Rm,Lm,Lp,E_Ystdev,kappa,Lb,mu,rb,transient,tlimit,perturbation) % 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 [S, R, G, D, meshpts] = BigMatrixShrinkSens(MatrixSize,Lmin,Lm,Lp,Rm,Lm,E_Y,E_Ystdev,rb,kappa,Lb,mu,perturbation); mat1 = G*S + D*R; vec1 = mat1 * vec1; growth1 = sum(vec1); vec1= vec1/growth1; if( i > transient) %only store the ones after the transients i1 = i - transient; growth(i1) = growth1; end end loglam_S=mean(log(growth)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%