In previous labs we have focused on analyses within a population. Today we are going to begin comparing populations, including those of different species. In ecology, a community is a group two or more populations consisting of at least two different species that occupy the same geographical area at the same time. Therefore, community ecology is the study of these interactions.
In the first portion of lab today we will explore “Island Biogeography Theory” as a simple introduction. We will then continue our discussion on population growth but expand it to a community setting.
Island Biogeography Theory
Simple model of “Island Biogeography Theory”
Island Biogeography Theory is a fundamental concept in ecology, developed by Robert MacArthur and E. O. Wilson in 1967 (and reprinted in 2001). This theory describes the number of species on an island as a dynamic balance between immigration and extinction rates.
Theory
The “Island Biogeography Theory” posits that the number of species on an island reflects a balance between two processes:
Immigration: The arrival of new species on the island.
Extinction: The loss of species from the island.
These processes are influenced by two main factors:
Island Size: Larger islands support more species due to greater habitat diversity and resources. Additionally extinction rates will be lower on larger islands because of more resources.
Distance from Mainland: Islands closer to the mainland have higher immigration rates due to proximity.
A classic real-world example of this process was the re-colonization of birds to certain Indonesian islands during the decades following the eruption of the Krakatoa volcano in 1883.
Equations
Equilibrium Condition: The number of species on an island \(s\) is at equilibrium when the colonization rate \(C(s)\) is equal to the extinction rate \(E(s)\):
\[
C(s) = E(s) \tag{1}
\]
Colonization Rate: The colonization rate \(C(s)\) depends on the total number of possible species minus the number of species \((p - s)\) and the average colonization rate \(c\):
\[
C(s) = c(p - s) \tag{2}
\]
Extinction Rate: The extinction rate \(E(s)\) depends on the number of existing species \(s\) and the average extinction rate \(h\):
\[
E(s) = hs \tag{3}
\]
Equilibrium Number of Species: By substituting equations (2) and (3) into equation (1) and solving for \(s\), we can find the equilibrium number of species \(s^*\):
\[
s^* = \frac{cp}{c + h} \tag{4}
\]
R
This equilibrium reflects the balance between colonization and extinction processes on an island, influenced by the total species pool \(p\), colonization rate \(c\), and extinction rate \(h\).
Now let’s plot this relationship in R with potential number of species of 100, average colonization rate of 0.1, and average extinction rate of 0.05
# Define parametersp <-100c <-0.1h <-0.05# Define the colonization `C_s` and extinction rate `E_s` functionsC_s <-function(s) { c * (p - s) }E_s <-function(s) { h * s }# calculate the equilibrium value `s_star`s_star <- (c * p) / (c + h)# Create a sequence of species numbers for plottings_values <-seq(0, p, by=1)# Calculate colonization and extinction rates for each number of speciesC_values <-C_s(s_values)E_values <-E_s(s_values)# Create a data frame for plottingdf <-data.frame(s = s_values, C = C_values, E = E_values)head(df)
s C E
1 0 10.0 0.00
2 1 9.9 0.05
3 2 9.8 0.10
4 3 9.7 0.15
5 4 9.6 0.20
6 5 9.5 0.25
Now plot your results. Note that we can plot more than one line on a single plot by defining y within geom_line()
ggplot(data = df, aes(x = s)) +geom_line(aes(y = C, color ="Colonization Rate")) +geom_line(aes(y = E, color ="Extinction Rate")) +labs(title ="Island Biogeography: Colonization and Extinction Rates",x ="Number of Species (s)", y ="Rate", color ="Rate Type")
Where is the equilibrium point?
Where the colonization and extinction rates meet, around 66 species.
Extended model of “Island Biogeography Theory”
As mentioned above, MacArthur and Wilson extended their basic model by considering how colonization and extinction rates might be influenced by geographic properties of islands, specifically:
Distance from the mainland (\(d\)), which affects the colonization rate.
Island size (\(a\)), which affects the extinction rate.
Extended Model Equations
The colonization rate \(C(s)\), initially defined as \(c(p - s)\), is modified to account for distance from the mainland:
\[
C(s) = c(p - s)e^{-\phi d} \tag{5}
\]
where \(\phi\) is a fit parameter governing the distance decay of the colonization rate.
The extinction rate \(E(s)\), initially \(hs\), is modified to reflect the effect of island size:
\[
E(s) = se^{-\epsilon a} \tag{6}
\]
where \(\epsilon\) is a fit parameter governing the effect of area on extinction.
Solving for the equilibrium number of species \(s^*\) considering both modifications yields:
Incorporate these modifications into our model with \(\phi\) of 0.001, \(\epsilon\) of 0.001, distance from mainland of 100, and island area of 1,000.
# Set the parametersp <-100c <-0.1phi <-0.0001epsilon <-0.001d <-100a <-1000# Define the colonization and extinction rate functionsC_s <-function(s) { c * (p - s) *exp(-phi * d) }E_s <-function(s) { s *exp(-epsilon * a) }# Create a sequence of species numbers for plottings_values <-seq(0, p, by =1)# Calculate colonization and extinction rates for each number of speciesC_values <-sapply(s_values, C_s)E_values <-sapply(s_values, E_s)# Calculate the equilibrium number of species for the given distance and areas_star <- (c * p *exp(epsilon * a)) / (c *exp(epsilon * a) +exp(phi * d))# Create a data frame for plottingdf <-data.frame(s = s_values, C = C_values, E = E_values)
Plot your results
ggplot(data = df, aes(x = s)) +geom_line(aes(y = C, color ="Colonization Rate")) +geom_line(aes(y = E, color ="Extinction Rate")) +labs(title ="Island Biogeography: Colonization and Extinction Rates",x ="Number of Species (s)", y ="Rate", color ="Rate Type")
What would we need to do to decrease colonization rate?
Increase distance from mainland
What would we need to do to decrease extinction rate?
Increase island size
Real-World Data: Island Biogeography of Coral Reefs
# Select the presence-absence and `Tax. Unit 1` data columnsst2 <- st2 %>%select(1:16, `Tax. Unit 1`)# Perform the analysis for the community overallrates.st2 <-irregular_single_dataset(dataframe = st2, vector =3:16, c =0.001, e =0.001, assembly =TRUE, jacobian =TRUE)# Perform the analysis for each taxonomic grouprates.groups <-irregular_single_dataset(dataframe = st2, vector =3:16, c =0.001, e =0.001, column ="Tax. Unit 1", n =5, assembly =TRUE, jacobian =TRUE)
Plot the results
ggplot() +geom_text(data = rates.st2, aes(x = c, y = e, label ="Whole Community")) +geom_text(data = rates.groups, aes(x = c, y = e, label = Group)) +labs(title ="Colonization and Extinction Rates", x ="Colonization Rate", y ="Extinction Rate")
Whole Community: This point represents the average colonization and extinction rates for the entire community. Its position gives a baseline for comparing individual groups. If we assume our island community is stable, this point would be where the extinction and colonization rates meet.
Based on the whole commmunity as a baseline, what could we say about the other communities?
Acarina: This group has a very low colonization rate and a very low extinction rate, suggesting a stable group with few species arriving or becoming locally extinct. Its position near the origin could imply it’s a mature community or one with a small pool of potential colonizers.
Coleoptera: With a higher extinction rate than colonization rate, this group may be losing species over time, indicating it’s potentially less stable or experiencing some stressors that increase extinction rates.
Araneae and Psocoptera: These groups have relatively high colonization and extinction rates, suggesting dynamic communities with frequent species turnover. Such communities could be undergoing ecological succession or experiencing variable environmental conditions that allow many new species to colonize but also result in higher rates of species loss.
Lepidoptera and Hymenoptera: These groups have a high colonization rate and a low extinction rate. They are probably not a threatened group. Note that many members of this group are flying insects which may be why colonization rate is so high.
Metapopulations
In previous labs, we considered closed populations, with population size governed by birth and death. Now we are going to consider all important terms of the BIDE model (Pulliam 1988): birth, immigration, death, and emigration. Therefore, we are going to expand our models to include the interaction among populations.
First described by H. G. Andrewartha and Charles Birch in 1954 (Andrewartha, Birch, et al. 1954) and later being named by Richard Levins in 1969 (Levins 1969) in the context of pest management, metapopulations is the term for these “populations of populations.”
Currently, in the context of conservation biology metapopulation theory has largely replaced the theory of island biogeography. In particular, the Levin’s Metapopulation Model helps us understand how habitat loss may affect species persistence.
Levin’s Metapopulation Model
Recall the simplest Lotka-Volterra model we covered:
Instead of carrying capacity limiting population growth as seen in the Lotka-Volterra model, in the Levin’s Metapopulation model the self-limiting portion \((1 - P)\) is the proportion of empty patches. \(P + (1 - P)\) must always equal what?
\[
P + (1 - P) = 1
\]
Because of the inclusion of the term \((1 - P)\), this equation represents a closed metapopulation:
Recall from past labs this semester, the general layout for partial derivative functions with the deSolve()(Soetaert, Petzoldt, and Setzer 2010) package:
function_name <- function(y, x, parms) {
X <- x[1]
with(as.list(parms), {
dX.dy <- X # complete formula here
return(list(c(dX.dy)))
})
}
function_name.out <- deSolve::ode(x = X0, times = times, func = function_name, parms = parms)
Write a function levins() for the Levin’s metapopulation model with parms <- c(c = 0.15, m = 0.06), P0 <- 0.01, and times <- 1:100
library(deSolve)parms <-c(c =0.15, m =0.06)P0 <-0.01times <-1:100levins <-function(t, y, parms) { P <- y[1]with(as.list(parms), { dP.dt <- c * P * (1- P) - (m * P)return(list(c(dP.dt))) })}levins_out <- deSolve::ode(y = P0, times = times, func = levins, parms = parms)
Now plot (remember to use ggplot2, we have to first convert our dataset to a dataframe)
levins_out_df <- levins_out %>%as.data.frame() %>%rename(time =1, V1 =2)ggplot(data = levins_out_df, aes(x = time, y = V1)) +geom_line() +labs(x ="Time", y ="Occupied Patches", title ="Occupied Patches over Time")
Habitat Destruction
Russel Lande (Lande 1987) found that habitat destruction \(D\) can be easily integrated into the simpler form of the Levin’s metapopulation model where \(D = 1\) is complete destruction and \(D = 0\) is the default Levin’s metapopulation model.
Integrate habitat destruction D = 0.2 into the Levin’s metapopulation model to create a function for the lande() model over period times <- 1:200.
parms <-c(c =0.15, m =0.06, D =0.2)P0 <-0.01times <-1:200lande <-function(t, P, parms) { P <- P[1]with(as.list(parms), { dP.dt <- ( c * P * (1- D - P) ) - ( m * P )return(list(c(dP.dt))) })}lande_out <- deSolve::ode(y = P0, times = times, func = lande, parms = parms)
Plot the Lande model output using ggplot2
lande_out_df <- lande_out %>%as.data.frame() %>%rename(time =1, V1 =2)ggplot(data = lande_out_df, aes(x = time, y = V1)) +geom_line() +labs(x ="Time", y ="Population Size", title ="Population Dynamics over Time")
What happens to the curve as you vary D?
lande_plots <-function(D) { parms <-c(c =0.15, m =0.06, D = D) lande_out <- deSolve::ode(y = P0, times = times, func = lande, parms = parms)# Convert output to a dataframe and rename columns df <-as.data.frame(lande_out) %>%rename(time =1, V1 =2) %>%mutate(D =as.factor(D))return(df)}D_values <-c(0, 0.2, 0.3, 0.5, 0.6)
Plot the results:
lande_df <-do.call(rbind, lapply(D_values, lande_plots))ggplot(lande_df, aes(x = time, y = V1, color = D, group = D)) +geom_line() +labs(x ="Time", y ="Population Size", title ="Population Dynamics for Different Disturbance Levels (D)")
What do you observe?
As D increases, the carrying capacity drops.
However, there is a threshold point, where the population goes extinct above a certain point. remember from lecture this is the “extinction threshold” where the equilibrium fraction of occupied patches is zero:
\[ D = 1 - \frac{m}{c}\]
This is an important reminder of the importance of empty patches.
Metapopulation Dynamics
In this lab, we have already worked with simple models of metapopulations and dispersal using the Levin’s Metapopulation Model/Lande Model. Recall from lecture the four perspectives to characterize how dispersal rates and environmental heterogeneity interact to influence species diversity in metacommunities.
Patch dynamics perspective
Species sorting perspective
Mass effects perspective
Neutral perspective
These perspectives cannot be fully represented in a Levin’s Metapopulation Model/Lande Model alone.
note that the exact code in the article was written in Julia, not R
This package allows you to easily model the following simultaneously and integrates stochaisticity:
Density independent responses to abiotic conditions
Density dependent competition
Dispersal
As the most basic implementation, you can vary the amount of patches and species with the following arguments:
simulate_MC(patches =5, species =4)
Warning in simulate_MC(patches = 5, species = 4): Environment is not spatially
autocorrelated in the current version of the package and so results will differ
from Thompson et al. 2020 Ecology Letters.
Review the “help page” for the function simulate_MC() using your console (?simulate_MC). With this information, simulate each of the four perspectives using the simulate_MC() function and discuss briefly what you observe.
# ?simulate_MC()
Patch dynamics perspective
High or lower dispersal, low response to environmental heterogeneity (note that a high number for env_niche_breadth means low response, because more of a generalist)
simulate_MC(patches =5, species =4, dispersal =0.9, env_niche_breadth =10)
Warning in simulate_MC(patches = 5, species = 4, dispersal = 0.9,
env_niche_breadth = 10): Environment is not spatially autocorrelated in the
current version of the package and so results will differ from Thompson et al.
2020 Ecology Letters.
landscape: Scatter plot that shows where the plots are located
disp_mat: Heatmap that shows dispersal rates and how they differ between close and distant plots
For the “patch dynamics perspective,” patches are homogeneous so we can either see high or low dispersal because all patches identical. In this case we see high dispersal (lowest heatmap values along the diagonal where to.patch and from.patch values equal).
dynamics.df: This timeseries plot show how populations grow and decline due to interactions with other species, environmental conditions, and possibly stochastic events.
Time series curves are asynchronous with one another. Overall, how patch dynamics work is that there is a trade-off between competition and colonization which we can see here.
env.df: This response curve shows how the environment changes in each patch over time
For the “patch dynamics perspective” there is no difference in patches, so the response to environmental differences is the same over time
int_mat: A competition matrix that visualizes the strength of competitive interactions between pairs of species
The default setting for this package is that 30% of species will be dominant, so we are seeing that affect here. Remember there is a trade-off between competition and colonization.
These time series plots represent the environmental optima over time for each species faceted by each patch
Because each species and each patch are equal in the “patch dynamics perspective,” each patch is taken over by a subset of species and the rest are excluded. That is why we see these “ranked” optimums of 3 > 2 > 1 > 5 > 4. This reflects the order of colonization (and exclusion) but not patch quality since all patches are equal.
Species sorting perspective
Low dispersal, high response to environmental heterogeneity
simulate_MC(patches =5, species =4, dispersal =0.1, env_niche_breadth =0.5)
Warning in simulate_MC(patches = 5, species = 4, dispersal = 0.1,
env_niche_breadth = 0.5): Environment is not spatially autocorrelated in the
current version of the package and so results will differ from Thompson et al.
2020 Ecology Letters.
High dispersal, high response to environmental heterogeneity
simulate_MC(patches =5, species =4, dispersal =0.8, env_niche_breadth =0.5)
Warning in simulate_MC(patches = 5, species = 4, dispersal = 0.8,
env_niche_breadth = 0.5): Environment is not spatially autocorrelated in the
current version of the package and so results will differ from Thompson et al.
2020 Ecology Letters.
Mid-dispersal, low response to environmental heterogeneity
simulate_MC(patches =5, species =4, dispersal =0.5, env_niche_breadth =10)
Warning in simulate_MC(patches = 5, species = 4, dispersal = 0.5,
env_niche_breadth = 10): Environment is not spatially autocorrelated in the
current version of the package and so results will differ from Thompson et al.
2020 Ecology Letters.
Pulliam, H Ronald. 1988. “Sources, Sinks, and Population Regulation.”The American Naturalist 132 (5): 652–61.
Soetaert, Karline, Thomas Petzoldt, and R. Woodrow Setzer. 2010. “Solving Differential Equations in r: Package desolve” 33. https://doi.org/10.18637/jss.v033.i09.
Stevens, Martin Henry Hoffman. 2009. A Primer of Ecology with r. Springer.
