Tutorial 10: Complex Experiments with the END

Author

Danet and Becks, based on originals by Delmas and Griffiths

Published

November 14, 2024

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

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

#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)
TaskLocalRNG()

Now, lets set up the collecting data frame.

df_collect = DataFrame(Z = [], K = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])
0×5 DataFrame
Row Z K FinalRichness FinalBiomass ShannonDiversity
Any Any Any Any Any

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.

Note

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?

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
 ***> This is iteration with Z = 0.1 and K = 0.1

 ***> This is iteration with Z = 0.1 and K = 1.0

 ***> This is iteration with Z = 0.1 and K = 10.0

 ***> This is iteration with Z = 0.1 and K = 100.0

 ***> This is iteration with Z = 1.0 and K = 0.1

 ***> This is iteration with Z = 1.0 and K = 1.0

 ***> This is iteration with Z = 1.0 and K = 10.0

 ***> This is iteration with Z = 1.0 and K = 100.0

 ***> This is iteration with Z = 10.0 and K = 0.1

 ***> This is iteration with Z = 10.0 and K = 1.0

 ***> This is iteration with Z = 10.0 and K = 10.0

 ***> This is iteration with Z = 10.0 and K = 100.0

 ***> This is iteration with Z = 100.0 and K = 0.1

 ***> This is iteration with Z = 100.0 and K = 1.0

 ***> This is iteration with Z = 100.0 and K = 10.0

 ***> This is iteration with Z = 100.0 and K = 100.0

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.

