10_complex_experiments_end.qmd

---
title: "Tutorial 10: Complex Experiments with the END"
date: last-modified
author: "Danet and Becks, based on originals by Delmas and Griffiths"
format:
    html:
        embed-resources: true
title-block-banner: true
engine: julia
---

::: {.callout-caution}
## Warning

Stability was replaced by Shannon diversity - need to review context as well as if we want that or rather re-integrate stability
:::

The previous tutorial focused on experiments where we manipulated the number of networks and various network parameters.  This is one set of things we can change/vary in an _in silico_ experiment.  The other set of things we can change are features of the model, such as the shape of the functional response (see Tutorial 7), features of the environment such as the carrying capacity, or even empirical relationships that drive trophic structure and interaction strengths, such as the predator-prey mass ratio.

In this tutorial, we are going to implement three experiments.  The first two will be 'simple' in that they vary only two things.  The final example will implement a large experiment changing five features of the model.

You may want to start a new script in the project.  We'll need the following packages (they are already installed... so we just need `using`).

```{julia}
using Random, Plots, Distributions, DataFrames, StatsPlots
using EcologicalNetworksDynamics
```

### Experiment 1: Carrying Capacity and the Predator Prey Mass Ratio

Now we are set for our first experiment.  Lets first establish the parameters we need to make the food web and do the experiment.  We fix `S` at 20 and `C` at 0.15.  We then create vectors of Z and K.  

