#Script for Karatayev, Kraft, and Zipkin 2015 Ecosphere


rm(list = ls())


###TWO CORE FUNCTIONS

#Define four functions (see Appendix A) for how maturation rate might change with abundance
min=0.2; #set minimum maturation rate (will vary this later)

mDD=function(x,max,DDtype,alpha,sj,sa,mDDstage,recruitType){
if(mDDstage==2){
	#Calculate equilibrium abundance of adults (K) at zero harvest
	#See Equations 5 and 9 in text and Appendix B, respectively (dependeing on Ricker or Beverton-Holt recruitment), with beta=1
	K=c(abs(log( (min*sj*alpha) / ((1-sj*(1-min))*(1-sa)) )),   abs(((min*sj*alpha) / ((1-sj*(1-min))*(1-sa))) -1))[recruitType]
} else {
	#For juvenile-dependent maturation - picked a value of K  which is slightly above equilibrium juvenile abundance at zero harvest
	#This value is valid for parameters in Fig. 2d (high survival, alpha=10, min=0.2, max=1), and will need to be re-adjusted for other parameters or model forms
	K=2.84985 
}

	y=matrix(0,length(x),4)
	for(i in 1:length(x)) ifelse((x[i]<K), (y[i,1]=max-(((max-min)/K)*x[i])), (y[i,1]=min))  ##Type I, Equation 4 in text
	for(i in 1:length(x)) ifelse((x[i]<K), (y[i,2]=max*K/(K+((max/min)-1)*x[i])), (y[i,2]=min))  ##Type II, Equation 6 in Appendix A
	for(i in 1:length(x)) ifelse((x[i]<K), (y[i,3]=min+(max-min)*sqrt(1-(((1/K)*x[i])^(2*(1-min)))) ), (y[i,3]=min))  ##Type III, Equation 7 in Appendix A
	y[,4]=max   ##Type 0 (constant maturation), with max representing mu
return(as.vector(y[,DDtype])) };





#Main simulation function
run=function(params){
ntrials=params[1]; runlength=params[2]; sj=params[3]; sa=params[4]; alpha=params[5]; beta=params[6:7]; phase=params[14]; #Core parameters
recruitType=params[8]; harvMod=params[9:10]; mDDstage=params[11]; max=params[15]; DDtype=params[16]; #Parameters exploring model variations
harvest = seq(from=0, to=1, length.out=ntrials)
hold = matrix(0,ntrials,4); hold[,1] = harvest;
	for(x in 1:ntrials){
		harvest=c(hold[x,1],hold[x,1])*harvMod #Set harvest level for simulation; adjust for stage-specific harvest (if present). hj=harvest[1], ha=harvest[2]
		J=double(runlength); J[1]=2;
		A=double(runlength); A[1]=1;
		mrec=double(runlength);
			for(t in 2:runlength){
				m=as.numeric( c(mDD(J[t-1],max,DDtype,alpha,sj,sa,mDDstage,recruitType), mDD(A[t-1],max,DDtype,alpha,sj,sa,mDDstage,recruitType))[mDDstage] )
				gamma=c( alpha*exp(-beta[1]*A[t-1]*(1-harvest[2])), alpha/(1+beta[2]*A[t-1]*(1-harvest[2])) )[recruitType]  #Note: in all simulations, both beta values =1
				Juves = A[t-1]*(1-harvest[2])*gamma + J[t-1]*(1-harvest[1])*sj*(1-m)
				Adults = J[t-1]*(1-harvest[1])*sj*m + A[t-1]*(1-harvest[2])*sa
				J[t]=ifelse(Juves>0.001,Juves,0); A[t]=ifelse(Adults>0.001,Adults,0); mrec[t] = m;
			}
	hold[x,2]=mean(A[phase:runlength]); hold[x,3]=mean(J[phase:runlength]); 
	hold[x,4]=mean(mrec[phase:(runlength-1)]);   #mean maturation rate (for reference)
	}
return(hold);  }


#Basic parameters (edited as needed later):
survival1=c(0.2,0.2); survival2=c(0.2,0.8); survival3=c(0.8,0.2); survival4=c(0.8,0.8);
#ntrials, runlength, [sj,sa], alpha, betaRickerDD,  betaBevertonHoltDD, DD_Type(1=Ricker|2=BevertonHolt), HarvestStrategy(J|A), mDDstage(J|A), [ignored parameters], #TransientStepsToIgnore 
params=c(30, 150,   survival4, 10,        1,                   1,                      1,                         c(1,1),              2,            1,      0,             100) 



#####################Figure 2
#This function just calls the run function twice (once for constant and once for density dependent maturation)
#Then standardizes equilibrium adult and juvenile abundance for each harvest rate to the total population size at zero harvest
#That allows us to relate change in abundance at a certain harvest level to abundance observed before population control was started

mDD2.Sim=function(params){  ntrials=params[1];
Stuffy = matrix(0, ntrials, 13); harvest=seq(0,1,length=ntrials); Stuffy[,1]=harvest;
type1.df=run(c(params, max=1, DDtype=1)); 
Stuffy[,2]=type1.df[,2]/sum(type1.df[1,2:3]); 
Stuffy[,3]=type1.df[,3]/sum(type1.df[1,2:3]); 
Stuffy[,4]=type1.df[,4];

const.df=run(c(params, max=0.5, DDtype=4));
Stuffy[,11]=const.df[,2]/sum(const.df[1,2:3]); 
Stuffy[,12]=const.df[,3]/sum(const.df[1,2:3]); 
Stuffy[,13]=const.df[,4];
Out.df=data.frame(Stuffy); colnames(Out.df)=c("h","type1_A","type1_J","type1_A.m","na","na","na","na","na","na","const_A", "const_J", "const_A.m");
return(Out.df) }