df_collect
16×5 DataFrame
Row Z K FinalRichness FinalBiomass ShannonDiversity
Any Any Any Any Any
1 0.1 0.1 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 5, 5, 5, 5, 5, 5, 5, 5, 5, 5] [8.15104, 7.35599, 6.86766, 6.20173, 5.62064, 4.93607, 4.31825, 3.63121, 3.00669, 2.39632 … 0.300001, 0.300001, 0.300001, 0.300001, 0.300001, 0.300001, 0.300001, 0.300001, 0.300001, 0.3] [16.6229, 16.7342, 16.5646, 16.1151, 15.6393, 15.1269, 14.7861, 14.5788, 14.548, 14.6324 … 3.00016, 3.00016, 3.00016, 3.00016, 3.00016, 3.00016, 3.00016, 3.00016, 3.00016, 3.00005]
2 0.1 1.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 11, 11, 11, 11, 11, 11, 11, 11, 11, 11] [12.4772, 11.4754, 10.7631, 9.95902, 9.14621, 8.41858, 7.61812, 6.84285, 6.08118, 5.33149 … 1.76616, 1.76616, 1.76616, 1.76616, 1.76616, 1.76616, 1.76616, 1.76616, 1.76616, 1.76616] [17.742, 17.8484, 17.6835, 17.2513, 16.4331, 15.4516, 14.5648, 14.1409, 13.9903, 13.8445 … 7.3824, 7.3824, 7.3824, 7.3824, 7.3824, 7.3824, 7.3824, 7.3824, 7.3824, 7.38235]
3 0.1 10.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 13, 13, 13, 13, 13, 13, 13, 13, 13, 13] [9.25026, 8.51092, 8.00559, 7.44123, 6.91671, 6.27631, 5.65377, 4.94222, 4.26539, 3.54824 … 2.25548, 2.2559, 2.25643, 2.25698, 2.25747, 2.25783, 2.25806, 2.25819, 2.25825, 2.25825] [14.6862, 14.6468, 14.2896, 13.6899, 12.9845, 12.0663, 11.3051, 10.7473, 10.5693, 10.7805 … 9.57129, 9.59447, 9.62346, 9.6528, 9.67659, 9.69307, 9.70369, 9.71014, 9.71404, 9.71427]
4 0.1 100.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [11.7083, 11.0424, 10.5769, 9.93713, 9.4234, 8.80592, 8.23005, 7.54626, 6.97782, 6.26384 … 3.48233, 3.48295, 3.48399, 3.48495, 3.48596, 3.48667, 3.48716, 3.48743, 3.48757, 3.48759] [18.2061, 17.7393, 17.1158, 15.8476, 14.5755, 13.1387, 12.2288, 11.7098, 11.6218, 11.8449 … 14.6644, 14.6766, 14.6955, 14.7124, 14.7303, 14.7429, 14.7523, 14.7579, 14.7612, 14.762]
5 1.0 0.1 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 14, 11, 11, 7, 7, 7, 7, 6, 6, 6] [9.94865, 9.31694, 8.91556, 8.37917, 7.98085, 7.48504, 7.05457, 6.53892, 6.03007, 5.45877 … 0.599995, 0.599995, 0.600001, 0.600001, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6] [16.9057, 16.7928, 16.5411, 15.9674, 15.4178, 14.6761, 13.9969, 13.1036, 12.1136, 10.9472 … 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0]
6 1.0 1.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 17, 17, 16, 16, 16, 16, 15, 15, 14, 14] [8.18822, 7.87201, 7.51513, 7.09849, 6.67715, 6.24736, 5.77386, 5.30732, 4.78624, 4.26858 … 2.1041, 2.1036, 2.1036, 2.10339, 2.10338, 2.10338, 2.10338, 2.1034, 2.1034, 2.10341] [17.8751, 17.9835, 17.9691, 17.7811, 17.4596, 17.0672, 16.6213, 16.2204, 15.8747, 15.676 … 9.70326, 9.70168, 9.70168, 9.70046, 9.7004, 9.7004, 9.7004, 9.70055, 9.70055, 9.70071]
7 1.0 10.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 18, 18, 18, 18, 17, 17, 17, 17, 17, 17] [11.9214, 11.5912, 11.2096, 10.7345, 10.2291, 9.65051, 9.08087, 8.50293, 7.95285, 7.4397 … 12.5097, 12.5097, 12.5097, 12.5098, 12.5098, 12.5098, 12.5098, 12.5098, 12.5098, 12.5098] [17.7631, 17.8918, 17.8487, 17.5309, 16.9389, 15.9573, 14.662, 13.1789, 11.9185, 11.0675 … 4.62639, 4.61727, 4.60955, 4.60285, 4.60285, 4.59699, 4.59175, 4.58698, 4.58255, 4.57863]
8 1.0 100.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [9.81854, 9.78941, 9.7508, 9.70287, 9.63915, 9.53889, 9.36204, 9.09799, 8.78533, 8.50462 … 7.83776, 7.83881, 7.83874, 7.83824, 7.83785, 7.8378, 7.83791, 7.83792, 7.83791, 7.83792] [15.8542, 15.9962, 16.1459, 16.2131, 16.1391, 15.9368, 15.6428, 15.2171, 14.6976, 14.3231 … 14.8658, 14.8748, 14.8756, 14.873, 14.8704, 14.8698, 14.8703, 14.8705, 14.8704, 14.8705]
9 10.0 0.1 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 13, 13, 12, 12, 10, 10, 10, 9, 9, 8] [12.3589, 11.1366, 10.528, 9.73548, 9.28095, 8.70753, 8.38507, 7.92468, 7.6241, 7.16456 … 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7] [17.9611, 18.047, 17.9192, 17.5077, 17.121, 16.4637, 16.0159, 15.3056, 14.8254, 14.13 … 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0]
10 10.0 1.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 19, 19, 19, 19, 19, 19, 19, 19] [8.42215, 8.20465, 7.99326, 7.76301, 7.53365, 7.30894, 7.04805, 6.82182, 6.50908, 6.21272 … 2.86069, 2.86431, 2.86431, 2.86744, 2.87126, 2.87422, 2.87704, 2.87913, 2.88089, 2.88177] [15.0882, 15.2155, 15.2234, 15.0786, 14.8185, 14.5613, 14.3405, 14.2066, 14.0631, 13.9396 … 13.2383, 13.2133, 13.2133, 13.1887, 13.1531, 13.1189, 13.0773, 13.0366, 12.9906, 12.9603]
11 10.0 10.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 17, 17, 17, 17, 17, 17, 16, 16, 16, 16] [8.76005, 8.6219, 8.45014, 8.27291, 8.03749, 7.83254, 7.54527, 7.26893, 6.89872, 6.5367 … 4.352, 4.37513, 4.3948, 4.41157, 4.4259, 4.43803, 4.43803, 4.44812, 4.45638, 4.45893] [16.3521, 16.3257, 16.2454, 16.1206, 15.8993, 15.6612, 15.275, 14.8714, 14.3179, 13.7872 … 5.58795, 5.46404, 5.36072, 5.2744, 5.20203, 5.14192, 5.14192, 5.09276, 5.0532, 5.04108]
12 10.0 100.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [8.25211, 8.29101, 8.4592, 8.87757, 9.76632, 11.7568, 16.001, 25.5672, 46.8751, 88.1647 … 243.87, 243.868, 243.863, 243.857, 243.852, 243.848, 243.845, 243.843, 243.841, 243.84] [15.4288, 15.4367, 15.3683, 14.9952, 14.0388, 12.1209, 9.40314, 6.56688, 4.55697, 3.55452 … 6.58682, 6.58612, 6.58472, 6.58295, 6.58118, 6.57958, 6.57821, 6.57708, 6.57618, 6.57556]
13 100.0 0.1 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 19, 19, 18, 18, 18, 18, 18, 18, 18, 18] [11.5095, 10.9511, 10.7168, 10.3916, 10.2152, 9.9237, 9.78424, 9.46529, 9.30356, 8.92472 … 0.472889, 0.463688, 0.463688, 0.455691, 0.448684, 0.442486, 0.436958, 0.431991, 0.427513, 0.425609] [17.231, 17.3743, 17.335, 17.1761, 17.0496, 16.8121, 16.6962, 16.4432, 16.324, 16.0738 … 5.95521, 5.74108, 5.74108, 5.5499, 5.37858, 5.22417, 5.08415, 4.95644, 4.83969, 4.78959]
14 100.0 1.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [8.94845, 8.84267, 8.70565, 8.58269, 8.46112, 8.33459, 8.2076, 8.07545, 7.93311, 7.80568 … 8.84776, 8.83373, 8.82279, 8.81646, 8.81702, 8.82562, 8.83778, 8.84947, 8.8573, 8.85948] [16.8594, 16.831, 16.7339, 16.6093, 16.4917, 16.4208, 16.4238, 16.5015, 16.6608, 16.8785 … 11.6325, 11.5115, 11.4211, 11.3514, 11.3021, 11.2653, 11.2375, 11.2081, 11.1754, 11.1534]
15 100.0 10.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [9.0416, 9.1737, 9.5697, 10.2502, 11.3574, 13.2685, 16.0394, 20.0685, 25.3442, 31.3042 … 79.4917, 79.5051, 79.498, 79.4849, 79.4774, 79.4781, 79.4801, 79.4797, 79.4785, 79.4782] [15.9684, 16.1206, 16.4693, 16.7484, 16.6086, 15.5647, 13.7709, 11.8335, 10.3515, 9.58567 … 11.1015, 11.0926, 11.0865, 11.0856, 11.0869, 11.0872, 11.0863, 11.0855, 11.0854, 11.0854]
16 100.0 100.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [11.5292, 11.6221, 11.8303, 12.2001, 12.9034, 14.2178, 16.6995, 21.3803, 29.102, 40.6004 … 466.619, 464.416, 463.06, 460.625, 457.236, 453.504, 450.831, 452.672, 465.621, 471.999] [18.3089, 18.2607, 18.15, 17.9446, 17.542, 16.8141, 15.6206, 13.9258, 12.0175, 9.86971 … 7.06769, 7.20399, 7.09213, 6.94422, 6.80881, 6.70295, 6.69553, 6.94315, 7.66919, 7.96389]

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

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

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.

