12 Consumer-resource Competition and Mutualisms

In Chapter 7, we used Lotka-Volterra models of competition and mutualism to learn about some fundamental properties of (-,-) and and (+,+) interactions. Here we add some model complexity to illustrate a more mechanistic approach to modeling greater biological reality.

In Chapter 8, on consumer-resource interactions, we learned a variety of approaches to representing different types of resources or prey, and different ways consumers consume them. Please read chapter 8 before proceeding here. There are a lot of terms and parameters we use here that are more fully explained in Chapter 8.

12.1 Competition for one resource

Let’s imagine we have competition between two species for the shared limiting inorganic resource, \[\frac{dN_1}{dt}= N_1\left(e_1a_1R - d_1\right)\] \[\frac{dN_2}{dt}= N_2\left(e_2a_2R - d_2\right)\] where $e$ is assimilation efficiency, $a$ is consumption rate, and $d$ is mortality. As we can see, these two species do not interact directly–$N_2$ is not in the expression for $dN_1/dt$ and vice versa.⁷⁶

The expression for $dR/dt$ includes a renewing resource, \[\frac{dR}{dt} = r(S-R) - e_1a_1N_1R - e_2a_2N_2R\] where $S$ is background resource supply rate, and $r$ is resource renewal rate.

Here we graph mass-specific or per capita growth rates as functions of resource concentrations. These curves are simply straight lines (Fig. @ref{fig:dN}) described by, \[y = e_ia_i x - d_i\] where $y$ is the per capita or mass-specific growth rate (dN/Ndt), and $x$ is resource concentration or density ($R$).

e1 <- 0.1; a1 <- 1; d1 <- 0.3
e2 <- 0.1; a2 <- .5; d2 <- 0.1

{
  # setting the margins of the plot
  par(mar=(c(5, 4,.5,.5)))
  
  # making the graph
  curve(e1*a1*x - d1, 0,10, ylab="dN/(Ndt)", xlab="Resource concentration")
  curve(e2*a2*x - d2, lty=2, col=2, add=TRUE)
  abline(h=0, lty=3)
  legend("topleft", c("Sp. 1", "Sp. 2"), lty=1:2, col=1:2, bty='n')
  
  #adding the arrows
  x1 <- d1/(e1*a1); x2 <- d2/(e2*a2); len <- 0.1
  arrows(c(x1,x2), c(-len, len), c(x1,x2), c(0,0), length=0.1)
  }

Mass-specific growth rates of two species across a range of resource concentrations. Note that where growth rates equal death rates, consumer growth is zero (arrows) and this determines each species' $R^*$s.

Figure 12.1: Mass-specific growth rates of two species across a range of resource concentrations. Note that where growth rates equal death rates, consumer growth is zero (arrows) and this determines each species’ $R^*$s.

We can calculate the equilibrium resource concentration that each species would create on its own, their $R^*$s. Recall that we find this by setting each species per capita growth rates to zero, and solving for $R$.

R.star.sp1 <- d1/(e1*a1)
R.star.sp2 <- d2/(e2*a2)
c(R.star.sp1, R.star.sp2)

## [1] 3 2

From Fig. 12.1, we see that at high resource concentrations, species 1 grows faster than species 2. However, at low resource concentrations, species 2 grows faster than species 1. Resource concentrations change as consumers take up resources, so when does the whole system reach an equilibrium?

Both species (Fig. 12.1) increase in biomass as long as resources are high enough to allow growth, whereas resources will decline as consumers grow more and more and take up increasing amounts of resources. Eventually consumer biomass is high enough to drive down resources so low that mass-specific growth rate is zero. This will happen at different resource concentrations for our two species. We see that species 2 has a positive growth rate when species 1 does not (Fig. 12.1). This type of resource competition for a single limiting resource was named the $R^*$ rule (“R-star”):

In the absence of other processes, the species with the lowest $R^*$ for the shared limiting resource will exclude all others.

We refer to this type of resource competition as $R^*$ competition.

We can create a model to integrate the dynamics of these three components.

cr_2plants  <- function(t,y,p){
  N1 <- y[1]; N2 <- y[2]; R <- y[3]
  with(as.list(p),{
    dN1 <- N1*(e1*a1*R - d1)
    dN2 <- N2*(e2*a2*R - d2)
    dR <- r*(S-R) - a1*N1*R - a2*N2*R
    return(list(c(dN1,dN2,dR)))
  } )
}

Now we run the model.

t <- seq(0, 100, by=.5)
p <- list(e1=0.1, a1=1, d1=0.3, 
          e2=0.1, a2=.5, d2=0.1,
          S=5, r=1)

y.initial <- c(N1=.1, N2=.1, R=4)
rout <- ode(y=y.initial, time=t, fun=cr_2plants, parms=p)
r.long <- rout %>% as.data.frame %>% pivot_longer(cols=-time, names_to = "State", values_to = "Grams")
ggplot(data=r.long, aes(time, Grams, colour=State, linetype=State)) + 
  geom_line()

Now if we like, we can find $N_2^*$ when $R^*=d_2/(e_2a_2)$, and $N_1 =0$. \[\frac{dR}{dt} = 0 = r\left(S-R\right) - a_2N_2R\] Substituting $R^*$ in for $R$ gives us this solution:

