8 Mutualisms

8.1 Background

Mutualisms occur in a wide variety of different species. Mutualisms form the foundations of ecosystems including coral reefs (coral and zooxanthellae), grasslands, and forests (plants and mycorrhizal fungi, pollinators, and dispersers). The natural history of mutualisms are wildly diverse and complex, and and Boucher, James, and Keeler (1982) list types of mutualisms:

Energetic
Nutritional
Defense
Transport
Obligate vs. facultative
Direct vs. indirect
Symbiotic vs. independent

One of the most common of these is indirect mutualism (Fig. 8.2). In this case, a series of negative direct interactions can yield net positive indirect interactions, as in the classic case where ``the enemy of my enemy is my friend.’’ For instance, consider a three plant species that all compete for the same limiting resource. When species A suppresses species B, it is indirectly helping species C (Miller 1994).

We will focus on direct mutualisms, and in this chapter, we’ll explore two approaches to modeling mutualisms, Lotka-Volterra models and consumer-resource models.

8.2 Lotka-Volterra mutualism

Lotka-Volterra mutualism is very straightforward, following the same template as competition. All that we have to do is change the sign associated with the other species. As usual, \(\alpha_{ij}\) is the per capita effect of species \(j\) on species \(i\).

\[\begin{align*} \frac{dN_1}{dt} &= r_1N_1\left( 1 - \alpha_{11}N_1 + \alpha_{12}N_2\right)\\ \frac{dN_2}{dt} &= r_2N_2 \left(1 + \alpha_{21}N_1 - \alpha_{22}N_2\right) \end{align*}\]

In the above equations we see the sign of the intraspecific terms \(\alpha_{ii}N_i\) is negative, while the sign of the interspecific terms \(\alpha_{ij}N_j\).

The zero net growth isoclines are found the same way we found those for interspecific competition, by solving \(dN_1/dt\) (for zero). We will graph these, so we solve each in terms of species 1. \[\begin{align*} 0 &= 1 + \alpha_{21}N_1 - \alpha_{22}N_2\\ N_2 &= 1/\alpha_{22} + \frac{\alpha_{21}}{\alpha_{22}}N_1 \end{align*}\] \[\begin{align*} 0 &= 1 - \alpha_{11}N_1 + \alpha_{12}N_2\\ N_2 &=\frac{\alpha_{11}}{\alpha_{12}}N_1 - 1/\alpha_{12} \end{align*}\]

The equilibria are \[N_i^* = \frac{\alpha_{jj} + \alpha_{ij}}{\alpha_{ii}\alpha_{jj} - \alpha_{ij}\alpha_{ji}}\] Note this assumes that both \(\alpha_{jj}\) and \(\alpha_{ji}\), per se, are greater than zero, because our equation assumed it was so.

Let’s use this to describe a facultative mutualism, that is, one in which neither species needs the other to persist, but in which they each do better. What we mean by that is that let \(r>0\). For now, we will also assume that the mutualism provides only a modest benefit, where species benefits to each is smaller than the negative effects on themselves (\(\alpha_{ii} > \alpha_{ij}\)). The dynamics are stable (Fig. 8.1)

To illustrate the dynamics of this population, we will parameterize an ODE model. Our parameters meet our assumptions stated above.

parameters <- c(alpha11 = .2, alpha22 = .2, alpha12 = .1, alpha21=.1,
                r1=.01, r2=.01)

Here we write an ODE model to help anaylze dynamics.

mutualism.LV <- function(t, y, params){
  n1 <- y[1]; n2 <- y[2]
  with(as.list(params), {
    dn1.dt = r1*n1*(1 - alpha11*n1 + alpha12*n2)
    dn2.dt = r2*n2*(1 + alpha21*n1 - alpha22*n2)
    return(list( c(dn1.dt, dn2.dt) ) )
  })
}