Tip

Remember that the logistic growth equation, with a carying capacity specified, is also a source of density dependent negative feedback

Tip

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.

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.

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 Stability 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
***> This is iteration with h = 1.0

***> This is iteration with h = 1.1

***> This is iteration with h = 1.25

***> This is iteration with h = 2.0
4×4 DataFrame
Row h FinalRichness FinalBiomass ShannonDiversity
Any Any Any Any
1 1.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 16, 16, 15, 15, 15, 15, 15, 15, 15, 15] [14.0154, 13.7814, 13.4658, 13.1783, 12.8834, 12.6007, 12.3015, 11.9503, 11.7107, 11.4361 … 6.42979, 6.44208, 6.44208, 6.44436, 6.43148, 6.44421, 6.44068, 6.43311, 6.44545, 6.44344] [18.7295, 18.6847, 18.5806, 18.4524, 18.2893, 18.0832, 17.7706, 17.2078, 16.6571, 15.8173 … 7.82943, 7.8319, 7.8319, 7.84614, 7.82841, 7.83428, 7.84109, 7.82605, 7.83731, 7.84078]
2 1.1 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 17, 17, 17, 17, 17, 17, 16, 16, 16, 16] [11.3827, 11.1345, 10.7624, 10.3913, 9.99559, 9.59857, 9.16844, 8.74251, 8.37833, 8.02815 … 5.77289, 5.7729, 5.77293, 5.77292, 5.7729, 5.77291, 5.77291, 5.77291, 5.77291, 5.77291] [17.2879, 17.2245, 17.0354, 16.7364, 16.3009, 15.7385, 14.957, 13.9498, 12.8428, 11.4674 … 8.05783, 8.05797, 8.05789, 8.05783, 8.0578, 8.05778, 8.05778, 8.05775, 8.05773, 8.05772]
3 1.25 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [7.95051, 7.71008, 7.36973, 7.05915, 6.81253, 6.60547, 6.33242, 6.20856, 5.90247, 5.79213 … 5.73092, 5.7267, 5.72288, 5.71981, 5.71723, 5.71558, 5.71443, 5.71385, 5.71357, 5.71354] [14.8364, 15.0814, 15.2094, 15.0455, 14.7772, 14.5307, 14.2246, 14.1041, 13.874, 13.816 … 12.4849, 12.445, 12.4101, 12.3868, 12.3692, 12.3593, 12.3531, 12.3504, 12.3493, 12.3492]
4 2.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [8.38822, 8.15315, 7.90496, 7.74148, 7.54947, 7.42629, 7.21904, 7.06338, 6.82622, 6.62508 … 5.63687, 5.63973, 5.64247, 5.64441, 5.64568, 5.64631, 5.64656, 5.64657, 5.64648, 5.64641] [13.7108, 13.8974, 13.9356, 13.9101, 13.8971, 13.9043, 13.9225, 13.9274, 13.9095, 13.8705 … 13.7105, 13.7429, 13.7652, 13.7748, 13.7762, 13.7723, 13.7654, 13.7572, 13.7488, 13.7434]