\[N^*_2 = r\left(S\frac{e_2}{d_2} - \frac{1}{a_2}\right)\] It is clear from this expression that increasing the death rate ($d$) will drive down $N^*_2$. However, increasing resource turnover rate, $r$, assimilation efficiency, $e$, or uptake rate, $a$, will all cause $N_2^*$ to increase.

In the above model, per capita growth rates were linear functions of resource concentration. It is more common to see a type II functional response for the consumers (D. Tilman 1982) (Fig. 12.2), which we described in Chapter 8. \[\frac{dN_i}{Ndt}= b_i\frac{R}{k_i + R} - d_i\]

b <-c(1, .5); k=c(5, .1); d = 0.3
{curve(b[1]*x/(k[1]+x) - d, 0, 10, ylab=bquote(dP/(Pdt)), xlab="R")
curve(b[2]*x/(k[2] + x) - d, lty=2, add=TRUE)
abline(h=0, lty=3); text(7.5, 0, "gain = loss", adj=c(.5, -.3))
}

The Monod function for per capita resource-dependent growth rates of two different populations. These are a type II functional responses.

Figure 12.2: The Monod function for per capita resource-dependent growth rates of two different populations. These are a type II functional responses.

As we saw earlier with linear mass-specific growth rates, one species will have a lower $R^*$ than the other, and this model predicts that it will out-compete the other.

12.2 Competition for Substitutable resources

For two species, $N_1$ and $N_2$, competing for two shared and substitutable⁷⁷ resources, $R_1,\,R_2$, we will identify $e_{ij}$ and $a_{ij}$, where $e$ is assimilation efficiency, $a$ is feeding rate. They now have two subscripts, where $i$ refers to the consumer species $N_{i}$, and $j$ refers to either resource 1 or 2 ($R_1$, $R_2$). With that, we have the following relations.

\[\frac{dN_1}{dt}= N_1\left(e_{11} a_{11} R_1 + e_{12} a_{12} R_2 - d_1\right)\] \[\frac{dN_2}{dt}= N_2\left(e_{21} a_{21} R_1 + e_{22} a_{22} R_2 - d_2\right)\]

\[\frac{dR_1}{dt}= r_1\left(S-R_1\right) - a_{11}N_1R_1 - a_{21}N_2 R_1\] \[\frac{dR_2}{dt}= r_2\left(S-R_2\right) - a_{12}N_1R_2 - a_{22}N_2R_2\]

where $S$ is the background supply rate of a resource, and $r(S-R)$ is the turnover rate of resources that are part of the available $R$ pool.

The differential feeding rates, $a_{i1}$ and $a_{i2}$, tell us a lot about the relative requirements and limitations experienced by the two consumers. If \[a_{i1} > a_{i2}\] then, all else being equal, species $i$ requires more $R_1$ than $R_2$ in the environment.

Here we create an R function for the 4-D system of ODEs.

## 
cr_2substitutable <- function(t, y, parameters){
  N1 <- y[1] ; N2 <- y[2] ; R1 <- y[3]; R2 <- y[4]; 
  with(as.list(parameters), {
    N1.dot <- N1 * ( e11*a11*R1 + e12*a12*R2  - d1 )
    N2.dot <- N2 * ( e21*a21*R1 + e22*a22*R2  - d2 )
    R1.dot <- r1 * (S1 - R1) - (a11*N1 + a21*N2)*R1
    R2.dot <- r2 * (S2 - R2) - (a12*N1 + a22*N2)*R2
    return( list(c(N1.dot, N2.dot, R1.dot, R2.dot)) )
  } )
}

Next we create a parameter set in which species 1 grows more with each unit of $R_1$ and also consumes more of it than $R_2$.

parameters <- c(# the environment supplies slightly more of resources 1
                r1 = 1, r2=1, S1=1.8, S2=1.2, 
                # sp 1 grows more with each unit of R1
                # sp 2 grows more with each unit of R2
                e11=.1, e12=.1, e21=.1, e22=.1,
                # sp 1 consumes more R1 than R2 and more than sp 2
                # sp 2 consumes more R2 than R1 and more than sp 1
                a11=2, a12=1, a21=1, a22=2, 
                # both spp die at the same rate
                d1=.1, d2=.1)

t <- seq(0, 200, by=.1); y.init <- c(N1=1, N2=1, R1=1.8, R2=1.2)
out <- ode(y = y.init, times=t, func=cr_2substitutable, 
           parms=parameters)
outL <- out %>% as.data.frame %>% pivot_longer(-time, names_to="State_vars", values_to="abun.")
ggplot(outL, aes(time, abun., colour=State_vars)) + geom_line()

Initial dynamics of two species competing for two shared substitutable resources. Species 1 uses resource 1 more efficiently than it does resource 2 and also consumes more of it, per capita. Conversely, species 2 uses resource 2 more efficiently than it does resource 1 and also consumes more of it per capita.

Figure 12.3: Initial dynamics of two species competing for two shared substitutable resources. Species 1 uses resource 1 more efficiently than it does resource 2 and also consumes more of it, per capita. Conversely, species 2 uses resource 2 more efficiently than it does resource 1 and also consumes more of it per capita.

It is

output <- cr_2subst_graph(parameters)

12.3 Tilman’s resource ratio with essential resources

