#This script simulates three models of the relationship between species' mean frequencies and the strength of negative frequency dependence (NFD) in their per-capita growth rates
#One model is a multispecies version of the annual plant competition model of Yenni et al. 2012
#This model incorporates intra- and interspecific competition and demographic stochasticity
#The primary question of interest is whether, for randomly-generated parameter sets, rare species tend to experience
#stronger negative frequency dependence (NFD) than common ones, as hypothesized by Yenni et al. 2012 for a 2-species
#version of the model and as suggested by the empirical data in Yenni et al. 2017
#The second model is a discrete time Lotka-Volterra competition model incorporating
#environmental rather than demographic stochasticity, and seasonality in the form of periodic density-independent morality ("winter")
#The third model considers purely mathematical and ecological constraints. Given that there is some finite
#maximum possible per-capita growth rate when rare,
#and that species' mean frequencies must sum to 1, what is the expected distribution of relationships between
#NFD and mean frequency?
#This script also considers species richness and evenness as mediators of the rarity-NFD relationship.
#************************************
#load package to calculate Gini coefficients as a measure of evennness
library(digest)
library(stringi)
library(reldist)
#***************************
#LOTKA-VOLTERRA MODEL WITH ENVIRONMENTAL STOCHASTICITY AND WINTER MORTALITY
rm(list=ls(all=TRUE)) #clear the decks
runs<-100 #number of simulations to run
tmax<-1000 #maximum number of timesteps
results<-matrix(-1,nrow=runs,ncol=11) #matrix to hold the results
for(x in 1:runs){
#NOTE: the outer loop over all the simulations may terminate if a simulation results in competitive exclusion of all but one species
#If the loop over all simulations terminates, just store the results of the simulations that ran, and then run additional simulations until the desired total number of simulation runs is achieved
S<-20 #initial species richness. The higher you set this, the more species will persist on average, up to some limit that depends on the other parameters
season.length<-10 #number of timesteps between each winter mortality event
winter.survival<-matrix(-1,nrow=tmax,ncol=S) #matrix to hold winter survival probabilities. Generating tmax winter survival probabilities for each species isn't actually necessary bc there aren't tmax winters, but it doesn't cause any problems.
for(i in 1:S){
winter.survival[,i]<-runif(tmax,min=0.5,max=0.8) #randomly generate sequence of winter survival probabilities for each species.
}
r<-matrix(-1,nrow=tmax,ncol=S) #matrix to hold species' intrinsic rates of increase at each timestep
for(i in 1:S){
r[,i]<-runif(tmax,min=0.5,max=2) #randomly generate sequence of r values for each species. They fluctuate over time due to envi. stochas.
}
a.ii<-runif(S,min=0.001,max=0.01) #randomly choose species' intraspecific competition coefficients
a<-matrix(data=runif(S*(S-1),min=0,max=0.003),nrow=S,ncol=S) #create matrix of intersp comp coeffs randomly chosen from a uniform distn
diag(a)<-a.ii #set the diagonal elements of the comp coef matrix equal to the intraspecific comp coefs
N0<-matrix(10,nrow=1,ncol=S) #initial abundances
N<-matrix(-1,nrow=tmax,ncol=S) #matrix to hold species abundances over time
N[1,]<-N0
for(t in 2:tmax){
for(i in 1:S){
N[t,i]<-ifelse(N[t-1,i]<0.00001,0,r[t,i]*N[t-1,i]-N[t-1,i]*sum(a[,i]*N[t-1,])) #Lotka-Volterra competition
N[t,i]<-ifelse(N[t,i]<0.00001,0,N[t,i]) #set abundances below 10^-5 to zero
N[t,i]<-ifelse(t%%season.length==0,N[t,i]*winter.survival,N[t,i]) #winter mortality
N[t,i]<-ifelse(N[t,i]<0.00001,0,N[t,i]) #set abundances below 10^-5 to zero
}
}
N.persist<-tail(N,500) #remove all but the last 500 timesteps
N.persist[N.persist==0]<-NA
N.persist <- N.persist[ , colSums(is.na(N.persist)) == 0] #drop species that went extinct
S.persist<-ncol(N.persist) #final species richness
sub1<-N.persist[-nrow(N.persist),]
sub2<-N.persist[-1,]
pcgr<-log(sub2/sub1) #per-capita growth rates
freq<-N.persist/rowSums(N.persist) #note that this is a fast but slightly risky way to calculate species' frequencies. See https://rpubs.com/bbolker/sweep_divide
freq<-head(freq,-1)
results_sim <- array(NA, dim = c(S.persist,2))
for(i in 1:S.persist){ #loop over persisting species to calculate strength of NFD
fit <- lm(pcgr[,i]~freq[,i])
slope <- (coef(fit)[[2]])
mean.freq <- mean(freq[,i])
results_sim[i,1] <- slope #Enters slope. Negative slope = NFD
results_sim[i,2] <- mean.freq #Enters mean frequency
}
FD <- as.numeric(results_sim[,1]) #strength of NFD
mean.frequency <- as.numeric(results_sim[,2]) #species' mean frequencies
plot(mean.frequency, FD, main="NFD vs. mean frequency",xlab="Mean frequency",ylab="NFD")
evenness<-gini(mean.frequency)
cov_sim <- cov(mean.frequency,FD)
cov_sim
results[x,]<-c(S,mean(r),sd(r),mean(a.ii),sd(a.ii),(sum(a)-sum(a.ii))/(S*(S-1)),max(mean.frequency)-min(mean.frequency),S.persist,mean(FD),cov_sim,evenness) #save key summary measures of each simulation run
}
results<-results[results[,1]>-1,] #drop rows of results corresponding to simulations that failed or weren't conducted at all
results
#****************************
#ANNUAL PLANT MODEL WITH DEMOGRAPHIC STOCHASTICITY
rm(list=ls(all=TRUE)) #clear the decks
runs<-100 #number of simulations to run
results<-matrix(-1,nrow=runs,ncol=11) #matrix to hold the results
for(x in 1:runs){
#NOTE: If you get the error "Error in -nrow(N.persist) : invalid argument to unitary operator", it probably means that all species but one went extinct in one of the simulations, which causes R to not do any further simulations.
#If the loop over all simulations terminates, just store the results of the simulations that ran, and then run additional simulations until the desired total number of simulation runs is achieved
S<-10 #initial species richness. The higher you set this, the more species will persist on average, up to some limit that depends on the other parameters
r<-runif(S,min=11,max=20) #randomly choose species' per-capita fecundities from a uniform distn
a.ii<-runif(S,min=0.7,max=3) #randomly choose species' intraspecific competition coefficients
a<-matrix(data=runif(S*(S-1),min=0,max=0.1),nrow=S,ncol=S) #create matrix of intersp comp coeffs randomly chosen from a uniform distn
diag(a)<-a.ii #set the diagonal elements of the comp coef matrix equal to the intraspecific comp coefs
N0<-matrix(50,nrow=1,ncol=S) #initial abundances
tmax<-500 #maximum number of timesteps
N<-matrix(-1,nrow=tmax,ncol=S) #matrix to hold species abundances over time
N[1,]<-N0 #initial abundances
for(t in 2:tmax){
for(i in 1:S){
N[t,i]<-ifelse(N[t-1,i]==0,0,r[i]*N[t-1,i]/(1+sum(a[,i]*N[t-1,])))
ifelse(N[t,i]<0,0,N[t,i]) #set negative abundances to zero
N[t,i]<-ifelse(N[t,i]==0,0,rpois(1,N[t,i])) #incorporate demographic stochasticity. Deterministic model defines the mean of a Poisson distn. of # of offspring in following year.
}
}
N.persist<-tail(N,100) #remove all but the last 100 timesteps
N.persist[N.persist==0]<-NA
N.persist <- N.persist[ , colSums(is.na(N.persist)) == 0] #drop species that went extinct
S.persist<-ncol(N.persist) #final species richness
sub1<-N.persist[-nrow(N.persist),]
sub2<-N.persist[-1,]
pcgr<-log(sub2/sub1) #per-capita growth rates
freq<-N.persist/rowSums(N.persist) #note that this is a fast but slightly risky way to calculate species' frequencies. See https://rpubs.com/bbolker/sweep_divide
freq<-head(freq,-1)
results_sim <- array(NA, dim = c(S.persist,2))
for(i in 1:S.persist){ #loop over persisting species to calculate strength of NFD
fit <- lm(pcgr[,i]~freq[,i])
slope <- (coef(fit)[[2]])
mean.freq <- mean(freq[,i])
results_sim[i,1] <- slope #Enters slope. Negative slope = NFD
results_sim[i,2] <- mean.freq #Enters mean frequency
}
FD <- as.numeric(results_sim[,1]) #strength of NFD
mean.frequency <- as.numeric(results_sim[,2]) #species' mean frequencies
evenness<-gini(mean.frequency)
cov_sim <- cov(mean.frequency,FD)
cov_sim
results[x,]<-c(S,mean(r),sd(r),mean(a.ii),sd(a.ii),(sum(a)-sum(a.ii))/(S*(S-1)),max(mean.frequency)-min(mean.frequency),S.persist,mean(FD),cov_sim,evenness) #save key summary measures of each simulation run
}
results<-results[results[,1]>-1,] #when <2 species persist until tmax, the code above fails. So we need to delete those simulation results. They're rows of -1.
results
#******************************************
#CONSTRAINTS-BASED MODEL
rm(list=ls(all=TRUE)) #clear the decks
runs<-1000 #number of simulations to run per value of species richness
Smax<-10 #maximum species richness to consider
results<-matrix(0,nrow=runs,ncol=Smax-1) #matrix to hold NFD-frequency covariance for each simulation; number of columns equals number of species richness values to try
richness<-matrix(0,nrow=runs,ncol=Smax-1) #matrix to hold species richness values
gini.coeff<-matrix(0,nrow=runs,ncol=Smax-1) #matrix to hold Gini coefficients, a measure of evenness
for(i in 1:runs){
for(j in 1:Smax-1){
S<-j+1 #species richness
rmax<-1 #maximum possible r value
N<-runif(S,min=0,max=1) #species' absolute abundances
N<-N/sum(N) #species' frequencies
r<-runif(S,min=0,max=rmax) #r values (per-capita growth rate at ~0 freq.)
NFD<--r/N #implied slope of relationship between per-capita growth rate and frequency for each species
covariance<--cov(NFD,N)
results[i,j]<-covariance
richness[i,j]<-j+1
gini.coeff[i,j]<-gini(t(N))
}
}