#Calculate how equilibrium population size changes with harvest intensity for density-dependent and constant maturation
params[1:8]=c(50,150,survival1,10,1,1, 1);  params[9:11]=c(1,1, 2); s1=mDD2.Sim(params);
params[1:8]=c(50,150,survival2,10,1,1, 1);  params[9:11]=c(1,1, 2); s2=mDD2.Sim(params);
params[1:8]=c(50,150,survival3,10,1,1, 1);  params[9:11]=c(1,1, 2); s3=mDD2.Sim(params);
params[1:8]=c(50,150,survival4,10,1,1, 1);  params[9:11]=c(1,1, 2); s4=mDD2.Sim(params);
params[1:8]=c(50,150,survival4,20,1,1, 1);  params[9:11]=c(1,1, 2); s4_20=mDD2.Sim(params);  #Case with alpha=20, Fig 2b

#Adding maturation dependence on juveniles and stage-specific harvest
params[1:8]=c(50,150,survival4,10,1,1, 1);  params[9:11]=c(0,1, 2); s4ha=mDD2.Sim(params);
params[1:8]=c(50,150,survival4,10,1,1, 1);  params[9:11]=c(1,0, 2); s4hj=mDD2.Sim(params);
params[1:8]=c(50,150,survival4,10,1,1, 1);  params[9:11]=c(1,1, 1); s4dj=mDD2.Sim(params);
params[1:8]=c(50,150,survival4,10,1,1, 1);  params[9:11]=c(0,1, 1); s4djha=mDD2.Sim(params);
params[1:8]=c(50,150,survival4,10,1,1, 1);  params[9:11]=c(1,0, 1); s4djhj=mDD2.Sim(params);

#Find the biggest population size in each simulation for density-dependent and constant maturation, and assemble into a matrix
DDm.A=round(matrix(c(   max(s4[,11]+s4[,12]),max(s4[,2]+s4[,3]), max(s4ha[,11]+s4ha[,12]),max(s4ha[,2]+s4ha[,3]), max(s4hj[,11]+s4hj[,12]),max(s4hj[,2]+s4hj[,3])), nr=2),2)
DDm.J=round(matrix(c(   max(s4dj[,11]+s4dj[,12]),max(s4dj[,2]+s4dj[,3]), max(s4djha[,11]+s4djha[,12]),max(s4djha[,2]+s4djha[,3]), max(s4djhj[,11]+s4djhj[,12]),max(s4djhj[,2]+s4djhj[,3])), nr=2),2)



#Plotting:
dev.new(width=10,height=5.5); par(mfrow=c(1,2)); ymax=2; xmax=0.8;
plot((s1[,11]+s1[,12])~s1[,1],    type="b",col=colors()[228],lwd=3,lty=1,pch=1, cex=1.25, xlim=c(0,xmax), ylim=c(0,2),ylab="Equilibrium population size",xlab="Proportion of total population harvested",cex.lab=1.35); #text(0.05,1.9, "A:",cex=2);
points((s2[,11]+s2[,12])~s2[,1],type="b",col=colors()[203],lwd=3,lty=1,pch=1, cex=1.25); points((s3[,11]+s3[,12])~s3[,1],type="b",col=colors()[178],lwd=3,lty=1,pch=1, cex=1.25); points((s4[,11]+s4[,12])~s4[,1],type="b",col=colors()[153],lwd=3,lty=1,pch=1, cex=1.25);
abline(1,0,col=1,lty=2,lwd=3)
legend("topright", seg.len=2.5,lty=c(1),pch=c(1),lwd=3,pt.cex=1.5,cex=1.25,box.col=0,rev(c(expression("s"[J]~"="~"s"[A]~"="~"0.2"), expression("s"[J]~"="~"0.2,  s"[A]~"="~"0.8"), expression("s"[J]~"="~"0.8,  s"[A]~"="~"0.2"), expression("s"[J]~"="~"s"[A]~"="~"0.8"))), col=rev(colors()[c(228,203,178,153)]));

plot((s4[,11]+s4[,12])~s4[,1],    type="b",col=8,lwd=3,lty=1,pch=1, cex=1.25, xlim=c(0,xmax), ylim=c(0,ymax),ylab=" ",xlab="Proportion of total population harvested",cex.lab=1.35); #text(0.05,1.9, "B:",cex=2);
points((s4[,2]+s4[,3])~s4[,1],    type="b",col=8,lwd=2,lty=3,pch=17,cex=1.25); points((s4_20[,11]+s4_20[,12])~s4_20[,1],type="b",col=1,lwd=3,lty=1,pch=1, cex=1.25); points((s4_20[,2]+s4_20[,3])~s4_20[,1],    type="b",col=1,lwd=2,lty=3,pch=17,cex=1.25);
abline(1,0,col=1,lty=2,lwd=3)
legend(x=.1,y=2.08, c("Density-dependent","Constant"), lty=c(3,1), pch=c(17,1), col=c(1,1), lwd=2.5, cex=1.2, pt.cex=1.35, box.col=0); 
text(0.45, 0.65, expression(alpha~"=10"), cex=1.5, col=colors()[203]); text(0.75, 1.13, expression(alpha~"=20"), cex=1.5);