Coefficients reflect resources 1 and 2, and species $i$.

$a_i$ = relative resource production rate
$e_{ij}$

\[ \dot{R_1} = a_1\left(S_1-R_1\right) - \left( \dot{N_1} + m_1 \right) - \left( \dot{N_2} + m_2 \right)\] \[ \dot{R_2} = a_2\left(S_2-R_2\right) - \left( \dot{N_1} + m_1 \right) - \left( \dot{N_2} + m_2 \right) \]

\[\dot{N_1} = r_1 \,\mathrm{MIN}\left(e_{11}c_{11}\frac{R_1}{k_{11} + R_1} \, , \, e_{12}c_{21}\frac{R_2}{k_{21} + R_1}\right) - m_1\]

\[\dot{N_2} = r_2 \,\mathrm{MIN}\left(e_{21}c_{12}\frac{R_1}{k_{12} + R_1} \, , \, e_{22}c_{22}\frac{R_2}{k_{22} + R_1}\right) - m_2\]

## Something is not right here
resratio <- function(t, y, parameters){
  R1 <- y[1]; R2 <- y[2]; N1 <- y[3] ; N2 <- y[4] ;
  with(as.list(parameters), {
    N1.dot <- r1 * N1 * ( min( c( c11*R1/(k11 + R1), c21*R2/(k21 + R2) ) ) - m1 )
    N2.dot <- r2 * N2 * ( min( c( c12*R1/(k12 + R1), c22*R2/(k22 + R2) ) ) - m2 )
    ## R.dot must follow N.dot if we use N.dot in R.dot
    R1.dot <- a1 * (S1 - R1) - 
      r1 * N1 * ( min( c( c11*R1/(k11 + R1), c21*R2/(k21 + R2) ) ))/e11 - 
      r2 * N2 * ( min( c( c12*R1/(k12 + R1), c22*R2/(k22 + R2) ) ))/e12 
    R2.dot <- a2 * (S2 - R2) - 
      r1 * N1 * ( min( c( c11*R1/(k11 + R1), c21*R2/(k21 + R2) ) ))/e21 - 
      r2 * N2 * ( min( c( c12*R1/(k12 + R1), c22*R2/(k22 + R2) ) ))/e22 
    return( list(c(R1.dot, R2.dot, N1.dot, N2.dot)) )
  } )
}

parameters <- c(a1 = 1, a2=1, S1=10, S2=10, 
                r1=1, r2=1, m1=.1, m2=.1, 
                c11=1, c21=1, c12=1, c22=1, 
                k11=.1, k21=1, k12=1, k22=.1,
                e11=.1, e21=.1, e12=.1, e22=.1)

with(as.list(parameters),{
  curve(r1*c11*x/(k11 + x)-m1, 0, 3, xlab=expression("R"["j"]), 
        ylab=expression("r"["i"]*"c"["i"]*"R"["j"] /( "k"["ji"] + "R"["j"]) - "m"["i"] ) )
  curve(r1*c21*x/(k21 + x)-m1, 0, 3, add=T, lty=2)
  curve(r2*c12*x/(k12 + x)-m2, 0, 3, col=4, add=T)
  curve(r2*c22*x/(k22 + x)-m2, 0, 3, add=T, col=4, lty=2)
  abline(h=0, lty=3)
  legend('right', c('N1(R1)', 'N1(R2)', 'N2(R1)', 'N2(R2)'), 
         lty=c(1,2,1,2), col=c(1,1,4,4), lwd=2, bty='n')
})

t <- 0:100
y0 <- c(R1=.5, R2=.5, N1=10, N2=10)
run1 <- ode(y=y0, times=t, func=resratio, parms=parameters)
parameters; plot(run1)

##   a1   a2   S1   S2   r1   r2   m1   m2  c11  c21  c12  c22  k11  k21  k12  k22 
##  1.0  1.0 10.0 10.0  1.0  1.0  0.1  0.1  1.0  1.0  1.0  1.0  0.1  1.0  1.0  0.1 
##  e11  e21  e12  e22 
##  0.1  0.1  0.1  0.1

run1[101,]

##        time          R1          R2          N1          N2 
## 100.0000000   0.1111048   0.1111048   4.9447001   4.9447001

\[R_{j}^* = \frac{k_{ji}m_i}{r_i-m_i}\]

12.4 Consumer-resource mutualism

In this section, we will model mutualisms using a framework that, unlike our Lotka-Volterra approach, is intended specifically to describe resource consumption by one or both partners (R. H. MacArthur 1972).

Sometimes mutualists share resources, as with mycorrhizal associations, where species share a resource or a product that they can acquire or manufacture in surplus and trade it for a resource or product that is limiting. This is a bi-drectional mutualism (Fig. 12.4) because the flow of are bidirectional: both species provide a resource to the other (Holland and Deangelis 2010).

Bidirection, unidirectional, and indirect mutualisms. Normal arrowhead is a positive effect, whereas an open dot arrowhead is a negative effect. A dashed line indicates a service, and a solid line indicates a consumable resource.

Figure 12.4: Bidirection, unidirectional, and indirect mutualisms. Normal arrowhead is a positive effect, whereas an open dot arrowhead is a negative effect. A dashed line indicates a service, and a solid line indicates a consumable resource.

