September 25, 2017, Monday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links


MBW:Habitat heterogentiy-area tradeoffs and species diversity

From MathBio

Jump to: navigation, search

Overview

  • Mathematics used: agent-based model, individual-based model, discrete time Markov process
  • Biological system studied: local and regional community dynamics

Executive Summary

There is an inherent tradeoff between habitat diversity and the amount of habitat of a particular type that can exist in a finite spatial area. Furthermore, habitat diversity and the amount of usable habitat for a species (effective area) are expected to have different effects on species richness. Allouche et al. (2012) provide a rigorous empirical and theoretical treatment of the area-heterogeneity tradeoff with respect to species diversity. They find that the relationship between habitat heterogeneity and species richness (the number of species in a locality) can be unimodal. Here we reiterate their results and extend the concept of habitat homogeneity to account for different kinds of nestedness.

Background

Figure 1. Illustration of the habitat heterogeneity-area tradeoff. Each of the nine boxes represent different localities, each consisting of 100 cells. Habitat types are represented by different colors. As more habitat types are introduced, the number of area that can be allocated to each decreases.

Anthropogenic habitat loss and homogenization is happening at accelerating rates, contributing to biodiversity loss and major ecosystem changes (Hooper et al. 2012). Understanding the mechanistic relationship between habitat area, habitat homogeneity, and biodiversity is crucial for interpreting and predicting current and future changes in ecosystems and the services that they provide. One framework for building this understanding, presented by Allouche et al. (2012), begins with an understanding that there is an inherent trade-off between habitat area and homogeneity. This is important from a community ecology perspective, because habitat heterogeneity and effective area available for a species have different effects.

Habitat heterogeneity facilitates colonization of a locality: if a diversity of resources are available, then it is likely that more species will be able to find a usable resource in the locality. This idea is often framed in terms of local and regional processes. Specifically, local areas are colonized by individuals that occur in a regional pool of potential colonists. For example, if there is a regional pool of species that each has a color preference for habitat type in the figure to the right, then as more habitat types are available, more species can colonize the locality from the regional pool.

Area affects carrying capacity of populations. Populations occuring in large areas with high carrying capacities experience lower risk of local extinction due to demographic stochasticity (Lande 1993). Thus, the trade-off between habitat heterogeneity and effective area is expected to have consequences for the colonization of new species in a locality and the persistence (or local extinction) of species in the locality once established.

Model

To theoretically explore the consequences of the habitat heterogeneity-area tradeoff, Allouche et al. (2012) built an individual based model that formalizes the conceptual model explained above. They used a continuous time Markov processes, which is discretized in this example.

This model is spatially implicit, with A equally connected sites. Each site falls on an environmental condition axis, receiving some value E that characterizes local conditions. The environmental conditions for each site are uniformly distributed between two values that dictate the range of environmental conditions in a focal area. The local range of environmental conditions is a subset of some global range. There are N species in the regional pool that can colonize habitat patches. Each species has some environmental optimum \mu _{i}, and some niche width \sigma _{i}, which together define a Gaussian function for the probability of establishment of species i given a colonization attempt in patch j and a habitat patch environmental condition E_{j}.

Pr(establish)=Z_{i}e^{{-2\sigma _{i}^{2}(E_{j}-\mu _{i})^{2}}}

Z_{i}=(\int _{{E_{{min}}^{R}}}^{{E_{{max}}^{R}}}e^{{-2\sigma _{i}^{2}(E-\mu _{i})^{2}}}dE)^{{-1}}

where Z_{i} is a normalization constant that makes the overall probability of establishing within the global range of environmental conditions equal for all species in the regional pool. E_{{min}}^{R} and E_{{max}}^{R} are the global minima and maxima for environmental conditions. Species idea conditions, \mu _{i} are uniformly distributed across the global range of environmental conditions. Two methods are used to assign niche width values: uniform niche width is created by fixing \sigma _{i} for all species equal to one value, and habitat specialism and generalism occurs when niche widths are uniformly distributed between two values \sigma _{{min}} and \sigma _{{max}}.

Niche space when species have identical niche widths
Niche space when species have different niche widths, so that some are habitat specialists, and others generalists

The model as presented here uses discrete time steps. At each time step, individuals die with some probability, equal across species, and if they produce some number of offspring that is a Poisson random variable, with some expected value that is the same across all species. Offspring attempt to colonize sites at random. Only unoccupied sites can be successfully colonized. If an attempt is made, a Bernoulli trial with probability from the above equation is made to determine whether the individual can successfully establish. One individual from every species in the regional pool has some probability of attempting to colonize each patch at each time step. If multiple individuals attempt to colonize the same patch simultaneously, a successful individual is randomly selected among those attempting.

