Data and tools

We want to build a species-rich metacommunity and a local, small community that is a sample of the metacommunity, according to the neutral model.

Species

We will use letters to build the species names, e.g. “AAC”. Combining 3 letters, we’ll have \(26^3 = 17576\) species names, i.e. more than the number of species in Amazonia (ter Steege et al. 2013).

To build the names, we use this code:

library("tidyverse")
# LETTERS is the vector of capital letters A:Z
# outer() combines the first two args by FUN
LETTERS %>% 
  outer(LETTERS, FUN = paste0) %>% 
  # add a third letter
  outer(LETTERS, FUN = paste0) ->
  species_names

Metacomunity

We’ll draw random log-series distributions.

#' Random draw of a log-series distribution
#' The complete set of functions (including density, distribution function and quantiles)
#' can be found in package _sads_ but this implementation of the random generation is much faster.
#'
#' @param n The Number of species
#' @param size The size of the distribution.
#' @param alpha Fisher's alpha, that links the number of species to that of individuals
#' @param show_progress Show a progress bar during long computations
#'
#' @return The number of individuals per species
rlseries <- function(
    n, 
    size, 
    alpha, 
    show_progress = TRUE) {
  # adapted from Dan Lunn, http://www.stats.ox.ac.uk/~dlunn/BS1_05/BS1_Rcode.pdf
  if (size + 1 == size || size - 1 == size) {
    # Too large values to deal with as integers
    stop("size is too large to simulate the distribution.")
  }
  # Fisher's x
  x <- size / (size + alpha)
  # Draw random numbers between 0 and 1 that are quantiles of the cumulative function
  u <- sort(stats::runif(n))
  # Prepare the corresponding number of abundances
  abd <- numeric(n)
  # Prepare a progress bar
  if (show_progress & interactive()) {
    cli::cli_progress_bar(
      total = n,
      format = "{cli::pb_spin} Drawing log-series values {cli::pb_current}/{cli::pb_total}"
    )
  }
  # Number of drawn values +1
  next_value <- 1
  # k is the abundance
  k <- 1
  # Calculate the probability at k=1
  P <- -x / log(1 - x)
  # Store it in the cumulative function
  F <- P
  # Repeat while all values are not drawn
  while (next_value <= n) {
    if (F > u[next_value]) {
      # Retain k as the next value
      abd[next_value] <- k
      # Increment the next value
      next_value <- next_value + 1
      # Update the progress bar
      if (show_progress & interactive()) cli::cli_progress_update()
    } else {
      # Probability at k+1 obtained from that at k
      P <- P * k * x / (k + 1)
      # Increment the cumulated probability
      F <- F + P
      # Increment k
      k <- k + 1
    }
  }
  return(abd)
}

Barro Colorado Island

We build a dataset according to BCI characteristics (Hubbell 2001): \(\theta=50\) and \(m=0.1\). The metacommunity follows a log-series distribution (Fisher, Corbet, and Williams 1943). We build it with a billion trees. This size is arbitrary: the important parameter is \(\theta\).

mc_size <- 1E9
# theta is also Fisher's alpha
alpha <- 50
# Number of species (Fisher, 1943)
(-alpha * log(alpha / (mc_size + alpha))) %>% 
  # integer needed
  as.integer() -> 
  n_species

The metacommunity contains 840 species. We simulate their abundances.

library("cli")
# Prepare a progress bar
mc_abundances <- rlseries(n = n_species, size = mc_size, alpha = alpha)

# Set the species names
names(mc_abundances) <- species_names[seq_along(mc_abundances)]

# Declare it is a vector of abundances (package entropart)
library("entropart")
mc_abundances %>% 
  as.AbdVector() ->
  mc_abundances

# Range
summary(mc_abundances)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##        1       44     4660  1246385   224242 61304615

# Rank-Abundance Curve
autoplot(mc_abundances, Distribution = "lseries")

Local community

We draw a local community from the metacommunity.

# Choose the size of the community
c_size <- 256
# Draw a random community according to the abundances of the metacommunity
rmultinom(1, size = c_size, prob = mc_abundances) %>% 
  # We need a vector, not a matrix
  as.vector() ->
  c_abundances
# Give the species a name
names(c_abundances) <- names(mc_abundances)
library("entropart")
# Clean-up the data
c_abundances %>% 
  # Eliminate species with 0 individual
  subset(c_abundances > 0) %>% 
  # Declare a vector of abundances
  as.AbdVector() ->
  c_abundances
# Show the RAC
autoplot(c_abundances)

Demographic drift

We set an arbitrary position to the trees of the local community just to visualize it.