In other cases, mutualists may share services, as with moray eel-grouper foraging. In this particular case, grouper and the giant moray eel have different foraging strategies, with the eel going into reef crevices and the grouper waiting outside to pick off prey that escape the eel (Bshary et al. 2006). Sometimes these mutualisms are called by-product mutualisms, in cases where there is no evidence of coordinated behavior, but nonetheless a mutualism arises via the independent activities of each species.

In the following, we explore the type of mutualism in which one species provides a resource and the other species provides a service (Fig. 12.4). This is a uni-directional mutualism because the flow of resources is only unidrectional, from one species to the other (Holland and Deangelis 2010).

A good example of a uni-directional of mutualism is seed dispersal via a frugivore that disperses seeds. A plant species provides a fruit that is nutritious for an animal. A primary benefit for the plant is that seeds are carried away from the parent plant. This may help seeds escape enemies that accumulate near parent plants (Connell-Janzen hypothesis). In addition, it provides a bet hedging opportunity wherein seeds are dispersed to a variety of habitats any one of which may be better for survival than the current location. It is even possible that dispersal is directed to preferred habitat that is consistently better than under the parent plant due to better resource levels or environmental conditions (Wenny and Levey 1998). The animal benefits by consuming the fleshy part of the fruit surrounding the seed. An interesting case arises when seed predators collect seeds but do not kill them before having buried, cached, or moved seeds.

Another very important example example of this type of mutualism is pollination. Pollinators (bats, insects, birds) receive resources (nectar, pollen) and provide services (pollination, gene flow).

12.4.1 A model of uni-directional mutualism

Here we describe a model of a facultative unidirectional mutualism that uses type II foraging model for both consumption of the resource and provision of the service.

The growth equation for the plant species, $M_1$, includes logistic growth, and has terms for the benefit it receives from its mutualist frugivore and the cost paid to that frugivore.

\[\frac{dM_1}{dt} = r_1 M_1\left(1-d_1M_1\right) - aM_2\left(\frac{M_1}{h_1 + M_1}\right) + sM_1\left(\frac{M_2}{h_2 + M_2}\right)\]

The plant species ($M_1$) grows logistically in the absence of fruit consumption or seed dispersal $r_1 M_1\left(1-d_1M_1\right)$. Consumption of fruits are governed by a type II functional response by the frugivore $aM_2\left(\frac{M_1}{h_1 + M_1}\right)$, so that the consumption rate of a single frugivore has an upper limit (Fig. 12.5). Successful seed dispersal and establishment is also governed by a type II functional response. The success of a single seed reaches an upper limit as the number of frugivores continues to increase, $c_1M_1\left(\frac{M_2}{h_2 + M_2}\right)$ (Fig. 12.5).

ggplot(data.frame(x=c(0, 1)), aes(x=x) ) + 
  stat_function(fun=function(x) { 10*x/(1+10*x) }) +
  labs( x="Fruit abundance", y="Consumption Rate of a frugivore")


ggplot(data.frame(x=c(0, 1)), aes(x=x) ) + 
  stat_function(fun=function(x) { 10*x/(1+10*x) }) +
  labs( x="Frugivore abundance", y="Dispersal rate of a seed")

Fruit consumption by a single frugivore is governed by a type II functional response as fruit abundance increases. Successful seed dispersal and establishment is governed by a type II functional response as frugivores increase.

Figure 12.5: Fruit consumption by a single frugivore is governed by a type II functional response as fruit abundance increases. Successful seed dispersal and establishment is governed by a type II functional response as frugivores increase.

The growth equation for the frugivore species, $M_2$, includes logistic growth, and one term for the benefit gained from the plant.

\[\frac{dM_2}{dt} = r_2M_2\left(1-d_2M_2\right) + eaM_2\left(\frac{M_1}{h_1+M_1}\right) \] The animal seed disperser grows logistically in in the absence of consuming this particular plant species $r_2M_2\left(1-d_2M_2\right)$. The animal benefits from consuming the fruit which it attacks at the same rate as for $M_1$, $aM_2\left(\frac{M_1}{h_1+M_1}\right)$, but converts the fruit into new individuals with efficiency $e$.

consumption > service = net predation
service > consumption = net mutualism

Here we provide code for the ODE for the complete bidirectional mutualism. When we define the parameters, we will set some of them equal to zero so that it represents a unidirectional consumer-resource mutualism.

# The derivative 
## A two-species consumer-resource mutualism
## Holland and DeAngelis 2009, Ecology Letters Ecology Letters, 12: 1357–1366.
## The parameterization follows that in this article. It differs from those 
## in the above equations. Can you see the differences?

cr_bimut <- function(t,y,parameters) {
  M1 <- y[1] # 
  M2 <- y[2] # 
  with( as.list( parameters ), {
    dM1 <- M1 * (r1 - d1*M1 + a12*M2/(h2 + M2) - B1*M2/(e1+M1) )
    dM2 <- M2 * (r2 - d2*M2 + a21*M1/(h1 + M1) - B2*M1/(e2+M2) )
    list( c( dM1, dM2 ))
  }
  )
}

And here are the parameters, and we graph the time series.

### These parameters create a unidirectional, resource-service mutualism (B2 = 0)
## Translating, a12 = a above, and a21 = ea above, and B1 = s above
## therefore e = a21/a12