Richness time series iterated 1 and 36 times for a set of parameters. ERmin=-50, ERmax=50, Emin=-50, Emax=50, pM=0.9, sigmin=1, sigmax=30, R=1, N=100, A=100, I=.9, timesteps=300.

For any set of parameters the model gives a time series of the number of unique species in the locality - species richness (Figure 3).

To investigate the role of habitat heterogeneity on species richness we construct replicate localities that vary in thier degree of habitat homogeneity but not in total area. For each locality the model is run with identical parameters. As we would expect from the conceptual model, a unimodal relationship between richness and habitat heterogeneity emerges, with increasing niche space facilitating increases in richness at low heterogeneities, and stochastic population extirpations reducing richness in highly heterogeneous localities.

Mean richness for 200 localities that vary in heterogeneity, but not total area, with the same regional pool of species for each iteration. ERmin=-50, ERmax=50, Emin=varied, Emax=varied, pM=0.9, sigmin=1, sigmax=30, R=1, N=100, A=100, I=.9, timesteps=300.


Extensions

Habitat heterogeneity is not necessarily sequentially nested. For instance, some localities may vary in heterogeneity, but fall along different ranges of a resource or condition axis. This is a pattern that could emerge when considering islands, for example, that have habitat conditions that emerge from priority colonization effects. Some islands could be dominated by different conditions early on, leading to non-nested habitat conditions among them. Extending the concept of homogeneity from the original model, we can construct intervals that vary both in heterogeneity and condition by randomizing the position of the previous intervals across the global range of conditions. Repeating the process outlined above, we find that this non-sequentially nested homogeneity leads to cases where species richness is far lower than the expectation for sequentially nested localities.

Mean richness for 200 localities that vary in heterogeneity and conditions, but not total area, with the same regional pool of species for each iteration. ERmin=-50, ERmax=50, Emin=varied, Emax=varied, pM=0.9, sigmin=1, sigmax=30, R=1, N=100, A=100, I=.9, timesteps=300.

These cases with reduced richness coincide with reduced overlap between Emin and Emax (the range of local environmental conditions) and the niche space of species in the regional pool. If an interval falls along a part of the resource axis where few species have high colonization probabilities, then the pool of potential colonists is reduced.

Discussion

The model of Allouche et al. (2012) provides a unified framework for disentangling the effects of habitat loss and homogenization in the context of regional source pools of colonizing species. In the original article, an empirical evaluation of the habitat area-heterogeneity tradeoff in montane birds demonstrates support for a unimodal relationship between habitat heterogeneity and species diversity, conditional on niche width (elevational occurrence range). The model has applied value in that it highlights the importance of habitat area and the degree of heterogeneity as important for conservation. Many conservationists would tend to agree that larger areas are better for preserving biodiversity, and more diverse habitats do the same. However, it is important to consider how the degree of heterogeneity impacts population sizes due to the habitat area heterogeneity trade-off.

R code