dev.new(width=12,height=4.8);  par(mfrow=c(1,2));
barplot(DDm.A, beside=T, col=c(1,8), ylim=c(1,2),xpd=F,ylab="Strength of overcompensation", xlab="Removal strategy", names.arg=c("Both", "Adults only", "Juveniles only"), cex.names=1.1 ,cex.lab=1.35,main="Adult-limited maturation")
legend(c(1,2.05),c("Constant","Density-dependent"), pch=15, col=c(1,8), box.col=0, pt.cex=1.8,cex=1.25); #text(1.3,2.3, "C:",cex=2);
barplot(DDm.J, beside=T, col=c(1,8), ylim=c(1,2),xpd=F,ylab=" ", xlab="Removal strategy", names.arg=c("Both", "Adults only", "Juveniles only"),cex.names=1.1,cex.lab=1.35,main="Juvenile-limited maturation")







#####################Figure 3
#This function calls the run function steps=16 times for density dependent maturation, each time with an increasing level of maximum maturation rate
#Then for each trial we look for the harvest levels at which each kind of population (a) collapses, and (b) adult abundance declines by 50%
MvsExtinct.Sim=function(params){     ntrials=params[1];
harvest = seq(from=0, to=1, length.out=ntrials);
steps=16  #number of scalings done
max = seq(from=min, to=1, length.out=steps);
hold = matrix(0, steps, 5); hold[,1] = max;

for(y in 1:steps){
	type1.df=run(c(params, max[y], DDtype=1))
	hold[y,2]=max(type1.df[,4]);
	hold[y,3]=harvest[which.min(subset(type1.df[,3], type1.df[,3]>0))]
	hold[y,4]=harvest[ which.min(abs((type1.df[,2]/type1.df[1,2])-0.5)) ]
	hold[y,5]=max(rowSums(type1.df[,2:3])/sum(type1.df[1,2:3]))
}
Out.df=data.frame(hold)
colnames(Out.df)=c("mmax","I_mmax","I_hext","I_Ahalf","I_Npeak")
return(Out.df) }


#New minimum maturation rate (at lower min values, low-survival populations can't persist)
min=0.4; max=1;
params[1:11]=c(40,150,survival1,10,1,1,2,  1,1, 2); s1=MvsExtinct.Sim(params);
params[1:11]=c(40,150,survival1,10,1,1,2,  0,1, 2); s1ha=MvsExtinct.Sim(params);

#New minimum maturation rate (for populations with higher survival rates)
min=0.1
params[1:8]=c(40,150,survival2,10,1,1,2); params[9:11]=c(1,1, 2); s2=MvsExtinct.Sim(params);
params[1:8]=c(40,150,survival3,10,1,1,2); params[9:11]=c(1,1, 2); s3=MvsExtinct.Sim(params);
params[1:8]=c(40,150,survival4,10,1,1,2); params[9:11]=c(1,1, 2); s4=MvsExtinct.Sim(params);
params[1:8]=c(40,150,survival2,10,1,1,2); params[9:11]=c(0,1, 2); s2ha=MvsExtinct.Sim(params);
params[1:8]=c(40,150,survival3,10,1,1,2); params[9:11]=c(0,1, 2); s3ha=MvsExtinct.Sim(params);
params[1:8]=c(40,150,survival4,10,1,1,2); params[9:11]=c(0,1, 2); s4ha=MvsExtinct.Sim(params);


#Plotting
dev.new(width=9.91,height=7.8); par(mfrow=c(2,2),mai=c(1,1,1,1),oma=c(0,0,2,0),mar=c(5,6,0.25,2),lheight=1); clr=c(1,8,1,8); lt=c(1,1,2,2); #clr=c(1,2,4,3);
xlab1=expression(Minimum~age~at~maturity~(~italic(m[max])~""^{-1})); ylab1=expression("Harvest level at population collapse (h"["    c"]~"),"); ylab2=expression("or 50% adult decline (h"["   half"]~")");
plotlab=c(expression("A:  s"[J]~"="~"s"[A]~"="~"0.2"), expression("B:  s"[J]~"="~"0.2,  s"[A]~"="~"0.8"), expression("C:  s"[J]~"="~"0.8,  s"[A]~"="~"0.2"), expression("D:  s"[J]~"="~"s"[A]~"="~"0.8"));

plot(supsmu((1/s1[,2]), s1[,3]), col=clr[1], lty=lt[1],ylim=c(0,0.65),cex.lab=1.3,cex.axis=1.2, ylab="", xlim=c(1,10), xlab="", lwd=4, type="l");   text(3.2,.62,plotlab[1],cex=1.6); mtext(2,text=ylab1,line=4); mtext(2,text=ylab2,line=3);
points(supsmu((1/s1ha[,2]), s1ha[,3]), lty=lt[2],col=clr[2],lwd=4,type="l"); points(supsmu(1/s1[,2], s1[,4]), lty=lt[3],col=clr[3],lwd=4,type="l"); points(supsmu(1/s1ha[,2], s1ha[,4]), lty=lt[4], col=clr[4],lwd=4,type="l"); 
legend(x=2.25,y=0.55, seg.len=2.5, cex=1.15,box.col=0,c("Population collapse","Population collapse, Adult harvest", "50% adult decline", "50% adult decline, Adult harvest"), col=clr, lwd=4, lty=lt);