p.u <- list(r1 = 1, r2 = 1, d1=0.01, d2=0.01, 
            a12 = 0.4, a21 = 0.25, B1 = 0.3, B2 = 0, 
            h1=0.3, h2=0.3, e1=0.3, e2=0.3)

#####
# times series
t <- seq(0,25, by=.1)
y <- c(M1=50, M2=1)

out <- ode(y=y, times=t, func= cr_bimut, parms=p.u)
outdf <- data.frame(out)
outL <- pivot_longer(outdf, -time, 
                     names_to="State_var", values_to="N")
ggplot(outL, aes(time, N, colour=State_var)) + geom_line()

The dynamics can result in multiple basins of attraction or alternative stable states (Fig. 12.6).

A resource-service mutualism may result in alternative stable states. If both species achieve moderate abundance, then both species increase toward a stable mutualism. However, if the frugivore is too abundant relative to the plant, then its consumption dominates the interaction and drives the equilibrium mutualism.

Figure 12.6: A resource-service mutualism may result in alternative stable states. If both species achieve moderate abundance, then both species increase toward a stable mutualism. However, if the frugivore is too abundant relative to the plant, then its consumption dominates the interaction and drives the equilibrium mutualism.

Holland and DeAngelis (2009) vary $a_{21}$ between $(0.25,\ldots,\,0.4)$ to vary the relation from mutualism to predation. Try that now. Save a picture of your results.

12.5 Plant-soil feedbacks

Vascular plants form the basis of terrestrial ecosystems, and they live in intimate contact with the most diverse group of organisms on Earth–soil microbes. Ecologists have begun to focus on this intimate relationship, and Jim Bever (Bever 2003) has been one of ecologists helping lead the way (Fig. 12.7).

Plants and soil microbes can interact in negative or positive ways. Plants may grow more poorly with their own soil flora than with that of a competitor's.

Figure 12.7: Plants and soil microbes can interact in negative or positive ways. Plants may grow more poorly with their own soil flora than with that of a competitor’s.

Fig. 12.7 reflects the following Lotka-Volterra style model of plant species $N_A$ and $N_B$, and soil microbial floras $S_A$ and $S_B$: \[\begin{equation} \frac{dN_A}{dt} = r_A N_A \left( 1 + \alpha_A S_A + \beta_A S_B - \frac{N_A + c_B N_B}{K_A}\right) \end{equation}\] where the effects of microbes ($\alpha$, $\beta$) could be positive or negative (e.g. pathogenic). The plants each have their own carrying capacities ($K$), and negative effects on each other ($c$).

The corresponding equation for the other plant species is \[\begin{equation} \frac{dN_B}{dt} = r_B N_B \left( 1 + \alpha_B S_A + \beta_B S_B - \frac{c_A N_A + N_B}{K_B}\right). \end{equation}\]

We often refer to the soil flora associated with a single plant species as its home flora, and the other as the away flora. Each microbial flora comprises many, many species and we aggregate the net effects. Bever (2003) further simplified his model by placing constraints on soil microbes, so that $S_A + S_B=1$. This allows him to focus on the relations among microbes and plants. With these assumptions, he showed that the net effect of each flora as a single state variable for the home flora of each species: \[\frac{dS_A}{dt} = S_A(1-S_A) \left( \frac{N_A}{N_A + N_B} - \nu\frac{N_B}{N_A + N_B}\right)\] and that \[S_B=1-S_A\].

Coexistence criteria Bever states that (2003, p. 467, bottom of the second column) out that the interactions are more stable if \[\alpha_A + \beta_B < \alpha_B + \beta_A \] which means that the effects of the home floras own their home plants are less positive (or more negative) then on the other species.

Bever also argued that for either plant species to increase when rare (i.e. coexist), that the interspecific competition has to be less than the benefits of the interactions of the soil floras. However, his criterion in eqn 5 (p. 469) is not sufficient, and this is revealed by its lack of consistency with the parameter set and outcomes in Fig. 4. So, we will alter his equation a bit to reveal what probably would work. We use a particular parameterization of Lotka-Volterra interactions so that it is consistent with others (B. M. Bolker, Pacala, and Neuhauser 2003; Chesson 2000b). We reformulate the carrying capacities and $c_A$ and $c_B$ as

\[\alpha_{AA}=1/K_A;\,\alpha_{BB}=K_B;\,\alpha_{BA}= c_A/K_B;\,\alpha_{AB}=c_B/K_A\] If we recall our earlier work, we showed that each species $i$ can invade when \[\frac{\alpha_{ij}}{\alpha_{jj}} < 0\] This translates to the Lotka-Volterra invasion criterion for species A is \[\frac{\alpha_{AB}}{\alpha_{BB}} =\frac{c_b}{K_A}K_B\] The invasion criterion for species B is \[\frac{\alpha_{BA}}{\alpha_{AA}} =\frac{c_a}{K_B}K_A\] The coexistence criterion is that each species must be able to invade when rare, including when the soil flora has no effect. Therefore, we have

\[\begin{equation} \left(\frac{c_b}{K_A}K_B,\, \frac{c_a}{K_B}K_A\right) < \frac{(1+\alpha_B)(1+\beta_A)}{(1+\alpha_A)(1+\beta_B)} \tag{12.1} \end{equation}\]