## Define model function ##
bortzIBM <- function(A=100, N=100, ERmin=-30, ERmax=30, Emin=-30, Emax=30, 
                    sigmin=1, sigmax=1, timesteps=500, pM=.1, pR=1, R=2, I=0.1){
  E <- runif(A, Emin, Emax) # habitat patch environment type
  mu.i <- runif(N, ERmin, ERmax) # optimum environment
  sigma <- runif(N, sigmin, sigmax)  # niche width
  
  Z <- array(dim=c(N)) # normalization constant
  
  # now calculating all elements of Z coefficient 
  # (ensures all species have equal Pr(establishment) in regional pool)
  for (i in 1:N){
    integrand <- function(E) {
      exp(-((E - mu.i[i]) ^ 2) / (2 * sigma[i] ^ 2))
    }
    res <- integrate(integrand, lower=ERmin, upper=ERmax)
    Z[i] <- 1 / res$value
  }
  
  # probability of establishment
  Pcol <- array(dim=c(A, N))
  for (i in 1:A){
    for (j in 1:N){
      Pcol[i, j] <- Z[j] * exp(-((E[i] - mu.i[j]) ^ 2) / (2*sigma[j] ^ 2))
    }
  }
  
  # store niche data
  species <- rep(1:N, each=A)
  Sigma <- ifelse(length(sigma==1), rep(sigma, A*N), rep(sigma, each=A))
  E <- rep(E, N)
  Pr.estab <- c(Pcol)
  niche.d <- data.frame(species, E, Pr.estab, Sigma)
  niche.d <- niche.d[with(niche.d, order(species, E)),]
  
  # initialize output objects
  state <- array(0, dim=c(timesteps, A, N))
  richness <- rep(NA, timesteps)
  richness[1] <- 0
  pr.occ <- rep(NA, timesteps)
  pr.occ[1] <- 0
  
  pb <- txtProgressBar(min = 0, max = timesteps, style = 3)
  
  for (t in 2:timesteps){
    state[t,,] <- state[t-1,,]
    
    ## DEATHS ##
    deaths <- array(rbinom(A*N, 1, c(state[t,,])*pM), dim=c(A, N))
    state[t,,] <- state[t,,] - deaths 
    
    ## BIRTHS ##
    pot.fecundity <- array(rpois(A * N, lambda = c(state[t,,] * R)), dim=c(A, N)) # potential number of offspring
    repro <- array(rbinom(A*N, 1, pR), dim=c(A, N)) # whether reproduction actually occurs
    fecundity <- repro * pot.fecundity
    sum.fec <- apply(fecundity, 2, sum) # number of offspring per species
    
    ## OFFSPRING COLONIZE EMPTY SITES ##
    occupancy <- apply(state[t,,], 1, max)
    if (sum(occupancy) < A & sum(sum.fec) > 0){ # if empty sites & new offspring
      empty.sites <- which(occupancy == 0)
      occ.sites <- which(occupancy == 1)
      # randomly assign sites (empty & filled) to each offspring individual
      col.sites <- sample(1:A, sum(sum.fec), replace=T)
      col.spec <- rep(1:N, times=sum.fec) # how many of each species colonizing
      colonizing.offspring <- array(0, dim=c(A, N)) # how many of each species colonizing each site
      for(i in 1:length(col.sites)){
        colonizing.offspring[col.sites[i], col.spec[i]] <- colonizing.offspring[col.sites[i], col.spec[i]] + 1 
      }
      # offspring attempting to colonize occupied sites fail to displace
      colonizing.offspring[occ.sites,] <- 0
      
      # which colonizing offspring can actually establish?
      binom.mat <- ifelse(colonizing.offspring > 0, 1, 0)
      colonists <- array(rbinom(n = A * N, 
                                size = c(colonizing.offspring),
                                prob = c(binom.mat * Pcol)),
                         dim=c(A, N))
      
      # are there colonization conflicts (> 1 individual trying to colonize each site?)
      attempting <- apply(colonists, 1, sum) # number indiv. attempting to colonize each site
      if (any(attempting > 1)){
        # resolve colonization conflicts
        conflicts <- which(attempting > 1) # which sites have conflicts
        for (k in conflicts){ # for each conflict
          # how many of individuals of each species are attempting to simultaneously colonize?
          n.attempting <- rep(1:N, times = colonists[k,])
          # randomly select one successful from those attempting
          successful <- sample(n.attempting, size=1)
          new.row <- rep(0, length.out=N)
          new.row[successful] <- 1
          colonists[k,] <- new.row # individual becomes the only colonist
        }
      }
      # add successful colonists
      state[t,,] <- state[t,,] + colonists
    } 
    # end of offspring colonization
    
    ## IMMIGRANTS COLONIZE EMPTY SITES ##
    occupancy <- apply(state[t,,], 1, max)
    if(sum(occupancy) < A){
      empty.sites <- which(occupancy == 0)
      # which species immigrate to each site?
      immigration <- array(rbinom(length(empty.sites)*N,
                                  1, I), dim=c(length(empty.sites), N))
      # which immigrants establish?
      Pest <- immigration * Pcol[empty.sites, ]
      establishment <- array(rbinom(length(Pest), 1, c(Pest)),
                             dim=c(length(empty.sites), N))
      
      # resolve conflicts arising from simultaneous colonization
      col.attempts <- apply(establishment, 1, sum)
      if (any(col.attempts > 1)){ # if individuals are trying to simultaneously colonize
        conflicts <- which(col.attempts > 1) # which empty sites have conflicts
        for (k in conflicts){ # for each conflict
          # how many of individuals of each species are attempting to simultaneously colonize?
          attempting <- rep(1:N, times = establishment[k,])
          # successful individual randomly selected from those attempting
          successful <- sample(attempting, size=1)
          new.row <- rep(0, length.out=N)
          new.row[successful] <- 1
          establishment[k,] <- new.row
        }
      }
      # add establishing immigrants
      state[t, empty.sites, ] <- state[t, empty.sites,] + establishment
    }

    
    richness[t] <- sum(apply(state[t,,], 2, max))
    pr.occ[t] <- length(which(state[t,,] == 1)) / A
    setTxtProgressBar(pb, t)
  }
  # check to confirm there are not multiple indivduals in any sites
  stopifnot(all(apply(state[t,,], 1, sum) < 2))
  return(list(richness=richness, pr.occ = pr.occ, 
              state=state,
              niche.d=niche.d))
}