Now, we can visualise these data

# Visualize the results
p1_h = @df df_collect_h plot(:h, [:FinalStability],
    legend = :bottomright,
    ylabel = "Stability",
    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, [:FinalRichness],
    legend = :bottomright,
    ylabel = "Richness",
    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.

   # 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, 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 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()

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 Stability 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
***> This is iteration with z = 0.1

***> This is iteration with z = 1.0

***> This is iteration with z = 10.0

***> This is iteration with z = 100.0
4×4 DataFrame
Row z FinalRichness FinalBiomass ShannonDiversity
Any Any Any Any
1 0.1 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 14, 14, 14, 14, 14, 13, 13, 12, 12, 12] [14.0154, 13.1242, 12.5054, 11.7867, 11.1144, 10.4411, 9.65219, 8.81407, 7.85734, 6.85837 … 2.54144, 2.54144, 2.54144, 2.54144, 2.54144, 2.54144, 2.54144, 2.54144, 2.54144, 2.54144] [18.7295, 18.6017, 18.0408, 17.0915, 16.2886, 15.7357, 15.3028, 14.9264, 14.5273, 14.1848 … 9.01235, 9.01233, 9.01232, 9.01231, 9.01231, 9.01231, 9.01231, 9.01231, 9.0123, 9.0123]
2 1.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 19, 19, 17, 17, 16, 16, 16, 16, 16, 16] [11.3827, 11.0237, 10.6428, 10.1558, 9.68492, 9.08034, 8.44996, 7.73821, 7.12812, 6.50964 … 4.04905, 4.04904, 4.04904, 4.04904, 4.04904, 4.04903, 4.04903, 4.04902, 4.04902, 4.04902] [17.2879, 17.2885, 17.2049, 16.9595, 16.5593, 15.7827, 14.6068, 12.8519, 11.257, 9.97863 … 12.411, 12.411, 12.411, 12.411, 12.411, 12.4109, 12.4109, 12.4108, 12.4108, 12.4108]
3 10.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [7.95051, 7.71606, 7.43686, 7.21154, 7.03513, 6.86424, 6.73195, 6.58726, 6.478, 6.39076 … 6.19265, 6.19575, 6.19649, 6.19564, 6.19442, 6.19361, 6.19319, 6.19306, 6.19304, 6.19304] [14.8364, 15.1738, 15.4241, 15.5023, 15.5318, 15.5827, 15.6634, 15.816, 15.987, 16.1622 … 17.1402, 17.1371, 17.1269, 17.1157, 17.1054, 17.0988, 17.0953, 17.0943, 17.0942, 17.0942]
4 100.0 [20, 20, 20, 20, 20, 20, 20, 20, 20, 20 … 20, 20, 20, 20, 20, 20, 20, 20, 20, 20] [8.38822, 8.25788, 8.10381, 7.99055, 7.9029, 7.81458, 7.76517, 7.70146, 7.67839, 7.69968 … 10.4944, 10.5011, 10.5112, 10.5238, 10.5387, 10.553, 10.5648, 10.5725, 10.5769, 10.5779] [13.7108, 13.8452, 13.9135, 13.9026, 13.8919, 13.9322, 13.9969, 14.1639, 14.3154, 14.5913 … 12.0742, 12.1247, 12.1854, 12.2463, 12.3038, 12.3476, 12.3764, 12.3919, 12.3991, 12.4006]

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 and Reger et al 2017 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.

