
MBW:Habitat heterogentiyarea tradeoffs and species diversityFrom MathBioContentsOverview
Executive SummaryThere 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 areaheterogeneity 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. BackgroundAnthropogenic 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 tradeoff 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 tradeoff 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. ModelTo theoretically explore the consequences of the habitat heterogeneityarea 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 equally connected sites. Each site falls on an environmental condition axis, receiving some value 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 species in the regional pool that can colonize habitat patches. Each species has some environmental optimum , and some niche width , 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 .
where 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. and are the global minima and maxima for environmental conditions. Species idea conditions, 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 for all species equal to one value, and habitat specialism and generalism occurs when niche widths are uniformly distributed between two values and . 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. 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.
ExtensionsHabitat 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 nonnested 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 nonsequentially nested homogeneity leads to cases where species richness is far lower than the expectation for sequentially nested localities. 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. DiscussionThe 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 areaheterogeneity 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 tradeoff. 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[t1,,] ## 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(201int, 201int), 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(201int, 201int), col="white") } ReferencesAllouche O, Kalyuzhny M, MorenoRueda G, Pizarro M, Kadmon R. Early Edition. Areaheterogeneity tradeoff and the diversity of ecological communities. Proceedings of the National Academy of Sciences. 