where both conditions on the left must be satisfied. In the numerator on the right, we have the effect of each soil flora on their away plant species. In the denominator, we have the effects of each flora on the home plant species. Thus, the fraction on the right is the relative benefit of the soil flora.

Now we will create an ODE function that will let us see the dynamics play out. We need only model one soil flora because $S_B = 1-S_A$.

bever3 <- function(t,y,p){
  with(as.list(c(y,p)), {
    dNa <- ra * Na * (1+alphaA*SA + betaA*(1-SA) - (Na + cB*Nb)/Ka)
    dNb <- rb * Nb * (1+alphaB*SA + betaB*(1-SA) - (cA*Na + Nb)/Kb)
    dSA <- SA*(1-SA) * ( Na/(Na+Nb) - nu*Nb/(Na+Nb) )
    return(list(c(dNa, dNb, dSA), SB=1-SA))
  })
  }

We will parameterize the model in a way that reveals the importance of the soil floras to coexistence, according to criteria stated above (12.1). We will create a vector of parameters, and then test whether the inequality in (12.1) is true. We start with an absence of a soil flora effect.

p0 <- list(ra=0.7, rb=0.5, Ka = 10, Kb=12, cA=.98, cB=.98, 
           alphaA=0, alphaB=0, betaA=0, betaB=0, nu=0.8)

Use the criterion above, and paper and a pencil, and see what you discover. We could also do it in R.

with(p0, { cA/Kb*Ka  < (1+alphaB)*(1+betaA) / ( (1+alphaA)*(1+betaB)) })

## [1] TRUE

with(p0, { cB/Ka*Kb < (1+alphaB)*(1+betaA) / ( (1+alphaA)*(1+betaB)) })

## [1] FALSE

                                                                            This shows us we will not get coexistence in the absence of soil flora. Now we add the soil flora.

p <- list(ra=0.7, rb=0.5, Ka = 10, Kb=12, cA=.98, cB=.98, alphaA=-0.03, alphaB=0.1, betaA=0.1, betaB=-0.2, nu=0.8)

Try your hand again, with paper and pencil. And then do it in R.

with(p, { cA/Kb*Ka  < (1+alphaB)*(1+betaA) / ( (1+alphaA)*(1+betaB)) })

## [1] TRUE

with(p, { cB/Ka*Kb < (1+alphaB)*(1+betaA) / ( (1+alphaA)*(1+betaB)) })

## [1] TRUE

On this basis, we should see that each species can invade when rare. Now let’s see the dynamics play out.

Figure 12.8: Dynamics of competing plant species, subject to self-limiting negative feedbacks with soil pathogens.

12.6 Simulations for learning

Often, we learn a lot by describing interactions in simple enough terms that we can analyze the interaction analytically, such as the way we analyzed Lotka-Volterra competition and mutualism. In more complex cases, we often cannot find an analytical, or closed form, solution, or if we do, it is so complex that it is difficult to extract meaning. In those cases, we may be able to learn a lot about the consequences of our model’s assumptions through simulation. Using simulation, we examine the observed dynamics for a wide range of parameter values and then infer general rules about how parameters or parameter combinations influence dynamics.

In this study, we will use simulations to attempt to determine rules that govern the long-term dynamics of the mutualism model with linear density dependence in both intra- and interspecific interactions. We will also use simulation to confirm our earlier analytical solution and vice versa.

12.6.1 `lvg()`

The core of our simulations will use a general purpose Lotka-Volterra model that simpy takes an interaction matrix and $r$ as inputs.

lvg <- function (time, n, parameters) 
{
  r <- parameters[[1]]
  a <- parameters[[2]]
  dns.dt <- r * n * (1 + (a %*% n))
  return(list(c(dns.dt)))
}

In this model, the parameters are a list (a type of R data object) with two components. The first is a vector of r, one for each species we want to model. The second component of the list is the matrix of interaction coefficients. This matrix is square (same number of rows and columns). The elements of this matrix must have the correct sign. For instance, to represent competition, the relevant elements would have to be negative.

Here we create an interaction matrix for a two-species Lotka-Volterra mutualism, where each species suffers from intraspecific competition, and benefits from its mutualist.

a <- matrix(c( -.1, .02,
               .02, -.09), nrow = 2, byrow=TRUE)
a

##       [,1]  [,2]
## [1,] -0.10  0.02
## [2,]  0.02 -0.09

Now let’s numerically integrate, or simulate, the model and plot the outcome. Comment your code to help remind you what is going on.

# Let r be the same fpor each species
r <- c(1,1)
parms <- list(r=r, a=a)
parms

## $r
## [1] 1 1
## 
## $a
##       [,1]  [,2]
## [1,] -0.10  0.02
## [2,]  0.02 -0.09

t=seq(0,30, by=.1)
N0 <- c(N1=20, N2=5)
lvout <- as.data.frame( ode(N0, t, lvg, parms) )
lvL <- pivot_longer(lvout, -time, names_to="State_vars", 
            values_to="N")