plot(supsmu((1/s2[,2]), s2[,3]), col=clr[1], lty=lt[1],ylim=c(0,0.65),cex.lab=1.3,cex.axis=1.2, ylab="", xlim=c(1,10), xlab="", lwd=4, type="l");   text(3.9,.62,plotlab[2],cex=1.6); #mtext(2,text=ylab1,line=4); mtext(2,text=ylab2,line=3);
points(supsmu((1/s2ha[,2]), s2ha[,3]), lty=lt[2],col=clr[2],lwd=4,type="l"); points(supsmu(1/s2[,2], s2[,4]), lty=lt[3],col=clr[3],lwd=4,type="l"); points(supsmu(1/s2ha[,2], s2ha[,4]), lty=lt[4], col=clr[4],lwd=4,type="l");

plot(supsmu((1/s3[,2]), s3[,3]), col=clr[1], lty=lt[1],ylim=c(0,1),cex.lab=1.3,cex.axis=1.2, ylab="", xlim=c(1,10), xlab=xlab1, lwd=4, type="l");   text(3.9,.95,plotlab[3],cex=1.6); mtext(2,text=ylab1,line=4); mtext(2,text=ylab2,line=3);
points(supsmu((1/s3ha[,2]), s3ha[,3]), lty=lt[2],col=clr[2],lwd=4,type="l"); points(supsmu(1/s3[,2], s3[,4]), lty=lt[3],col=clr[3],lwd=4,type="l"); points(supsmu(1/s3ha[,2], s3ha[,4]), lty=lt[4], col=clr[4],lwd=4,type="l");

plot(supsmu((1/s4[,2]), s4[,3]), col=clr[1], lty=lt[1],ylim=c(0,1),cex.lab=1.3,cex.axis=1.2, ylab="", xlim=c(1,10), xlab=xlab1, lwd=4, type="l");   text(3.2,.95,plotlab[4],cex=1.6); #mtext(2,text=ylab1,line=4); mtext(2,text=ylab2,line=3);
points(supsmu((1/s4ha[,2]), s4ha[,3]), lty=lt[2],col=clr[2],lwd=4,type="l"); points(supsmu(1/s4[,2], s4[,4]), lty=lt[3],col=clr[3],lwd=4,type="l"); points(supsmu(1/s4ha[,2], s4ha[,4]), lty=lt[4], col=clr[4],lwd=4,type="l");
par(mfrow=c(1,1));






#########APPENDIX A
#Figure A1: uses a simplified version of the function mDD (with arbitrary K=1000) to illustrate 4 functional forms density-dependent maturation
min=0.2; K=1000; max=1
mDDsimple=function(x,max,type){  
y=matrix(0,length(x),5);
	for(i in 1:length(x)){ ifelse((x[i]<K),(y[i,1]=max-(((max-min)/K)*x[i])),(y[i,1]=min)) };  ##Type I, Equation 4 in text
	for(i in 1:length(x)){ ifelse((x[i]<K),(y[i,2]=max*K/(K+((max/min)-1)*x[i])),(y[i,2]=min)) };  ##Type II, Equation 6 in Appendix A
	for(i in 1:length(x)){ ifelse((x[i]<K),(y[i,3]=min+(max-min)*sqrt(1-(((1/K)*x[i])^(2*(1-min)))) ),(y[i,3]=min)) };  ##Type III, Equation 7 in Appendix A
	y[,4]=max;
return(as.vector(y[,type])) };



x=seq(1,1500); dev.new();
plot(mDDsimple(x,max,3)~x, ylim=c(0,1), xlim=c(0,1500), type="l", lwd=4, lty=4, cex.lab=1.25,cex.axis=1.35,xlab="Adult abundance",ylab="Maturation rate");
points(mDDsimple(x,max,2)~x,col=2, type="l", lty=2, lwd=4); points(mDDsimple(x,max,1)~x,col=3, lty=3, type="l", lwd=6); m=c(0,1500); n=c(.5,.5); points(n~m,type="l",lwd=4,col=4);
legend("topright", c("Type 0","Type I","Type II","Type III"), col=(4:1), lwd=4, lty=c(1,3,2,4), cex=1.25, seg.len=3, box.col=0);




#Figure A2: similar to the function for Figure 3 above, but with 2 more types of density-dependent maturation functions (and also a constant maturation case for reference)
MvsExtinct2.Sim=function(params){     ntrials=params[1];
harvest = seq(from=0, to=1, length.out=ntrials);
steps=16  #number of scalings done
max = seq(from=min, to=1, length.out=steps);
hold = matrix(0, steps, 16); hold[,1] = max;

const.df=run(c(params, max[1], DDtype=4))
hold[,14]=harvest[which.min(subset(const.df[,3], const.df[,3]>0))];
hold[,15]=harvest[ which.min(abs((const.df[,2]/const.df[1,2])-0.5)) ];
hold[,16]=max(rowSums(const.df[,2:3])/sum(const.df[1,2:3]));

for(y in 1:steps){
	type1.df=run(c(params, max[y], DDtype=1))
	hold[y,2]=max(type1.df[,4]);
	hold[y,3]=harvest[which.min(subset(type1.df[,3], type1.df[,3]>0))] 
	hold[y,4]=harvest[ which.min(abs((type1.df[,2]/type1.df[1,2])-0.5)) ] 
	hold[y,5]=max(rowSums(type1.df[,2:3])/sum(type1.df[1,2:3])) 

	type2.df=run(c(params, max[y], DDtype=2))
	hold[y,6]=max(type2.df[,4]);
	hold[y,7]=harvest[which.min(subset(type2.df[,3], type2.df[,3]>0))] 
	hold[y,8]=harvest[ which.min(abs((type2.df[,2]/type2.df[1,2])-0.5)) ] 
	hold[y,9]=max(rowSums(type2.df[,2:3])/sum(type2.df[1,2:3])) 

	type3.df=run(c(params, max[y], DDtype=3))
	hold[y,10]=max(type3.df[,4]);
	hold[y,11]=harvest[which.min(subset(type3.df[,3], type3.df[,3]>0))]  
	hold[y,12]=harvest[ which.min(abs((type3.df[,2]/type3.df[1,2])-0.5)) ]
	hold[y,13]=max(rowSums(type3.df[,2:3])/sum(type3.df[1,2:3]))
}
Out.df=data.frame(hold)
colnames(Out.df)=c("mmax","I_mmax","I_hext","I_Ahalf","I_Npeak", "II_mmax","II_hext","II_Ahalf","II_Npeak", "III_mmax","III_hext","III_Ahalf","III_Npeak", "C_hext","C_Ahalf","C_Npeak")
return(Out.df) }