Z is the predator - prey mass ratio, and defines how much bigger or smaller the predators are from their prey. The data suggest it is between predators are between 10 and 100 times bigger than their prey [see Brose et al 2006](https://doi.org/10.1890/0012-9658(2006)87[2411:CBRINF]2.0.CO;2). This value interacts with setting trophic levels in the model. 

The default setting for the models is 1 - i.e. all species are within the same order of magnitude, predators are not bigger than their prey.  Here, we create a vector of values to explore, from predators being smaller, to them being 10 or 100 x larger as the data suggests.


```{julia}

#Fixed Parameters
S = 20
C = 0.15

# Variable Parameters
Z_levels = [0.1, 1, 10, 100]
K_levels = [0.1, 1, 10, 100]

# run this to get same results as in the document
Random.seed!(123)
```

Now, lets set up the collecting data frame.

```{julia}
df_collect = DataFrame(Z = [], K = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])
```

Now, set up the loop to use these variables and generate outputs. Notice that we use `for z in Z_levels` - this is a clever trick of the looping method, where `z` simply iterates over the values of `Z_levels` without having to specify the index value (e.g. no use of `Z_levels[i]` etc).

The significant BIG thing here is the LogisticGrowth function which allows us to set things like the carrying capacity (K) of the resources.  Here we use it to define custom values of the K paramter for carrying capacity, drawing on the values in the `K_levels` above.  Here, `pg` stands for Producer Growth function, and the paramter set in the food web is K.

Note too our use of `println` and the values of `Z` and `K` to produce an informative _break_ between each combination.

::: {.callout-note icon=false}

Can you guess what increasing K will do to the biomass and richness of the community at equilibrium?  How about Z?  Will higher Z make things more or less stable?
:::



```{julia}
#| output: false

for z in Z_levels
    for k in K_levels

        println(" ***> This is iteration with Z = $z and K = $k\n")

        # Define the food web
        fw = Foodweb(:niche; S = S, C = C)
        # specify the K value of the producer growth function

        B0 = rand(S)
        # specify model to simulate logistic growth as well as BM ratio
        params = default_model(fw, BodyMass(; Z = z), LogisticGrowth(; K = k))
        
        # number of timestamps
        t = 300

        out = simulate(params, B0, t)

        # calculate metrics
        fin_rich = richness(out)
        fin_biomass = total_biomass(out)
        s_div = shannon_diversity(out)

        push!(df_collect, [z, k, fin_rich, fin_biomass, s_div])
    end
end
```

Wonderful.  Now we are in a position to learn about two new plotting methods.  First, let's look at the data frame we've created.

```{julia}
df_collect
```

#### Visualising the experiment

One option here is to plot one of our `Final` Objects as the response variable against the valuse of Z and K.  In R, we'd use ggplot2. Here we'll use `StatsPlots` as we learned about in Tutorial 5.  Can you make this work in the regular `Plots` syntax?

Let's first look at a single plot of stability
```{julia}
#| eval: false

@df df_collect plot(:K, [:FinalStability], group = :Z, 
    ylabel = "Stabilty", 
	xlabel = "Karrying Kapacity",
    seriestype = [:scatter, :line],
    legend = false)
```

Now some new ploting tricks... 3 plots in a layout.

```{julia}
#| eval: false

p1 = @df df_collect plot(:K, [:FinalStability], group = :Z, 
    legend = :bottomright,
    ylabel = "Stabilty", 
	xlabel = "Karrying Kapacity",
    seriestype = [:scatter, :line])

p2 = @df df_collect plot(:K, [:FinalBiomass], group = :Z, 
    legend = :bottomright,
    ylabel = "Biomass", 
	xlabel = "Karrying Kapacity",
    seriestype = [:scatter, :line])
    
p3 = @df df_collect plot(:K, [:FinalRichness], group = :Z, 
    legend = :bottomright,
    ylabel = "Richness", 
	xlabel = "Karrying Kapacity",
    seriestype = [:scatter, :line])

# create a layout of 3 graphs stacked on top of each other.
plot(p1, p2, p3, layout=(3,1), legend = false)
```

### Interpretation!

#### Challenge - can you get the number of extinctions into the data frame?

### Experiment 2: The Functional Response

The functional response is the relationship between how much a consumer eats and the 'density' of the prey items.  If you can recall from your ecology courses/classes/modules, there are three classic shapes: The Type I, Type II and Type III.

A predator feeding with a Type I delivers to the prey a 'constant mortality rate' (the slope of the Type I line).  This means that the effect of predation is density _independent_ because prey mortality rate does not vary by prey density.  Remember, density dependence (negative feedback that stabilises communities) is defined by survival decreasing with increasing density, or in this case, mortality rates _increasing_ with increasing density.

A predator feeding with the Type II delivers an _inverse density dependent_ mortality rate.  The slope of the Type II line actually goes down as density of the prey goes up meaning that mortality rates for the prey, caused by the predator, are going down with prey density.  This means that the effect of predation is _inverse density dependent_ in the Type II.  This is **destabilising**.

Finally, a predator feeding via a Type III can deliver a _density dependent_ mortality rate to the prey, but only at low prey densities.  This is an S shaped curve.  Below the inflection point, the slope is actually getting steeper.  This means that as prey density increases up to the inflection, their mortality rate from predation increases (survival goes down with density going up). This is the hallmark of density dependence and can **stabilise** consumer-resource interactions.

::: {.callout-tip icon=false}

Remember that the logistic growth equation, with a carying capacity specified, is also a source of _density dependent negative feedback_
:::

::: {.callout-tip icon=false}

The Type II is the MOST common.  Type I is rare and even non-existent because it suggests there are no limits to how much a consumer can eat.  Type III is also rare, but it is at least plausible and interesting.
:::

```{julia}

f_t1(n) = 0.5*n
f_t2(n) = 0.5*n/(0.2+0.01*n)
f_t3(n) = 0.5*n^2/(10 + 0.01*n^2)

plot(f_t1, 0, 100, label = "Type I")
plot!(f_t2, 0, 100, label = "Type II")
plot!(f_t3, 0, 100, label = "Type III")
```

#### How does the BEFW make a functional response?

There are two formulations of the functional response.  One of them is called the _Bioenergetic_ response and the other is called the _Classic_.  In both cases, we ignore the Type I.

The Bioenergetic functional response is deeply phenomenological in that the parameters that make the shapes move between Type II and III have no deliberate biological interpretation.  They function is defined by a 1/2 saturation point, an asymptote (which is nominally a maxiumum feeding rate) and an exponent, which is called the _hill exponent_.  The value of the exponent moves the model from Type II (h = 1) to Type III (h = 2).  The other variables define the overall shape.

The Classic functional less phenomenological in that the response is defined more by 'traits': the attack rate of a consumer on a prey and the handling time of that prey.  But it also moves between the Type II and Type III shape based on an exponent.

#### Creating Type II vs. Type III with the Bioenergetic response

Let's look at using the Bioenergetic functional response, and see here how we can vary the shape between Type II and Type III. We can do this by modifying the *hill_exponent* after we have specified the model (*i.e.,* after the `default_model` call). We will look at how Richness, Biomass and Shannon Diversity are affected by the hill exponent.

```{julia}
#| output: false

Random.seed!(12352)

# fixed parameters
S = 20
C = 0.15

# set the hill exponent to move from Type II to Type III)
h_levels = [1.0, 1.1, 1.25, 2.0]

# set collecting data frame 
# we will look at how Richness, Biomass and Shannon Diversity are affected by the hill exponent
df_collect_h = DataFrame(h = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])

# create look across values of h
for h in h_levels 
    println("***> This is iteration with h = $h\n")
    
    # make the network
    # Note that we specify the Environment, but there is no K or T set (using defaults)
    # Note the new BioenergeticResponse function
    fw_h = Foodweb(:niche; S = S, C = C)
    
    # set body sizes and parameters 
    B0 = rand(S)
    params = default_model(fw_h)

    # here we now update the exponent of the hill function
    params.hill_exponent = h

    # specify number of time steps
    t = 300

    # simulate
    sim_niche = simulate(params, B0, t)

    # collect data 
    fin_rich = richness(sim_niche)
    fin_bio = total_biomass(sim_niche)
    s_div = shannon_diversity(sim_niche)

    push!(df_collect_h, [h, fin_rich, fin_bio, s_div])
end

df_collect_h
```

Now, we can visualise these data

```{julia}
#| output: false

# Visualize the results
p1_h = @df df_collect_h plot(:h, [:FinalRichness],
    legend = :bottomright,
    ylabel = "Richness",
    xlabel = "Functional response",
    seriestype = [:scatter, :line])

p2_h = @df df_collect_h plot(:h, [:FinalBiomass],
    legend = :bottomright,
    ylabel = "Biomass",
    xlabel = "Functional response",
    seriestype = [:scatter, :line])

p3_h = @df df_collect_h plot(:h, [:ShannonDiversity],
    legend = :bottomright,
    ylabel = "Shannon Diversity",
    xlabel = "Functional response",
    seriestype = [:scatter, :line])

plot(p1_h, p2_h, p3_h, layout=(3,1), legend = false, size = (1000, 1000))
```

#### INTERPRETATION?

What can you see happening as we move away from the destabilising Type II functional response?

Can you modify this code to explore what happens at different values of K? You'll need to modify this section, and the collection data frame.

```{julia}
#| eval: false

   # make the network
    fw_h = Foodweb(:niche; S = S, C = C)

    # set body sizes and parameters 
    B0 = rand(S)
    params = default_model(fw_h, LogisticGrowth(; K = k))

    # update the exponent of the hill function
    params.hill_exponent = h
```

### Experiment 3: What is Z

One of the central features of the link between the Bioenergetic Food Web model and the structure of a foodweb created by models like the Niche Model is the organisation of trophic levels.  At the heart of this is a _data driven_ assumption about the ratio of predator size to prey size.  This is called the _Predator Prey Mass Ratio_, or `PPMR` for short.  

In 2006, Uli Brose and team collated hundreds of data [to reveal that](https://esajournals.onlinelibrary.wiley.com/doi/10.1890/0012-9658%282006%2987%5B2411%3ACBRINF%5D2.0.CO%3B2), on average, predators were between 10 and 100x bigger than their prey.

In our modelling framework, we use this ratio to help organise species into trophic levels.  This is done by organising the bodymass vector, and via a parameter called `Z`.  The body mass of consumers is a function of their mean trophic level (T), and it increases with trophic level when Z ≥ 1 and decreases when Z ≤ 1 via this relationship (see Delmas et al 2017 and Williams et al 2007):

$M_C = Z^(T-1)$

[Brose et al 2006](https://onlinelibrary.wiley.com/doi/10.1111/j.1461-0248.2006.00978.x) explored the impact of the _PPMR_ on stability and dynamics as part of their wider exploration of scaling and allometry in the bioenergetic model.  Here we show you how to manipulate `Z` and it's effect on stability.  `Z` is specified in the call to FoodWeb as the allocation of species with specific sizes is central to the trophic structure of the model.  This argument is interfaced with the bodysize vector in `model_parameters()`

```{julia}
#| output: false

Random.seed!(12352)

# fixed parameters
S = 20
C = 0.15

# set the PPRM
z_levels= [0.1, 1, 10, 100]

# set collecting data frame 
# we will look at how Richness, Biomass and Shannon Diversity are affected by the hill exponent
df_collect_z = DataFrame(z = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])

# create look across values of h
for z in z_levels 
    println("***> This is iteration with z = $z\n")
    
    # make the network
    # Note that we specify the Environment, but there is no K or T set (using defaults)
    # Note Z is specified when building the FoodWeb() network
    fw_z = Foodweb(:niche; S = S, C = C)
    
    # set body sizes and parameters 
    B0 = rand(S)
    params = default_model(fw_z, BodyMass(; Z = z))

    # specify number of time steps
    t = 300

    # simulate
    out_z = simulate(params, B0, t)

    # collect data 
    fin_rich = richness(out_z)
    fin_bio = total_biomass(out_z)
    s_div = shannon_diversity(out_z)

    push!(df_collect_z, [z, fin_rich, fin_bio, s_div])
end

df_collect_z
```

As with the variation in `h`, we can create a set of figures too!  Perhaps it's worth your time to consult [Brose et al 2006](https://onlinelibrary.wiley.com/doi/10.1111/j.1461-0248.2006.00978.x) and [Reger et al 2017](https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.12713) to make sure you understand how Z works and particularly how stability is expected to vary with Z!  One of the most important things to understanding is why the stability metric is negative and what values of stability close, or far away, from zero mean.

```{julia}
#| output: false

# Visualize the results
p1_z = @df df_collect_z plot(:z, [:FinalRichness],
    legend = :bottomright,
    ylabel = "Richness",
    xlabel = "Z value (PPMR)",
    seriestype = [:scatter, :line])

p2_z = @df df_collect_z plot(:z, [:FinalBiomass],
    legend = :bottomright,
    ylabel = "Biomass",
    xlabel = "Z value (PPMR)",
    seriestype = [:scatter, :line])

p3_z = @df df_collect_z plot(:z, [:ShannonDiversity],
    legend = :bottomright,
    ylabel = "Shannon Diversity",
    xlabel = "Z value (PPMR)",
    seriestype = [:scatter, :line])

plot(p1_z, p2_z, p3_z, layout=(3,1), legend = false, size = (1000, 1000))
```

#### Challenge

Could you modify this to ask about the interaction between Z and K?  This would be asking the question of whether the effect of PPMR on stability varies by system productivity.  Or you could ask about the interaction between the functional response, which we know also has a direct effect on stability by the assumption we make of a Type II or Type III shape, and the value of Z, which we also know impacts stability from Brose et al's work.

### Experiment 4: Manipulating Competition among producers

Our final experiment for this section involves manipulating how the basal producers compete with each other.  The default paramterisation of the model has each producer growing via the logistic equation and competing with itself via density dependence.  There is only intraspecific competition, no interspecific competition.

We can modify this assumption by invoking another function called `ProducerCompetition`.  This function acts like `LogisticGrowth` that we use to set `K` and `BioenergeticResponse` that we used to modify the functional response between Type II and Type III.

The theory to recall is that coexistence among species is mediated by the balance between intraspecific and interspecific competition.  When intraspecific competition is greater than interspecific competition, there is coexistence.  However, when interspecific competition is greater than intraspecific competition, there will be competitive exclusion and no coexistence.  

We call the competition parameters $\alpha$. $\alpha~ii$ defines intraspecific competition and $\alpha~ij$ defines interspecific competition.  The $\alpha~ij$ defines how the species $j$ reduces the carrying capacity (equilibrium) of species $i$.  

What we can do is set $\alpha~ii = 1$ and then vary $\alpha~ij$ from $<1$ to $>1$.  We can expect that there will a dramatic change in biomass and species richness as we move from $alpha~ii > alpha~ij$ to $alpha~ii < alpha~ij$.


```{julia}
#| output: false

S = 20 # define the number of species
C = 0.2 # define the connectance (complexity) of the network
Z = 100 # Predator Prey Mass Ratio
t = 300 # time steps

# here we set the 
interspecific_vals = 0.8:0.05:1.2 # a set of (9) values between 0.8 and 1.2 in steps of 0.05
# collect(0.8:0.05:1.2) # see them if you want to


# set collecting data frame 
# we will look at how Richness, Biomass and Shannon Diversity are affected by interspecific competition
df_collect_comp = DataFrame(InterComp = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])

for alpha_ij in interspecific_vals
    println("***> This is iteration with alpha_ij = $alpha_ij\n")
    
    # this will make always the same network and body mass
    Random.seed!(123)
    foodweb = Foodweb(:niche; S = S, C = C)

    # enhance detail in the network
    #We fix intraspecific = 1 (diag) and vary interspecific (off)
    p = ProducersCompetition(foodweb; diag = 1.0, off = alpha_ij)
    #Specify K and competition in LogisticGrowth
    r = LogisticGrowth(K = 10, producers_competition = p)

    #Set the hill exponent to 1 (type II Functional Response)
    b = BioenergeticResponse(h = 1)

    #Define parameters with extras
    params_comp = default_model(foodweb, BodyMass(; Z = Z), b, r)

    # set initial biomasses
    B0 = rand(S)

    # simulate
    # note verbose = false ; remove this to see extinction detail for each sim
    out_comp = simulate(params_comp, B0, t, verbose = false)

    # generate stats and add to collector
    # collect data 
    fin_rich = richness(out_comp)
    fin_bio = total_biomass(out_comp)
    s_div = shannon_diversity(out_comp)

    push!(df_collect_comp, [alpha_ij, fin_rich, fin_bio, s_div])
end

df_collect_comp
```

Let's review the assumptions above.  We've set the `Predator Prey Mass Ratio` to 100.  We've set carrying capacity `K` to 10.  We've set the functional response `h` value to 1, so it's a Type II functional response.  Finally, we've set a range of interspecific competition to be 0.8 to 1.2 around the fixed intraspecific effect of 1.


```{julia}
#| output: false

p1 = @df df_collect_comp plot(:InterComp, [:FinalRichness],
    ylabel = "Richness",
    xlabel = "alpha_ij",
    seriestype = [:scatter, :line])
p2 = @df df_collect_comp plot(:InterComp, [:FinalBiomass],
    ylabel = "Biomass",
    xlabel = "alpha_ij",
    seriestype = [:scatter, :line])
p3 = @df df_collect_comp plot(:InterComp, [:ShannonDiversity],
    ylabel = "Shannon Diversity",
    xlabel = "alpha_ij",
    seriestype = [:scatter, :line])

plot(p1, p2, p3, layout = (3,1), legend = false)
```

#### Challenge - Competition

Perhaps consider expanding the code above to assess one of these?

Is this pattern sensitive to species richness?
Is this pattern sensitive to the functional response?
Is this pattern sensitive to the PPMR?
Is this pattern sensitive to values of K?

## Experiment 5: Multiple networks (replicates)

To Do: run S (3 values), C (3 values) and h (3 values) where there are 5 replicate networks per combination.  Note we need 45 networks...