The Lotka-Volterra approach predicts exclusion of one species in two realistic scenarios, obligate mutualisms and mutualisms with large effects. An obligate mutualism is one in which growth rate is negative in the absence of the mutualism. The way to describe that in Lotka-Volterra models is to set \(r < 0\). This results in an unstable equilibrium which is a global repellor which results in the loss of one species. Second, if the benefit is greater than the negative intraspecific effect, we get an uncontrolled positive feedback loop leading to infinite population sizes. In this sense, it is like exponential growth: a potentially very useful idea, but it must be moderated somehow in order to reflect long term trajectories.

Only weak facultative mutualisms are stable.

Figure 8.1: Only weak facultative mutualisms are stable.

8.3 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 (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. 8.2) 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 8.2: 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. 8.2). 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).

8.3.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. 8.3). 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. 8.3).

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 8.3: 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. 8.4).

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 8.4: 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.

8.4 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. 8.5).

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 8.5: 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. 8.5 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 the parameterization from Chapter ?? so that it is consistent with others (Bolker, Pacala, and Neuhauser 2003; Chesson 2000 b). 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{8.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 (8.1). We will create a vector of parameters, and then test whether the inequality in (8.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 8.6: Dynamics of competing plant species, subject to self-limiting negative feedbacks with soil pathogens.

8.5 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.

8.5.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 8.7: Dyanmics of a Lotka-Volterra mutualism.

8.5.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)

8.5.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()

Figure 8.8: 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.

8.5.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 8.9: 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      N2.eq
## 43 -0.1 0.18 0.06 -0.1 3.241398 2.434044e+103 1.390882e+103 -350.00000 -200.00000
## 44 -0.1 0.20 0.06 -0.1 2.399838  1.164783e+99  6.212179e+98 -150.00000  -80.00000
## 52 -0.1 0.14 0.08 -0.1 2.889456  1.856660e+65  1.392495e+65 -200.00000 -150.00000
## 53 -0.1 0.16 0.08 -0.1 2.139629 1.175804e+109 8.140182e+108  -92.85714  -64.28571
## 54 -0.1 0.18 0.08 -0.1 1.774181 9.054385e+146 5.820676e+146  -63.63636  -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.

8.5.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

What do these results mean?

8.5.6 In Fine

So, how can you digest and communicate what you have found? You should be able to answer these questions:

Did the simulation perform as you thought and how do you know?
Did the simulation confirm something we first assessed analytically? If so, how?
How do you explain the differences between the simulations and analytical equilibrium (extremely large vs. negative population sizes)?

References

Bever, James D. 2003. “Soil community feedback and the coexistence of competitors: Conceptual frameworks and empirical tests.” New Phytologist 157 (3): 465–73.

Bolker, B M, S W Pacala, and C Neuhauser. 2003. “Spatial dynamics in model plant communities: {W}hat do we really know?” American Naturalist 162: 135–48.

Boucher, D. H., S. James, and K. H. Keeler. 1982. “The ecology of mutualism.” Annual Review of Ecology and Systematics 13: 315–47.

Bshary, Redouan, Andrea Hohner, Karim Ait-el-Djoudi, and Hans Fricke. 2006. “Interspecific communicative and coordinated hunting between groupers and giant moray eels in the red sea.” PLoS Biology 4 (12): 2393–8.

Chesson, P L. 2000b. “Mechanisms of maintenance of species diversity.” Annual Review of Ecology and Systematics 31: 343–66.

Holland, J. Nathaniel, and Donald L. Deangelis. 2010. “A consumer-resource approach to the density-dependent population dynamics of mutualism.” Ecology 91 (5): 1286–95.

MacArthur, R H. 1972. Geographical Ecology: Patterns in the Distribution of Species. New York: Harper & Row.

Miller, T E. 1994. “Direct and indirect interactions in an early old-field plant community.” The American Naturalist 143: 1007–25.

Wenny, D G, and D J Levey. 1998. “Directed seed dispersal by bellbirds in a tropical cloud forest.” Proceedings of the National Academy of Sciences of the United States of America 95 (11): 6204–7.