# Predator-Prey Interactions

## Agent-Based Models

Explore the agent-based model below:

::: callout-tip
## Agent-based Models

Compared to other simulation types that attempt to simulate systems as a whole, agent-based models are a type of simulation where the behavior of individual units (agents) are simulated, to observe the additive result of how individual behaviors and/or interactions might affect the system as a whole.
:::

```{=html}
<iframe src="Wolf_Sheep_Grass_Predation.html" data-external="1" style="width:100%; width:100vw; width:100dvw; height:100%; height:100vh; height:100dvh;"></iframe>
```
```{html include = FALSE}
[Link to ABM Predation Model](https://charleslehnen.github.io/BISC_499_Ecology_and_Biodiversity/Lab_4/Wolf_Sheep_Grass_Predation.html)
```

Select the model version `sheep-wolves`. Click `setup` to set-up the model based on your parameters. Use the `go` button to start and stop the simulation.

<details>

<summary>What do you notice about the synchronicity between the sheep and wolf populations?</summary>

This is seen throughout predator-prey dynamics, including the famous lynx--snowshoe hare cycle, based on data from the Hudson Bay Trading.

![](images/image-1335811125.png)

</details>

<details>

<summary>Challenge: Modify the parameters to stabilize the system so that neither population becomes extinct. What do you find?</summary>

</details>

<details>

<summary>

Change the model version to `sheep-wolves-grass`. Before running the simulation, what do you expect to happen?

Now run the simulation and modify the parameters to stabilize the system. What do you find?

</summary>

</details>

<details>

<summary>Why is this the case?</summary>

</details>

------------------------------------------------------------------------

## Lotka-Volterra Models

In R, we can construct a simple Lotka-Volterra model to simulate single predator/prey interactions.

$$
 \begin{array}{l}
\frac{dN}{dt} \ =\ rN\ -\ aNP\\
\\
\frac{dP}{dt} \ =\ faNP\ -\ qP
\end{array}
$$

<details>

<summary>What do the terms for the above equation refer to?</summary>

$$
 \begin{array}{l}
r\ =\ prey\ per\ capita\ growth\ rate\\
N\ =\ number\ of\ prey\ species\ individuals\\
a\ =\ per\ capita\ attack\ rate\\
P\ =\ number\ of\ predator\ species\ individuals\\
f\ =\ predator\ efficiency\\
q\ =\ predator\ per\ capita\ mortality\ rate
\end{array}
$$

Remember from last lab, $r$ is generally used for continuously reproducing populations

</details>

<details>

<summary>In the absence of a predator, would the terms for prey population growth be density dependent or independent?</summary>

</details>

Now let's create a function for the Lotka-Volterra model in R. We can write differential equations in R using the ordinary differential equation (ODE) package `deSolve`[@deSolve] following the format:

``` r
function_name <- function(y, x, parameters) {
  X <- x[1]
  with(as.list(params), {
    dX.dy <- X # complete formula here
    return(list(c(dX.dy)))
  })
}
```

<details>

<summary>

R's ordinary differential equation (ODE) solver deSolve() can be used to analyze population growth. ODEs are commonly used to calculate rates of change in order to model dynamic systems.

Now you try to use this format to define an ODE function for the prey (`dN.dt`) and predator (`dP.dt`) Lotka-Volterra models named `lotka_volterra`

</summary>

```         
```

</details>

Now let's run our model and observe the outcome:

```{r}

library(deSolve)
library(ggplot2)

parms <- c(
  r = 0.5,
  a = 0.01,
  f = 0.1,
  q = 0.2  
)

Time <- seq(0, 200, by = 0.1)

lotka_volterra.out <- deSolve::ode(c(N0 = 25, P0 = 5), times = Time, func = lotka_volterra, parms = parms)

df <- data.frame(Time, lotka_volterra.out)

ggplot(df, aes(x = Time)) +
  geom_line(aes(y = P0), color = "red") +
  geom_line(aes(y = N0), color = "black") +
  labs(x = "Time", y = "Population Size")

```

::: callout-tip
## Namespaces in R

In the above example we call an example using this format `package::function().`

Namespaces are used to organize functions, data, and other objects within a package. The double colon operator allows us to explicitly access a namespace within a specified package.

We are implicitly accomplishing this same task whenever we run a package function or access a package dataset after loading that package using `library()` or `require()`. The double colon allows us to do this same tasks explicitly either to avoid conflicts with a function/dataset from another package with the same name (as above) or to be more explicit programmers.
:::

<details>

<summary>What do you notice about the plot?</summary>

</details>

<details>

<summary>Now take some time to alter the parameters, including the starting population size, and observe what happens. Also note if you are able to have populations fully stabilize to a node or if they continue to oscillate?</summary>

</details>

### Functional Responses

The term $aN$ describes the functional response of the predator species, the rate at which a prey kills/consumes its prey (population size $N$) at an attack rate $a$.

So $aN$ can be used to describe the functional response of a single individual predator. We can plot the number of prey available by the number of prey killed per unit time to visualize the functional response.