min=0.4
params[1:8]=c(30,150,survival1,10,1,1, 1); params[9:11]=c(1,1, 2); s1=MvsExtinct2.Sim(params);
params[1:8]=c(30,150,survival1,10,1,1, 1); params[9:11]=c(0,1, 2); s1ha=MvsExtinct2.Sim(params);

min=0.2;
params[1:8]=c(30,150,survival2,10,1,1, 1); params[9:11]=c(1,1, 2); s2=MvsExtinct2.Sim(params);
params[1:8]=c(30,150,survival3,10,1,1, 1); params[9:11]=c(1,1, 2); s3=MvsExtinct2.Sim(params);
params[1:8]=c(30,150,survival4,10,1,1, 1); params[9:11]=c(1,1, 2); s4=MvsExtinct2.Sim(params);
params[1:8]=c(30,150,survival2,10,1,1, 1); params[9:11]=c(0,1, 2); s2ha=MvsExtinct2.Sim(params);
params[1:8]=c(30,150,survival3,10,1,1, 1); params[9:11]=c(0,1, 2); s3ha=MvsExtinct2.Sim(params);
params[1:8]=c(30,150,survival4,10,1,1, 1); params[9:11]=c(0,1, 2); s4ha=MvsExtinct2.Sim(params);


#plotting; survival3 (sj=0.8, sa=0.2) is very similar to survival4 (sj=sa=0.8) and commented out here for clarity
dev.new(width=14.8,height=6.3); par(mfrow=c(1,2),mai=c(1,1,1,1)); ymax=0.85;
plot(supsmu(1/s1[,2], s1[,4]), type="l",pch=1,cex=1,lwd=4,lty=3,cex.lab=1.3,cex.axis=1.2,ylim=c(0,ymax), ylab=expression("Harvest level at 50% adult decline (h"["   half"]~")"), xlim=c(1,5), xlab=expression(Minimum~age~at~maturity~(~italic(m[max])~""^{-1})));
legend("topleft", box.col=0,c("Type II","Type III"), pch=c(1,18),cex=1.8); text(3.2,0.02,expression("s"[J]~"="~"s"[A]~"="~"0.2"),cex=1.35);  text(4.4,0.4,expression("s"[J]~"="~"s"[A]~"="~"0.8"),cex=1.3);  text(4.18,0.18,expression("s"[J]~"="~"0.2,  s"[A]~"="~"0.8"),cex=1.25);
points(supsmu(1/s1[,2], s1[,8]), col=1,lwd=2,cex=1.25,type="p",pch=1); points(supsmu(1/s1[,2], s1[,12]), col=1,lwd=3,type="p",pch=18,cex=1.5);
points(supsmu(1/s4[,2], s4[,4]), col=3,lwd=4,type="l"); points(supsmu(1/s4[,2], s4[,8]), col=3,lwd=2,cex=1.25,type="p",pch=1); points(supsmu(1/s4[,2], s4[,12]), col=3,lwd=3,type="p",pch=18,cex=1.5); 
points(supsmu(1/s2[,2], s2[,4]), col=2,lwd=4,type="l",lty=2); points(supsmu(1/s2[,2], s2[,8]), col=2,lwd=2, cex=1.25,type="p",pch=1); points(supsmu(1/s2[,2], s2[,12]), col=2,lwd=3,type="p",pch=18,cex=1.5); 
#points(supsmu(1/s3[,2], s3[,4]), col=4,lwd=4,type="l",lty=4); points(supsmu(1/s3[,2], s3[,8]), col=4,lwd=3,type="p",pch=15); points(supsmu(1/s3[,2], s3[,12]), col=4,lwd=3,type="p",pch=18,cex=1.5); 

plot(supsmu(1/s1ha[,2], s1ha[,4]), type="l",pch=1,cex=1,lwd=4,lty=3,cex.lab=1.3,cex.axis=1.2,ylim=c(0,ymax), ylab=expression("Harvest level at 50% adult decline (h"["   half"]~")"), xlim=c(1,5), xlab=expression(Minimum~age~at~maturity~(~italic(m[max])~""^{-1})));
text(3.2,0.02,expression("s"[J]~"="~"s"[A]~"="~"0.2"),cex=1.35);  text(4,0.7,expression("s"[J]~"="~"s"[A]~"="~"0.8"),cex=1.35);  text(4.2,0.25,expression("s"[J]~"="~"0.2,  s"[A]~"="~"0.8"),cex=1.25);
points(supsmu(1/s1ha[,2], s1ha[,8]), col=1,cex=1.25,type="p",pch=1,lwd=2.5); points(supsmu(1/s1ha[,2], s1ha[,12]), col=1,type="p",pch=18,cex=1.5);
points(supsmu(1/s2ha[,2], s2ha[,4]), col=2,lwd=4,type="l",lty=2); points(supsmu(1/s2ha[,2], s2ha[,8]), col=2,cex=1.25,type="p",pch=1,lwd=2); points(supsmu(1/s2ha[,2], s2ha[,12]), col=2,type="p",pch=18,cex=1.5); 
points(supsmu(1/s4ha[,2], s4ha[,4]), col=3,lwd=4,type="l"); points(supsmu(1/s4ha[,2], s4ha[,8]), col=3,lwd=2,cex=1.25,type="p",pch=1); points(supsmu(1/s4ha[,2], s4ha[,12]), col=3,lwd=3,type="p",pch=18,cex=1.5); 
#points(supsmu(1/s3ha[,2], s3ha[,4]), col=4,lwd=4,type="l",lty=4); points(supsmu(1/s3ha[,2], s3ha[,8]), col=4,lwd=3,type="p",pch=15); points(supsmu(1/s3ha[,2], s3ha[,12]), col=4,lwd=3,type="p",pch=18,cex=1.5);