ggplot(lvL, aes(time, N, colour=State_vars, linetype=State_vars)) + 
  geom_line() + 
  geom_hline(yintercept=c(-1/a[1,1], -1/a[2,2]), linetype=c(1,3)) +
  annotate(geom="text", 
           x=c(15, 20),y=c(-1/a[1,1], -1/a[2,2]), label=c("K1", "K2"), vjust=-0.2 )

Figure 12.9: Dyanmics of a Lotka-Volterra mutualism.

12.6.2 Systematic simulations

When we use a model to explore the consequences of our assumptions, we want to use a range of plausible assumptions. In a case like this one where we are exploring a model more than the natural phenomenon it represents, we should use a wide range of parameter values that include cases where it is likely to fail. In our case, we will hold intraspecific density dependence constant and vary interspecific interaction coefficients from zero to twice the value of intraspecific interactions.

We start with building the bits and pieces that the simulation will use, and then use a for-loop to cycle through all the parameter combinations we are interested in.

# Use expand.grid to create systematically combinations of parameters
df <- expand.grid(a11 = -0.1, a12=seq(0, .2, by =.02),
                  a21 = seq(0, .2, by=.02),  a22=-0.1 )
dim(df) # The number of rows and columns

## [1] 121   4

head(df) # the top of the data frame

##    a11  a12 a21  a22
## 1 -0.1 0.00   0 -0.1
## 2 -0.1 0.02   0 -0.1
## 3 -0.1 0.04   0 -0.1
## 4 -0.1 0.06   0 -0.1
## 5 -0.1 0.08   0 -0.1
## 6 -0.1 0.10   0 -0.1

# set intrinsic rates of increase...
r <- c(1,1)

# and starting values...
N0 <- c(N1 = 10, N2= 10.001)

# and times for which to return output.
t <- 0:200

# Set values that will be useful later on in our for-loop.
n <- length( t + 1 )
npar <- nrow( df ) 

# create an empty data frame to hold our simulation output.
output <- data.frame( t = numeric(npar), N1 = numeric( npar ), 
                      N2 = numeric( npar ) )

### USE A  FOR-LOOP TO INTEGRATE THE MODEL FOR EACH COMBINATION OF PARAMETERS
for( i in 1:npar ) {
  coefs.true <- matrix( as.numeric( df[i,] ), 
                        nrow = 2, byrow = TRUE)
  invisible(
    capture.output( 
    # capture.output is a function that allows us to avoid seeing warnings...
    # not usually a good idea to ignore warnings!!!
    # but it is OK here....
    lvout <- ode(N0, t, lvg, parms = list(r, coefs.true) )
  ) )
  
  # hang on to only the last row of output, and  
  # put it in the i-th row of our output data frame.
  last <- nrow(lvout)
   output[i,] <- lvout[last,]
}

# Combine columns of parameters and output.
simdat <- cbind(df, output)

12.6.3 Checking simulation results

Let’s look for anomalous results, such as negative population sizes, or infinite or “too” large sizes.

# Do simplistic numerical summaries.
summary(simdat[,5:7])

##        t                  N1                   N2            
##  Min.   :  0.6931   Min.   : 1.000e+01   Min.   : 1.000e+01  
##  1st Qu.:  1.5576   1st Qu.: 2.100e+01   1st Qu.: 2.100e+01  
##  Median :200.0000   Median : 7.000e+01   Median : 7.000e+01  
##  Mean   :122.9019   Mean   :3.355e+146   Mean   :3.355e+146  
##  3rd Qu.:200.0000   3rd Qu.:8.858e+145   3rd Qu.:8.858e+145  
##  Max.   :200.0000   Max.   :8.187e+147   Max.   :8.187e+147

If we run this, we see some strange population sizes and time was not always 200, showing that the numerical integration failed before completing. Also, $N$ included unreasonably large values.

Let’s look for all rows in which either $N_1$ or $N_2$ are huge ($> 10\,000$).

# the vertical bar | means "or"
sim.N <- subset(simdat, N1 > 10000 | N2 > 10000)
head(sim.N)

##     a11  a12  a21  a22        t            N1            N2
## 43 -0.1 0.18 0.06 -0.1 3.241398 2.434044e+103 1.390882e+103
## 44 -0.1 0.20 0.06 -0.1 2.399838  1.164783e+99  6.212179e+98
## 52 -0.1 0.14 0.08 -0.1 2.889456  1.856660e+65  1.392495e+65
## 53 -0.1 0.16 0.08 -0.1 2.139629 1.175804e+109 8.140182e+108
## 54 -0.1 0.18 0.08 -0.1 1.774181 9.054385e+146 5.820676e+146
## 55 -0.1 0.20 0.08 -0.1 1.541959 1.041579e+147 6.249472e+146

Notice also that values for $t$ indicate that the integration couldn’t find values beyond a certain time point. Hmm. For now, let’s assume our simulation is doing what we think it is and that these results indicate a true lack of equilibrium.

What could we do to learn the emerging rules? What if we just plot the values of $\alpha_{ij},\,\alpha_{ji}$ where $N$ is really, really big?

ggplot(sim.N, aes(a12, a21)) + geom_point()

$Values of $\alpha$ when N is extraordinarily large.*$

Figure 12.10: Values of $\alpha$ when N is extraordinarily large.*

What do we learn from just this figure? Write down a possible set of conditions necessary to get these anomalous results.

12.6.4 Comparing simulation results to analytical solutions