library("dbmss")
# Build a dataframe for the marks of a point pattern
marks <- data.frame(
  # Species. Each line of the dataframe is an individual
  PointType = rep(names(c_abundances), times = c_abundances), 
  # Size of the points
  PointWeight = rep(1, c_size)
)
# Draw random locations in a square window to obtain a point pattern
c_size %>% 
  runifpoint() ->
  c_ppp
# Add the marks to the point pattern
marks(c_ppp) <- marks
# Make it a wmppp object (package dbmss) to make nice plots
c_wmppp <- as.wmppp(c_ppp)
# Set the factor levels of the point types identical to those of the metacommunity
# so that species of the metacommunity can be used in the local community
marks(c_wmppp)$PointType <- factor(
  marks(c_wmppp)$PointType, 
  levels = names(mc_abundances)
)
# Map the point pattern
autoplot(c_wmppp) +
  labs(colour = "Species", size = "Size")

The drift model is as follows: at each step, a tree dies and is replaced by a new tree from the species of the community, with respect to their probabilities.

#' One-step evolution of a community without migration
#'
#' @param community A wmppp object (package dbmss) where points are trees.
#' @param show If TRUE, a map of the evolved community is produced.
#'
#' @return A wmppp object with the modified community.
evolve_drift <- function(community, show = FALSE) {
  # Choose a tree
  focus_tree <- sample(seq_len(community$n), size = 1)
  # Species of the survivors
  survivors_sp <- marks(community)$PointType[-focus_tree]
  # Draw the species of the new tree
  marks(community)$PointType[focus_tree] <- 
    survivors_sp[sample(seq_along(survivors_sp), size = 1)]
  if (show) {
    # Copy the community
    community_to_plot <- community
    # Increase the size of the modified tree
    marks(community_to_plot)$PointWeight[focus_tree] <- 2
    # Plot the community (print is mandatory inside a function)
    print(autoplot(community_to_plot))
  }
  # return
  return(community)
}

# Test the function
c_wmppp %>% 
  evolve_drift(show = TRUE)

## Marked planar point pattern: 256 points
## Mark variables: PointWeight, PointType 
## window: rectangle = [0, 1] x [0, 1] units

Drift causes the local extinction of all species but one.

library("cli")
# Copy the community to modify it
the_community <- c_wmppp
# Number of species
richness <- length(unique(marks(the_community)$PointType))
# Prepare a long vector to save the richness along time
richness_time <- numeric(1E6)
# Count the steps
step <- 0
# Prepare a progress bar
if (interactive()) cli_progress_bar("Drifting")
while (richness > 1) {
  step <- step + 1
  # Replace a tree
  the_community <- evolve_drift(the_community)
  # Update the richness
  richness <- length(unique(marks(the_community)$PointType))
  # Save it
  richness_time[step] <- richness
  # Show the progress
   if (interactive()) cli_progress_update()
}
# Finalise the progress bar
if (interactive()) cli_progress_update(force = TRUE)
# Save the necessary time to reach a single species
drift_time <- step
# Plot the community
autoplot(the_community)

# Plot the decreasing richness
data.frame(Time = seq_len(step), Richness = richness_time[seq_len(step)]) %>% 
  ggplot() +
    geom_line(aes(x = Time, y = Richness))

Immigration

Immigration is the ability to replace the dead tree in the local community by a tree of a species of the metacommunity.

The parameter m is the probability for that to happen. The evolution function is rewritten.

#' One-step evolution of a community with migration
#'
#' @param community A wmppp object (package dbmss) where points are trees.
#' @param metacommunity A named vector of abundances or probability that describes the metacommunity.
#' @param m The migration parameter, that is the probability for a dead tree to be replaced by an immigrant species.
#'
#' @return A wmppp object with the modified community.
evolve_migrate <- function(community, metacommunity, m = 0) {
  # Choose a tree
  focus_tree <- sample(seq_len(community$n), size = 1)
  # Migration or local parent? Draw in a uniform distribution.
  u <- runif(1)
  if (u > m) {
    # Local parent: Species of the survivors
    survivors_sp <- marks(community)$PointType[-focus_tree]
    # Draw the species of the new tree
    marks(community)$PointType[focus_tree] <- 
      survivors_sp[sample(seq_along(survivors_sp), size = 1)]
  } else {
    # Migration : Draw the species of the new tree in the metacommunity
    marks(community)$PointType[focus_tree] <- 
      names(metacommunity)[which(rmultinom(1, size = 1, prob = metacommunity) == 1)]
  }
  # return
  return(community)
}

We make the community evolve the same amount of time as needed for complete drift and see what happens.