```{r}

library(ggplot2)

a <- 0.1
df <- data.frame(N = c(0, 50))

ggplot(data = df, aes(x = N)) +
  stat_function(fun = function(N) a * N) +
  labs(x = "Prey Density", y = "Prey Killed Per Predator")

```

<details>

<summary>Do you think the trend seen in this curve is realistic in nature? Why or why not?</summary>

</details>

$$
\frac{aN}{1 + ahN}
$$ What are each of the terms in this equation?

When are Type II responses seen in nature?

Now plot a Type II functional response:

```{r}


```

In contrast, when do we typically see Type III responses?

$$
\frac{\frac{h}{1}N^{2}}{(\frac{ah}{1})^{2} + N^{2}} = \frac{cN^{2}}{d^{2} + N^{2}}
$$

Where $c = \frac{h}{1}$ is the maximum consumption rate and $d = \frac{ah}{1}$ is the half-saturation constant.

Plot a Type III functional response:

```{r}


```

You will see in your homework how these functional curves are actually edge cases of the same equation, so they are all closely interrelated.

------------------------------------------------------------------------

## Rosenzweig-MacArthur Models

Michael Rosenzweig (a graduate student at the time) and his adviser Robert MacArthur had issues with the applicability of the Lotka-Volterra in natural systems.

Especially the often inaccurate assumptions we have discusses in the past two labs: the Type I response of predators and the density-independence of prey populations. Because of these, Rosenzweig and MacArthur made the following changes to the model:

1.  Predators would be unlikely to follow a Type I response due to:

    -   Predators would become satiated *and/or*

    -   Handling time (`h`) would physically limit the number of prey a predator could consume

<!-- -->

2.  In the absence of predators, prey would become self limiting (density-dependent)

These changes can be integrated into the predator portion of the Lotka-Volterra formula as:

$$
\frac{dP}{dt} \ =\ \frac{faNP}{1\ -\ ahN} \ -\ qP
$$

As for the prey portion, we can add the handling time as follows:

$$
\frac{dN}{dt} \ =\ rN\ -\ \frac{aNP}{1\ +\ ahN}
$$

<details>

<summary>Remember from last lab we modeled density-dependent growth. How would we integrate density-dependence into the Lotka-Volterra model formula?</summary>

We need to integrate the carrying capacity, $K$, into the model. From last we saw that we can represent density-dependent population growth as:

$$
\frac{dN}{dt} \ =\ rN\left( 1\ -\ \frac{N}{K}\right)
$$ Therefore, to integrate into our model we would write:

$$
\frac{dN}{dt} \ =\ rN\left( 1\ -\ \frac{N}{K}\right) \ -\ \frac{aNP}{1\ +\ ahN}
$$

</details>

<details>

<summary>Now put this all together and integrate these additional parameters into the Lotka-Volterra R function we constructed above to write the function `rosenzweig_macarthur`.</summary>

```{r}

```

</details>

<details>

<summary>Then, call this function using the `deSolve::ode` function like we did above with `h <- 0.1` and `K <- 1000` and plot with `ggplot2`.</summary>

```{r}

```

</details>

<details>

<summary>Modify the new parameters of handling time and carrying capacity that we just introduced into our model. What do you observe?</summary>

```{r, echo = FALSE}

library(patchwork)

rosenzweig_macarthur_plot <- function(r = 0.5, a = 0.01, q = 0.2, f = 0.1, h = 0.1, K = 1000, Time = seq(0, 2000, by = 0.1)) {
  
  parms = c(r = r, a = a, q = q, f = f, h = h, K = K)
  
  rosenzweig_macarthur.out <- deSolve::ode(c(N0 = 25, P0 = 5), times = Time, func = rosenzweig_macarthur, parms = parms)
  
  df <- data.frame(Time, rosenzweig_macarthur.out)
  ggplot(df, aes(x = Time)) + 
    geom_line(aes(y = P0), color = "red") +
    geom_line(aes(y = N0), color = "black") +
    labs(y = "Population Size") +
    ggtitle(paste0("f = ", f, "; K = ", K))
}

```

```{r, echo = FALSE}

plots <- list(
  rosenzweig_macarthur_plot(),
  rosenzweig_macarthur_plot(K = 5000),
  rosenzweig_macarthur_plot(K = 20),
  rosenzweig_macarthur_plot(h = 0.001),
  rosenzweig_macarthur_plot(h = 1)
)

patchwork::wrap_plots(plots) + 
  plot_layout(ncol = 3)
```

-   Increasing carrying capacity increases the height of the prey population size each cycle before the predator interaction drives the prey population size downwards again. This destabilizes the system because the prey population is not trending towards a node. This concept is called the "paradox of enrichment."

-   Decreasing the carrying capacity causes the prey population size to drop in a density-independent manner towards stabilization at `K`.

-   Decreasing the handling time causes prey population to stabilize towards a node. The predator is said to be "controlling" the prey population.

-   Increasing the handling time allows the prey population to "escape" the control of the predator. Note that this escape is self-limiting, the prey population size converges at `K`.

</details>

<details>

<summary>What happens if we increase predator efficiency?</summary>

</details>