#########APPENDIX B
#Figure B1 - NOTE: must have run function mDD2.Sim earlier
min=0.2
params[1:8]=c(50,150,survival1,10,1,1, 1);  params[9:11]=c(1,1, 2); s1=mDD2.Sim(params);
params[1:8]=c(50,150,survival2,10,1,1, 1);  params[9:11]=c(1,1, 2); s2=mDD2.Sim(params);
params[1:8]=c(50,150,survival3,10,1,1, 1);  params[9:11]=c(1,1, 2); s3=mDD2.Sim(params);
params[1:8]=c(50,150,survival4,10,1,1, 1);  params[9:11]=c(1,1, 2); s4=mDD2.Sim(params);

params[1:8]=c(50,150,survival1,10,1,1, 2);  params[9:11]=c(1,1, 2); s1bh=mDD2.Sim(params);
params[1:8]=c(50,150,survival2,10,1,1, 2);  params[9:11]=c(1,1, 2); s2bh=mDD2.Sim(params);
params[1:8]=c(50,150,survival3,10,1,1, 2);  params[9:11]=c(1,1, 2); s3bh=mDD2.Sim(params);
params[1:8]=c(50,150,survival4,10,1,1, 2);  params[9:11]=c(1,1, 2); s4bh=mDD2.Sim(params);


dev.new(); par(mfrow=c(1,2)); ymax=1.3; xmax=0.8;
plot((s1[,11]+s1[,12])~s1[,1],    type="b",col=colors()[228],lwd=3,lty=1,pch=1, cex=1.25, xlim=c(0,xmax), ylim=c(0,ymax),ylab="Equilibrium population size",xlab="Proportion of total population harvested",cex.lab=1.35); text(0.05,1.9, "A:",cex=2);
points((s2[,11]+s2[,12])~s2[,1],type="b",col=colors()[203],lwd=3,lty=1,pch=1, cex=1.25); points((s3[,11]+s3[,12])~s3[,1],type="b",col=colors()[178],lwd=3,lty=1,pch=1, cex=1.25); points((s4[,11]+s4[,12])~s4[,1],type="b",col=colors()[153],lwd=3,lty=1,pch=1, cex=1.25);

plot((s1bh[,11]+s1bh[,12])~s1bh[,1],    type="b",col=colors()[228],lwd=3,lty=1,pch=1, cex=1.25, xlim=c(0,xmax), ylim=c(0,ymax),ylab="Equilibrium population size",xlab="Proportion of total population harvested",cex.lab=1.35); text(0.05,1.9, "A:",cex=2);
points((s2bh[,11]+s2bh[,12])~s2bh[,1],type="b",col=colors()[203],lwd=3,lty=1,pch=1, cex=1.25); points((s3bh[,11]+s3bh[,12])~s3bh[,1],type="b",col=colors()[178],lwd=3,lty=1,pch=1, cex=1.25); points((s4bh[,11]+s4bh[,12])~s4bh[,1],type="b",col=colors()[153],lwd=3,lty=1,pch=1, cex=1.25);
legend("topright", seg.len=2.5, lty=c(1), pch=c(1), lwd=3, pt.cex=1.5,  cex=1.25,box.col=0,rev(c(expression("s"[J]~"="~"s"[A]~"="~"0.2"), expression("s"[J]~"="~"0.2,  s"[A]~"="~"0.8"), expression("s"[J]~"="~"0.8,  s"[A]~"="~"0.2"), expression("s"[J]~"="~"s"[A]~"="~"0.8"))), col=rev(colors()[c(228,203,178,153)]));





#Figure B2 - NOTE: must have run function MvsExtinct.Sim earlier
min=0.4
params[1:11]=c(40,300,survival1,10,1,1,1, 1,2, 2); s1=MvsExtinct2.Sim(params);
params[1:11]=c(40,300,survival1,10,1,1,1, 1,1, 2); s1ha=MvsExtinct2.Sim(params);

min=0.2
params[1:8]=c(40,300,survival2,10,1,1,2); params[9:11]=c(1,1, 2); s2=MvsExtinct2.Sim(params);
params[1:8]=c(40,300,survival3,10,1,1,2); params[9:11]=c(1,1, 2); s3=MvsExtinct2.Sim(params);
params[1:8]=c(40,300,survival4,10,1,1,2); params[9:11]=c(1,1, 2); s4=MvsExtinct2.Sim(params);

