R code accompanying Vellend (2016), ‘The Theory of Ecological Communities’ (Princeton University Press)

## specify parameters, initial conditions, and output vector
num.years <- 10000
J <- 10000 
sp <- 1 # a counter to keep track of added species
COM <- vector(length = J); COM[] <- sp 

nu <- 5e-04 # speciation rate

year <- 2 

COM.out <- matrix(nrow = (num.years/100+1), ncol = J); COM.out[1,] <- COM 

## run simulation
for (i in 1:(J*(num.years))) {
  dead.indiv <- ceiling(J*runif(1))
  
  ## speciation (replace dead individual with new species)
  if (runif(1) < nu) {
    COM[dead.indiv] <- sp+1
    sp <- sp + 1
  } else {

  ## not speciation (replace dead individual from community)
    COM[dead.indiv] <- COM[ceiling(J*runif(1))]
  }
  
  ## record data every 100 years
  if (i %% (J*100) == 0) {
    COM.out[year,] <- COM
    year <- year+ 1   
  } } 

################################
# Calculations to allow plotting relative abundance distributions

## First calculate the number of species (SR) at each time

SR.out <- vector(length = dim(COM.out)[1])

for (t in 1:dim(COM.out)[1]) {
  COM.hist <- hist(COM.out[t,],c(1:max(COM.out[t,])),plot=FALSE)
  abund.vec <- COM.hist$counts[COM.hist$counts>0]
  SR.out[t] <- sum(abund.vec>0)
}

## Calculate species abundances
abund.out <- matrix(nrow = (num.years/100+1), ncol = max(SR.out))
abund.out[] <- 0

for (t in 1:dim(COM.out)[1]) {
  COM.hist <- hist(COM.out[t,],c(1:max(COM.out[t,])),plot=FALSE)
  abund.vec <- COM.hist$counts[COM.hist$counts>0]
  abund.out[t,1:SR.out[t]] <- abund.vec
}

abund.out <- abund.out/J

## The final rank-order of relative abundances
final.rad <- sort(abund.out[(num.years/100+1),1:SR.out[(num.years/100+1)]],decreasing=TRUE)

## Write the final relative abundance distribution to a file for use in Online Box 8 (island biogeography model)
write.csv(final.rad, file = "RAD1.csv")

## graph the results
plot(log(final.rad),type="l",
     xlab="Species rank order", ylab="log(Relative abundance)")

Description of the code.

1. num.years, J, and COM are defined as in previous Boxes. nu is the speciation rate. With probability nu, the dead individual will be replaced by a species new to the metacommunity. Since we will have an undetermined number of species, sp will keep track of species identities. The first species is species 1, and each new species will be numbered in sequence.

num.years <- 100
J <- 100 
sp <- 1 # a counter to keep track of added species
COM <- vector(length = J); COM[] <- sp 

nu <- 1e-04 # speciation rate

2. With many different species, we can no longer just keep track of the frequency of species 1. ncol = J because we will keep track of the identity of each individual in each year of recorded data. In order to reduce the duration of the simulation, we will record data only every 100 years, which is why nrow = (num.years/100+1)

COM.out <- matrix(nrow = (num.years/100+1), ncol = J); COM.out[1,] <- COM 

3. Since we don’t have just one line in the simulation where COM can be updated, we select an individual to die before figuring out if it will be replaced by a speciation even or by a local reproducer.

  dead.indiv <- ceiling(J*runif(1))

4. Here’s where speciation happens. If the new individual is produced by speciation (with probaiblity nu), we update the sp counter, keeping track of how many species have entered the community.

  if (runif(1) < nu) {
    COM[dead.indiv] <- sp+1
    sp <- sp + 1
  } else {

5. Stop and record data every 100 years (after J*100 deaths) instead of every year (after J deaths) as in previous Boxes).

  if (i %% (J*100) == 0) {
    COM.out[year,] <- COM
    year <- year+ 1   
  } } 

6. We need to calculate the number of species present (Species Richness) at each time step (when data were recorded) in order to avoid outputting zero abundances later. dim(COM.OUT)[1] returns the first dimension (number of rows) of COM.OUT, which is the number of time steps at which we recorded data. The hist (histogram) function returns, in the $counts object, a vector of abundances of each species in the community, most of which will be zero (because most species entering via speciation go quickly extinct). We specify plot=FALSE so that R doesn’t produce a histogram graph each time. abund.vec will contain only the non-zero elements, and abund.vec>0 returns a vector with a ‘1’ for each species with non-zero abundance, the sum of which is species richness.

SR.out <- vector(length = dim(COM.out)[1])

for (t in 1:dim(COM.out)[1]) {
  COM.hist <- hist(COM.out[t,],c(1:max(COM.out[t,])),plot=FALSE)
  abund.vec <- COM.hist$counts[COM.hist$counts>0]
  SR.out[t] <- sum(abund.vec>0)
}

7. Using the same functions as in point 6 above, we calculate a matrix of species relative abundances at teach time point. This is done separately from (and after) the species richness calculations because we don’t need an abundance matrix with any more columns than there were species present at at least one of the data collection time points. Many species arise and disappear within a 100-year window of time. Dividing abundances by J returns relative abundances that sum to one.

abund.out <- matrix(nrow = (num.years/100+1), ncol = max(SR.out))
abund.out[] <- 0

for (t in 1:dim(COM.out)[1]) {
  COM.hist <- hist(COM.out[t,],c(1:max(COM.out[t,])),plot=FALSE)
  abund.vec <- COM.hist$counts[COM.hist$counts>0]
  abund.out[t,1:SR.out[t]] <- abund.vec
}

abund.out <- abund.out/J

8. final.rad is the rank-order relative abundance vector at the final time step. This needs to be saved if we want to use it as the metacommunity relative abundance distribution in Online Box 8. Note that you will need to create/set the directory on your own (i.e., replaced c:/temp/RAD1.csv with whatever directory and filename you want).

final.rad <- sort(abund.out[(num.years/100+1),1:SR.out[(num.years/100+1)]],decreasing=TRUE)

## Write the final relative abundance distribution to a file for use in Online Box 8 (island biogeography model)
write.csv(final.rad, file = "RAD1.csv")