# Choose the migration probability
m <- 0.1

# Copy the community to modify it
the_community <- c_wmppp
# Number of species
richness <- length(unique(marks(the_community)$PointType))
# Prepare a long vector to save the richness along time
richness_time <- numeric(1E6)
# Prepare a progress bar
if (interactive()) cli_progress_bar(paste("Migration", m, ":"), total = drift_time)
for (step in seq_len(drift_time)) {
  # Replace a tree
  the_community <- evolve_migrate(the_community, mc_abundances, m)
  # Save the richness
  richness_time[step] <- length(unique(marks(the_community)$PointType))
  # Show the progress
  if (interactive()) cli_progress_update()
}
autoplot(the_community)

# Plot richness
data.frame(Time = seq_len(drift_time), Richness = richness_time[seq_len(drift_time)]) %>% 
  ggplot() +
    geom_line(aes(x = Time, y = Richness))

# Plot the rank-abundance curve
the_community %>% 
  # Consider the dataframe of marks
  marks() %>% 
  # Group by species
  group_by(PointType) %>% 
  # Count the numnber of individuals by species
  summarise(Abundance = n()) %>% 
  # Extract the column 
  pull(Abundance) %>% 
  # Make it an abundance vector (package entropart)
  as.AbdVector() %>% 
  # Plot the rank-abundance curve and try to fit a log-normal distribution
  autoplot(Distribution = "lnorm")

The average number of species can be calculated from the data, after eliminating the burn-in steps (i.e. before the pattern is stable).

mean(richness_time[10000:drift_time])

## [1] 42.87037

The equation for the expected number of species is in Etienne and Alonso (2007, eq. 4) but it requires knowing the size of the metacommunity.

The minimum richness is 32.

Exercises

Re-use the code above to run a few experiments about drift time and the shape of the local community.

The immigration parameter

From the BCI metacommunity, draw a local 1-ha community. Choose two immigration rates in turn: 0.1, as observed at BCI, and 0.01. Run 50,000 steps and observe the final RAC. Compare the number of species with the actual value.

BCI data is available in package vegan.

library("vegan")
data(BCI)
# Number of trees per ha
sum(BCI)/50

## [1] 429.14

# Number of trees in each 1-ha plot
rowSums(BCI)

##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
## 448 435 463 508 505 412 416 431 409 483 401 366 409 438 462 437 381 347 433 429 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
## 408 418 340 392 442 407 417 387 364 475 421 459 436 447 601 430 435 447 424 489 
##  41  42  43  44  45  46  47  48  49  50 
## 402 414 407 409 444 430 425 415 427 432

# Number of species in each 1-ha plot
rowSums(BCI > 0)

##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##  93  84  90  94 101  85  82  88  90  94  87  84  93  98  93  93  93  89 109 100 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##  99  91  99  95 105  91  99  85  86  97  77  88  86  92  83  92  88  82  84  80 
##  41  42  43  44  45  46  47  48  49  50 
## 102  87  86  81  81  86 102  91  91  93

# RAC of plot 1
BCI[1, ] %>%
  as.AbdVector() %>% 
  autoplot()

The effect of the size of the local community to drift

Draw small communities with 10 species, in log-normal distributions, ignoring the metacommunity. Make their size vary: 16, 32, 64, 128, 256 individuals. Run the drift model and note the necessary time for complete drift. Compare your results with figure 4.12 (page 106) of Hubbell (2001).

# Hint
rCommunity(1, size = 256) %>% 
  autoplot()

References

Etienne, Rampal S., and David Alonso. 2007. “Neutral Community Theory: How Stochasticity and Dispersal-Limitation Can Explain Species Coexistence.” Journal of Statistical Physics 128 (1-2): 485–510. https://doi.org/10.1007/s10955-006-9163-2.

Fisher, Ronald Aylmer, A. Steven Corbet, and Carrington Bonsor Williams. 1943. “The Relation Between the Number of Species and the Number of Individuals in a Random Sample of an Animal Population.” Journal of Animal Ecology 12: 42–58. https://doi.org/10.2307/1411.

Hubbell, Stephen P. 2001. The Unified Neutral Theory of Biodiversity and Biogeography. Princeton University Press.

ter Steege, Hans, Nigel C. A. Pitman, Daniel Sabatier, Christopher Baraloto, Rafael P. Salomão, Juan Ernesto Guevara, Oliver L. Phillips, et al. 2013. “Hyperdominance in the Amazonian Tree Flora.” Science 342 (6156): 1243092. https://doi.org/10.1126/science.1243092.

Simulation of the Neutral Model

Eric Marcon

2025-09-23