params[1:8]=c(40,300,survival2,10,1,1,1); params[9:11]=c(1,1, 2); s2ha=MvsExtinct2.Sim(params);
params[1:8]=c(40,300,survival3,10,1,1,1); params[9:11]=c(1,1, 2); s3ha=MvsExtinct2.Sim(params);
params[1:8]=c(40,300,survival4,10,1,1,1); params[9:11]=c(1,1, 2); s4ha=MvsExtinct2.Sim(params);


dev.new(); par(mfrow=c(1,2),mai=c(1,1,1,1)); ymax=0.6;
plot(supsmu(1/s1ha[,2], s1ha[,4]), type="l",cex=1,col=8,lwd=4,lty=3,cex.lab=1.3,cex.axis=1.2,ylim=c(0,ymax), ylab=expression("Harvest level at 50% adult decline (h"["   half"]~")"), xlim=c(1,5), xlab=expression(Minimum~age~at~maturity~(~italic(m[max])~""^{-1})));
legend("topright", box.col=0,c("Type II","Type III"), pch=c(1,18),cex=1.6); text(3.2,0.02,expression("s"[J]~"="~"s"[A]~"="~"0.2"),cex=1.35);  text(4,0.4,expression("s"[J]~"="~"s"[A]~"="~"0.8"),cex=1.35);  text(4.2,0.18,expression("s"[J]~"="~"0.2,  s"[A]~"="~"0.8"),cex=1.25);
points(supsmu(1/s1ha[,2], s1ha[,8]), col=8,cex=1.25,type="p",pch=1,lwd=2.5); points(supsmu(1/s1ha[,2], s1ha[,12]), col=8,type="p",pch=18,cex=1.5);
points(supsmu(1/s2ha[,2], s2ha[,4]), col=8,lwd=4,type="l",lty=2); points(supsmu(1/s2ha[,2], s2ha[,8]), col=8,cex=1.25,type="p",pch=1,lwd=2); points(supsmu(1/s2ha[,2], s2ha[,12]), col=8,type="p",pch=18,cex=1.5); 
points(supsmu(1/s4ha[,2], s4ha[,4]), col=1,lwd=4,type="l"); points(supsmu(1/s4ha[,2], s4ha[,8]), col=1,lwd=2,cex=1.25,type="p",pch=1); points(supsmu(1/s4ha[,2], s4ha[,12]), col=1,lwd=3,type="p",pch=18,cex=1.5); 

plot(supsmu(1/s1[,2], s1[,4]), type="l",col=8,cex=1,lwd=4,lty=3,cex.lab=1.3,cex.axis=1.2,ylim=c(0,.4), ylab=expression("Harvest level at 50% adult decline (h"["   half"]~")"), xlim=c(1,5), xlab=expression(Minimum~age~at~maturity~(~italic(m[max])~""^{-1})));
points(supsmu(1/s1[,2], s1[,8]), col=8,lwd=2,cex=1.25,type="p",pch=1); points(supsmu(1/s1[,2], s1[,12]), col=8,lwd=3,type="p",pch=18,cex=1.5);
points(supsmu(1/s2[,2], s2[,4]), col=8,lwd=4,type="l",lty=2); points(supsmu(1/s2[,2], s2[,8]), col=8,lwd=2, cex=1.25,type="p",pch=1); points(supsmu(1/s2[,2], s2[,12]), col=8,lwd=3,type="p",pch=18,cex=1.5);  
points(supsmu(1/s4[,2], s4[,4]), col=1,lwd=4,type="l"); points(supsmu(1/s4[,2], s4[,8]), col=1,lwd=2,cex=1.25,type="p",pch=1); points(supsmu(1/s4[,2], s4[,12]), col=1,lwd=3,type="p",pch=18,cex=1.5); 
text(3.2,0.02,expression("s"[J]~"="~"s"[A]~"="~"0.2"),cex=1.35);  text(4.4,0.16,expression("s"[J]~"="~"s"[A]~"="~"0.8"),cex=1.3);  text(4.18,0.07,expression("s"[J]~"="~"0.2,  s"[A]~"="~"0.8"),cex=1.25);
par(mfrow=c(1,1));








################Appendix C
#For this part we used a simplified version of mDD (need to run mDDsimple from Appendix A beforehand) and a slightly modified version of the run function (both changes noted in comments)
rm(list = ls())
min=0.5; K=1000; max=1;
mDDsimple=function(x,max,type){  
y=matrix(0,length(x),5);
	for(i in 1:length(x)){ ifelse((x[i]<K),(y[i,1]=max-(((max-min)/K)*x[i])),(y[i,1]=min)) };  ##Type I, Equation 4 in text
	for(i in 1:length(x)){ ifelse((x[i]<K),(y[i,2]=max*K/(K+((max/min)-1)*x[i])),(y[i,2]=min)) };  ##Type II, Equation 6 in Appendix A
	for(i in 1:length(x)){ ifelse((x[i]<K),(y[i,3]=min+(max-min)*sqrt(1-(((1/K)*x[i])^(2*(1-min)))) ),(y[i,3]=min)) };  ##Type III, Equation 7 in Appendix A
	y[,4]=max;
return(as.vector(y[,type])) };