## Analysis ##
#----------------------------------------------------------------------#
# Investigate the effect of resource heterogeneity on species richness #
#----------------------------------------------------------------------#
# 1. generate random intervals that are subsets of the global interval
# Uniform interval boundaries
ERmin <- -50
ERmax <- 50
global.median <- median(c(ERmin, ERmax))
n.intervals <- 200
lower.limits <- seq(ERmin, global.median-.5, length.out=n.intervals)
upper.limits <- seq(ERmax, global.median, length.out=n.intervals)
hab.het <- upper.limits - lower.limits # habitat heterogeneity width
hist(hab.het)

# Run sequentially nested localities
require(doMC)
n.iter <- 4 # number of iterations per core
ts <- 300
registerDoMC(8) # register however many cores you want to use
output <- foreach(i=1:getDoParWorkers(), .combine="c") %dopar% {
  rich <- array(NA, dim=c(n.intervals, n.iter, ts))
  hrvec <- NULL
  timesteps <- NULL
  interv <- NULL
  iteration <- NULL
  core <- NULL
  for (int in 1:n.intervals){
    for (iter in 1:n.iter){
      Emin <- lower.limits[int]
      Emax <- upper.limits[int]
      result <- bortzIBM(pM=.9, Emin=Emin, Emax=Emax, sigmin=1, sigmax=30, 
                        pR=1, R=1, N=100, A=100, I=.9, timesteps=ts)
      rich[int, iter, ] <- result$richness
      hrvec <- c(hrvec, result$richness)
      timesteps <- c(timesteps, 1:ts)
      interv <- c(interv, rep(int, ts))
      iteration <- c(iteration, rep(iter, ts))
      core <- c(core, rep(i, ts))
    }
  }
  return(data.frame(hrvec, timesteps, interv, iteration, core))
}

require(ggplot2)
ggplot(result$niche.d, aes(x=E, y=Pr.estab, color=factor(species))) + 
  geom_line(size=1.5) +
  theme_classic() + xlab("Environmental conditions")+ 
  ylab("Probability of establishing") + theme(legend.position="none") +
  ggtitle("Ecological niches of species in regional pool")

# gather data in dataframe
rich <- NULL
times <- NULL
intvl <- NULL
iteration <- NULL
core <- NULL
for(i in 1:getDoParWorkers()){
  rich <- c(rich, output[[which(names(output) == "hrvec")[i]]])
  times <- c(times, output[[which(names(output) == "timesteps")[i]]])
  intvl <- c(intvl, output[[which(names(output) == "interv")[i]]])
  iteration <- c(iteration, output[[which(names(output) == "iteration")[i]]])
  core <- c(core, output[[which(names(output) == "core")[i]]])
}

DATA <- data.frame(rich, times, intvl, iteration, core)

# for each interval * iteration * core combination, calculate a mean & sd of richness
runs <- n.intervals*n.iter*getDoParWorkers()
hrich <- NULL
srich <- NULL
interval <- NULL
itertn <- NULL
for (int in 1:n.intervals){
  for (iter in 1:n.iter){
    for (core in 1:getDoParWorkers()){
      rows <- which(DATA$iteration == iter & DATA$intvl == int & DATA$core==core)
      mhrich <- mean(DATA$rich[rows[-(1:100)]])
      hrich <- c(hrich, mhrich)
      sdrich <- sd(DATA$rich[rows[-(1:100)]])
      srich <- c(srich, sdrich)
      interval <- c(interval, int)
      itertn <- c(itertn, iter)
    }
  }
}


## Now what happens when the intervals are not centered around one value? #
dev <- ERmax - upper.limits # deviation from center
adj <- runif(n.intervals, 0, dev)
signs <- sample(c(-1, 1), size=n.intervals, replace=T)
adj <- adj*signs
intervals2 <- cbind(lower.limits, upper.limits)
intervals2 <- intervals2  + adj
lower.limits2 <- intervals2[,1]
upper.limits2 <- intervals2[,2]
cbind(lower.limits2, upper.limits2)
hab.het2 <- upper.limits2 - lower.limits2 # habitat heterogeneity width