# Visualize the results
p1_z = @df df_collect_z plot(:z, [:FinalStability],
    legend = :bottomright,
    ylabel = "Stability",
    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, [:FinalRichness],
    legend = :bottomright,
    ylabel = "Richness",
    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 Environment 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 compeitive 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\).

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

# 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 Stability are affected by the hill exponent
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, Z = Z)

    # enhance detail in the network
    # specify K
    LogisticGrowth(foodweb; K = 10)
    # we fix intraspecific = 1 and vary interspecific
    ProducerCompetiton(foodweb; αii = 1.0, αij = 1)
    # set the hill exponent to 1 (type II Functional Response)
    BioenergeticResponse(foodweb, h = 1)

    # define parameters with extras
    params_comp = default_model(foodweb, BodyMass(; Z = z))

    # set bodymass
    B0 = rand(S)

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

    # generate stats and add to collector
    # collect data 
    fin_rich = richness(out_comp)
    fin_bio = biomass(out_comp).total
    stab = community_cv(out_comp)

    push!(df_collect_comp, [alpha_ij, fin_rich, fin_bio, stab])
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.

p1 = @df df_collect_comp plot(:InterComp, [:FinalRichness],
    ylabel = "Richness",
    xlabel = "alpha_ij")
p2 = @df df_collect_comp plot(:InterComp, [:FinalBiomass],
    ylabel = "Biomass",
    xlabel = "alpha_ij")
p3 = @df df_collect_comp plot(:InterComp, [:FinalStability],
    ylabel = "Stability",
    xlabel = "alpha_ij")

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 specie 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…