run2=function(params){
ntrials=params[1]; runlength=params[2]; sj=params[3]; sa=params[4]; alpha=params[5]; beta=params[6:7]; recruitType=params[8]; harvMod=params[9:10]; mDDstage=params[11]; ntimes=params[12]; varmod=params[13]; phase=params[14];
max=params[15]; DDtype=params[16];
harvest = seq(from=0, to=1, length.out=ntrials)
hold = matrix(0,ntrials,4); hold[,1] = harvest;

for(x in 1:ntrials){
		harvest=c(hold[x,1],hold[x,1])*harvMod #Set harvest level for simulation; adjust for stage-specific harvest (if present). hj=harvest[1], ha=harvest[2]
		J=double(runlength); J[1]=1500;
		A=double(runlength); A[1]=1000;
		mrec=double(runlength);
			for(t in 2:runlength){
				m=as.numeric( c(mDDsimple(J[t-1],max,DDtype), mDDsimple(A[t-1],max,DDtype))[mDDstage] )    ## <-- Uses mDDsimple
				alpha=params[5]*(1-params[7]*((m-min)/(max-min)));                                         ## <-- FOR m-d alpha vs h_c
				gamma=c( alpha*exp(-beta[1]*A[t-1]*(1-harvest[2])), alpha/(1+beta[2]*A[t-1]*(1-harvest[2])) )[recruitType]  #Note: in all simulations, both beta values =1
				Juves = A[t-1]*(1-harvest[2])*gamma + J[t-1]*(1-harvest[1])*sj*(1-m)
				Adults = J[t-1]*(1-harvest[1])*sj*m + A[t-1]*(1-harvest[2])*sa
				J[t]=ifelse(Juves>0.001,Juves,0); A[t]=ifelse(Adults>0.001,Adults,0); mrec[t] = m;
			}

hold[x,2]=mean(A[(runlength-phase):runlength]); hold[x,3]=mean(J[(runlength-phase):runlength]); 
hold[x,4]=mean(mrec[(runlength-phase):(runlength-1)]);   #mean maturation rate
}
return(hold);  }





MdAlphaVsExtinct.Sim=function(params){     ntrials=params[1];
harvest = seq(from=0, to=1, length.out=ntrials);
steps=16  #number of scalings done
varmod = seq(from=0, to=0.75, length.out=steps);
hold = matrix(0, steps, 10); hold[,1] = varmod;

params[7]=0; const.df=run2(c(params, (min+0.01), DDtype=4))
hold[,5]=harvest[which.min(subset(const.df[,3], const.df[,3]>0))];
hold[,6]=harvest[ which.min(abs((const.df[,2]/const.df[1,2])-0.5)) ];
hold[,7]=max(rowSums(const.df[,2:3])/sum(const.df[1,2:3]));

for(y in 1:steps){
	params[7]=varmod[y]; type1.df=run2(c(params, 1, DDtype=1))
	hold[y,2]=harvest[which.min(subset(type1.df[,3], type1.df[,3]>0))]    /    hold[1,5];
	hold[y,3]=harvest[ which.min(abs((type1.df[,2]/type1.df[1,2])-0.5)) ]    /    hold[1,6];
	hold[y,4]=max(rowSums(type1.df[,2:3])/sum(type1.df[1,2:3]))    /    hold[1,7];

	hold[y,8]=harvest[which.min(subset(type1.df[,3], type1.df[,3]>0))];
	hold[y,9]=harvest[ which.min(abs((type1.df[,2]/type1.df[1,2])-0.5)) ];
	hold[y,10]=max(rowSums(type1.df[,2:3])/sum(type1.df[1,2:3]));
}

Out.df=data.frame(hold)
colnames(Out.df)=c("da/dm", "I_hextRelC","I_AhalfRelC","I_NpeakRelC", "C_hext","C_Ahalf","C_Npeak", "I_hext","I_Ahalf","I_Npeak")
return(Out.df) }








survival1=c(0.2,0.2); survival2=c(0.2,0.8); survival3=c(0.8,0.2); survival4=c(0.8,0.8);
#ntrials, runlength, [sj,sa], alpha, betaOCdd,  betaBHdd, DD_Type(1=OC|2=BH), harvMod(J|A), mDDstage(J|A), ntrials, SurvVarMod, phase(def=3)
params=c(30, 300,   survival1, 10,  0.0031,    0.00064,          1,            c(1,1),         2,            1,      0.,        60) 
params[1:8]=c(55,300,survival1,10,0.00044,0,1); s1=MdAlphaVsExtinct.Sim(params);
params[1:8]=c(55,300,survival4,10,0.0031,0,1); s4=MdAlphaVsExtinct.Sim(params);

dev.new()
plot(supsmu(s1[,1],s1[,8]), ylim=c(0,1),xlim=c(0,0.8),xlab=expression("Proportional decrease in"~alpha~"at m"[max]~(sigma)),ylab="Harvest intensity",col=1,type="l",lwd=3.5,cex.lab=1.25); points(supsmu(s1[,1],s1[,9]), col=1,type="p",pch=16); 
points(supsmu(s4[,1],s4[,8]), col=3,type="l",lwd=3.5); points(supsmu(s4[,1],s4[,9]), col=3,type="p",pch=16);
points(s1[,1],s1[,5], col=1,type="l",lty=2,lwd=2); points(s1[,1],s1[,6], col=1,type="l",lty=3,lwd=2); 
points(s4[,1],s4[,5], col=3,type="l",lty=2,lwd=2); points(s4[,1],s4[,6], col=3,type="l",lty=3,lwd=2);
legend("topright",box.col=0, c(expression("Type I  h"[c]), expression("Type I  h"[half])), col=1, lwd=3.5, lty=c(1,NA), pch=c(NA,16));
legend("top",box.col=0, c(expression("Type 0  h"[c]), expression("Type 0  h"[half])), col=1, lwd=2.5, lty=c(2,3));
text(0.6,0.7,expression("s"[J]~"="~"s"[A]~"="~"0.8"),cex=1.25,col=3); text(0.6,0.28,expression("s"[J]~"="~"s"[A]~"="~"0.2"),cex=1.25,col=1);