require(doMC)
n.iter <- 4 # number of iterations per core
ts <- 300
output2 <- foreach(i=1:getDoParWorkers(), .combine="c") %dopar% {
  rich <- array(NA, dim=c(n.intervals, n.iter, ts))
  hrvec <- NULL
  timesteps <- NULL
  interv <- NULL
  iteration <- NULL
  core <- NULL
  for (int in 1:n.intervals){
    for (iter in 1:n.iter){
      Emin <- lower.limits2[int]
      Emax <- upper.limits2[int]
      result <- bortzIBM(pM=.9, Emin=Emin, Emax=Emax, sigmin=100, sigmax=100, 
                         pR=1, R=1, N=100, A=100, I=.9, timesteps=ts)
      rich[int, iter, ] <- result$richness
      hrvec <- c(hrvec, result$richness)
      timesteps <- c(timesteps, 1:ts)
      interv <- c(interv, rep(int, ts))
      iteration <- c(iteration, rep(iter, ts))
      core <- c(core, rep(i, ts))
    }
  }
  return(data.frame(hrvec, timesteps, interv, iteration, core))
}

# gather data in dataframe
rich2 <- NULL
times2 <- NULL
intvl2 <- NULL
iteration2 <- NULL
core2 <- NULL
for(i in 1:getDoParWorkers()){
  rich2 <- c(rich2, output2[[which(names(output2) == "hrvec")[i]]])
  times2 <- c(times2, output2[[which(names(output2) == "timesteps")[i]]])
  intvl2 <- c(intvl2, output2[[which(names(output2) == "interv")[i]]])
  iteration2 <- c(iteration2, output2[[which(names(output2) == "iteration")[i]]])
  core2 <- c(core2, output2[[which(names(output2) == "core")[i]]])
}

DATA2 <- data.frame(rich2, times2, intvl2, iteration2, core2)

# for each interval * iteration * core combination, calculate a mean & sd of richness
runs <- n.intervals*n.iter*getDoParWorkers()
hrich2 <- NULL
srich2 <- NULL
interval2 <- NULL
itertn2 <- NULL
for (int in 1:n.intervals){
  for (iter in 1:n.iter){
    for (core in 1:getDoParWorkers()){
      rows <- which(DATA2$iteration2 == iter & DATA2$intvl2 == int & DATA2$core2==core)
      mhrich <- mean(DATA2$rich2[rows[-(1:100)]])
      hrich2 <- c(hrich2, mhrich)
      sdrich <- sd(DATA2$rich2[rows[-(1:100)]])
      srich2 <- c(srich2, sdrich)
      interval2 <- c(interval2, int)
      itertn2 <- c(itertn2, iter)
    }
  }
}

# Comparing the nested vs. dispersed results
par(mfrow=c(2,2))
plot(rep(hab.het, each=32), hrich, xlab="Habitat heterogeneity", ylab="Mean richness", col="white",
     col.axis="white", col.lab="white", col.main="white", col.sub="white", fg="white")
plot(rep(hab.het2, each=32), hrich2, xlab="Habitat heterogeneity", ylab="Mean richness", col="white",
     col.axis="white", col.lab="white", col.main="white", col.sub="white", fg="white")
plot(x=NULL, y=NULL, ylim=c(ERmin*1.05, ERmax*1.05), xlim=c(0, n.intervals), type="n",
     ylab="Environmental conditions", xlab="Interval",
     col.axis="white", col.lab="white", col.main="white", col.sub="white", fg="white")
for (int in n.intervals:1){
  lines(y=c(lower.limits[int], upper.limits[int]), x=c(201-int, 201-int), col="white")
}
plot(x=NULL, y=NULL, ylim=c(ERmin*1.05, ERmax*1.05), xlim=c(0, n.intervals), type="n",
     ylab="Environmental conditions", xlab="Interval",
     col.axis="white", col.lab="white", col.main="white", col.sub="white", fg="white")
for (int in n.intervals:1){
  lines(y=c(lower.limits2[int], upper.limits2[int]), x=c(201-int, 201-int), col="white")
}


References

Allouche O, Kalyuzhny M, Moreno-Rueda G, Pizarro M, Kadmon R. Early Edition. Area-heterogeneity tradeoff and the diversity of ecological communities. Proceedings of the National Academy of Sciences.