Any time we do something novel in our modeling, we need to find a way to evaluate the novel approach to make sure it is doing what we think it is. When we have an analytical solution for our model, it is great to compare the analytical solution to the simulations.

One thing we can do is to compare the long-term steady state values we get in simulation to the analytical equilibria of this mutualism model.

First, we find analytical solutions for all of the equilibria predicted by the coefficients. Consult the equilibrium solution in the first section of the chapter. We then compare these to numerical solutions.

Here we calculate the predicted equilibria for each row of our coefficients, using the above solution.

denom <- simdat[,1] * simdat[,4] - simdat[,2] * simdat[,3]
# note this solution is subtracting the negative intraspecific comp coef.
# making it consistent with our equations in section 1
num1 <- simdat[,2] - simdat[,4]
num2 <- simdat[,3] - simdat[,1]
N1.eq <- num1 / denom
N2.eq <-num2 / denom

# combine with simulation data set
sim.check1 <- cbind(simdat, N1.eq, N2.eq)
# sim.check1[order(sim.check1$a12, sim.check1$a21),]

Let’s find the subset with potential problems, such as for $N_1$ or $N_2$ less than than zero or really, really big ($>10000$).

problems1 <- subset(sim.check1, (N1 < 0 | N2 < 0 | N1 > 1e4 | N2 > 1e4 ) )
# look at the subset
ggplot(problems1, aes(N1, N1.eq)) + geom_point()

Figure 12.11: Analytically calculated predictions for $N^*$ vs. observed simulated $N$.

Fascinating! Both the analytical and the simulated results give us wacky, but very different results. Let’s examine the coefficients that generated these results.

problems1[1:5,]

##     a11  a12  a21  a22        t            N1            N2      N1.eq
## 43 -0.1 0.18 0.06 -0.1 3.241398 2.434044e+103 1.390882e+103 -350.00000
## 44 -0.1 0.20 0.06 -0.1 2.399838  1.164783e+99  6.212179e+98 -150.00000
## 52 -0.1 0.14 0.08 -0.1 2.889456  1.856660e+65  1.392495e+65 -200.00000
## 53 -0.1 0.16 0.08 -0.1 2.139629 1.175804e+109 8.140182e+108  -92.85714
## 54 -0.1 0.18 0.08 -0.1 1.774181 9.054385e+146 5.820676e+146  -63.63636
##         N2.eq
## 43 -200.00000
## 44  -80.00000
## 52 -150.00000
## 53  -64.28571
## 54  -40.90909

What do you notice about these coefficients?

What would you get for each term in the for the analytical equilibrium using the first row of data?
How do you reconcile these results?

Once we understand the anomalous results above, let’s compare the simulated and the analytical solutions for $N_1$ and $N_2$ both greater than zero and both less than a huge number.

check2 <- subset(sim.check1, (N1.eq > 0 & N2.eq > 0) & (N1 < 1e4 & N2 < 1e4 ))

Next we calculate the differences between the analytical solutions and the simulated values.

N1.diff <- check2$N1 - check2$N1.eq
N2.diff <- check2$N2 - check2$N2.eq
summary(N1.diff); summary(N2.diff)

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -1.421e-14  1.776e-15  3.553e-15  1.363e-13  7.105e-15  3.183e-12

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -2.842e-14  1.776e-15  3.553e-15  1.558e-13  7.105e-15  5.230e-12

These summaries show us that the differences between the analytical solutions and the numerical approximations are very, very, very small. That is a good thing. It gives us confidence that our simulations are doing what we think they are doing.

So, we can now tentatively assume that our simulations are doing what we think they are, and we could have confidence with results that we generate.

In practice, sometimes we simulate a very complex model and compare the results with a simpler, tractable analytical solution. Here, we get the hang of the approach by comparing simulations against known solutions.

12.6.5 One last thing: using the Jacobian matrix

Use eigenanalysis of the Jacobian matrix to assess the stability of these equilibria. Consult eq. 5.21 and Table 5.2 of Stevens (2009).

A sensible equilibrium.

subset(sim.check1, a12 == .02 & a21 == .02)

##     a11  a12  a21  a22   t   N1   N2 N1.eq N2.eq
## 13 -0.1 0.02 0.02 -0.1 200 12.5 12.5  12.5  12.5

# form remember that we set r = 1
J1 <- matrix(c(-.1*12.5, -.02*12.5,
               -.02*12.5, -.1*12.5), nrow=2, byrow=TRUE)
eigen(J1)

## eigen() decomposition
## $values
## [1] -1.0 -1.5
## 
## $vectors
##            [,1]      [,2]
## [1,] -0.7071068 0.7071068
## [2,]  0.7071068 0.7071068

What do these results mean?

A nonsensical equilibrium

subset(sim.check1, a12 == .2 & a21 == .2)

##      a11 a12 a21  a22         t            N1            N2 N1.eq N2.eq
## 121 -0.1 0.2 0.2 -0.1 0.6931166 8.389656e+145 8.389656e+145   -10   -10

J2 <- matrix( c( -.1*(-10), -.2*(-10), 
                 -.2*(-10), -.1*(-10)), nrow=2, byrow=TRUE) 
eigen(J2)

## eigen() decomposition
## $values
## [1]  3 -1
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068