diff --git a/10_complex_experiments_end.qmd b/10_complex_experiments_end.qmd index f77a6b2..a5552e5 100644 --- a/10_complex_experiments_end.qmd +++ b/10_complex_experiments_end.qmd @@ -378,9 +378,9 @@ Could you modify this to ask about the interaction between Z and K? This would 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. +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 compeitive exclusion and no coexistence. +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$. @@ -393,6 +393,7 @@ What we can do is set $\alpha~ii = 1$ and then vary $\alpha~ij$ from $<1$ to $>1 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 @@ -400,7 +401,7 @@ interspecific_vals = 0.8:0.05:1.2 # a set of (9) values between 0.8 and 1.2 in s # set collecting data frame -# we will look at how Richness, Biomass and Stability are affected by the hill exponent +# 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 @@ -411,30 +412,31 @@ for alpha_ij in interspecific_vals 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) + #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)) + #Define parameters with extras + params_comp = default_model(foodweb, BodyMass(; Z = Z), b, r) - # set bodymass + # 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, verbose = false) + out_comp = simulate(params_comp, B0, t, 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) + fin_bio = total_biomass(out_comp) + s_div = shannon_diversity(out_comp) - push!(df_collect_comp, [alpha_ij, fin_rich, fin_bio, stab]) + push!(df_collect_comp, [alpha_ij, fin_rich, fin_bio, s_div]) end df_collect_comp @@ -451,8 +453,8 @@ p1 = @df df_collect_comp plot(:InterComp, [:FinalRichness], p2 = @df df_collect_comp plot(:InterComp, [:FinalBiomass], ylabel = "Biomass", xlabel = "alpha_ij") -p3 = @df df_collect_comp plot(:InterComp, [:FinalStability], - ylabel = "Stability", +p3 = @df df_collect_comp plot(:InterComp, [:ShannonDiversity], + ylabel = "Shannon Diversity", xlabel = "alpha_ij") plot(p1, p2, p3, layout = (3,1), legend = false) @@ -462,11 +464,11 @@ plot(p1, p2, p3, layout = (3,1), legend = false) Perhaps consider expanding the code above to assess one of these? -Is this pattern sensitive to specie richness? +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... \ No newline at end of file +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... diff --git a/_freeze/10_complex_experiments_end/execute-results/html.json b/_freeze/10_complex_experiments_end/execute-results/html.json index e24f799..2e3d1d3 100644 --- a/_freeze/10_complex_experiments_end/execute-results/html.json +++ b/_freeze/10_complex_experiments_end/execute-results/html.json @@ -1,10 +1,10 @@ { - "hash": "93f6b735a5eb08313674e48cf128e6df", + "hash": "d06996d6456bea5a1dc5ee061072c5a0", "result": { - "engine": "jupyter", - "markdown": "---\ntitle: \"Tutorial 10: Complex Experiments with the END\"\ndate: last-modified\nauthor: \"Danet and Becks, based on originals by Delmas and Griffiths\"\nformat:\n html:\n embed-resources: true\ntitle-block-banner: true\n---\n\n\n\n\n\n::: {.callout-caution}\n## Warning\n\nStability was replaced by Shannon diversity - need to review context as well as if we want that or rather re-integrate stability\n:::\n\nThe 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.\n\nIn 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.\n\nYou 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`).\n\n::: {#6c9caa58 .cell execution_count=1}\n``` {.julia .cell-code}\nusing Random, Plots, Distributions, DataFrames, StatsPlots\nusing EcologicalNetworksDynamics\n```\n:::\n\n\n### Experiment 1: Carrying Capacity and the Predator Prey Mass Ratio\n\nNow 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. \n\nZ 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. \n\nThe 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.\n\n::: {#7aeb281d .cell execution_count=2}\n``` {.julia .cell-code}\n#Fixed Parameters\nS = 20\nC = 0.15\n\n# Variable Parameters\nZ_levels = [0.1, 1, 10, 100]\nK_levels = [0.1, 1, 10, 100]\n\n# run this to get same results as in the document\nRandom.seed!(123)\n```\n\n::: {.cell-output .cell-output-display execution_count=3}\n```\nTaskLocalRNG()\n```\n:::\n:::\n\n\nNow, lets set up the collecting data frame.\n\n::: {#37e1e107 .cell execution_count=3}\n``` {.julia .cell-code}\ndf_collect = DataFrame(Z = [], K = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])\n```\n\n::: {.cell-output .cell-output-display execution_count=4}\n```{=html}\n
0×5 DataFrame
RowZKFinalRichnessFinalBiomassShannonDiversity
AnyAnyAnyAnyAny
\n```\n:::\n:::\n\n\nNow, 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).\n\nThe 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.\n\nNote too our use of `println` and the values of `Z` and `K` to produce an informative _break_ between each combination.\n\n::: {.callout-note icon=false}\n\nCan 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?\n:::\n\n::: {#62075ef2 .cell execution_count=4}\n``` {.julia .cell-code}\nfor z in Z_levels\n for k in K_levels\n\n println(\" ***> This is iteration with Z = $z and K = $k\\n\")\n\n # Define the food web\n fw = Foodweb(:niche; S = S, C = C)\n # specify the K value of the producer growth function\n\n B0 = rand(S)\n # specify model to simulate logistic growth as well as BM ratio\n params = default_model(fw, BodyMass(; Z = z), LogisticGrowth(; K = k))\n \n # number of timestamps\n t = 300\n\n out = simulate(params, B0, t)\n\n # calculate metrics\n fin_rich = richness(out)\n fin_biomass = total_biomass(out)\n s_div = shannon_diversity(out)\n\n push!(df_collect, [z, k, fin_rich, fin_biomass, s_div])\n end\nend\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n ***> This is iteration with Z = 0.1 and K = 0.1\n\n ***> This is iteration with Z = 0.1 and K = 1.0\n\n ***> This is iteration with Z = 0.1 and K = 10.0\n\n ***> This is iteration with Z = 0.1 and K = 100.0\n\n ***> This is iteration with Z = 1.0 and K = 0.1\n\n ***> This is iteration with Z = 1.0 and K = 1.0\n\n ***> This is iteration with Z = 1.0 and K = 10.0\n\n ***> This is iteration with Z = 1.0 and K = 100.0\n\n ***> This is iteration with Z = 10.0 and K = 0.1\n\n ***> This is iteration with Z = 10.0 and K = 1.0\n\n ***> This is iteration with Z = 10.0 and K = 10.0\n\n ***> This is iteration with Z = 10.0 and K = 100.0\n\n ***> This is iteration with Z = 100.0 and K = 0.1\n\n ***> This is iteration with Z = 100.0 and K = 1.0\n\n ***> This is iteration with Z = 100.0 and K = 10.0\n\n ***> This is iteration with Z = 100.0 and K = 100.0\n\n```\n:::\n:::\n\n\nWonderful. Now we are in a position to learn about two new plotting methods. First, let's look at the data frame we've created.\n\n::: {#85639f0b .cell execution_count=5}\n``` {.julia .cell-code}\ndf_collect\n```\n\n::: {.cell-output .cell-output-display execution_count=6}\n```{=html}\n
16×5 DataFrame
RowZKFinalRichnessFinalBiomassShannonDiversity
AnyAnyAnyAnyAny
10.10.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]
20.11.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]
30.110.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]
40.1100.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]
51.00.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]
61.01.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]
71.010.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]
81.0100.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]
910.00.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]
1010.01.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]
1110.010.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]
1210.0100.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]
13100.00.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]
14100.01.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]
15100.010.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]
16100.0100.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]
\n```\n:::\n:::\n\n\n#### Visualising the experiment\n\nOne 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?\n\nLet's first look at a single plot of stability\n\n::: {#9159766f .cell execution_count=6}\n``` {.julia .cell-code}\n@df df_collect plot(:K, [:FinalStability], group = :Z, \n ylabel = \"Stabilty\", \n\txlabel = \"Karrying Kapacity\",\n seriestype = [:scatter, :line],\n legend = false)\n```\n:::\n\n\nNow some new ploting tricks... 3 plots in a layout.\n\n::: {#50bb5a5c .cell execution_count=7}\n``` {.julia .cell-code}\np1 = @df df_collect plot(:K, [:FinalStability], group = :Z, \n legend = :bottomright,\n ylabel = \"Stabilty\", \n\txlabel = \"Karrying Kapacity\",\n seriestype = [:scatter, :line])\n\np2 = @df df_collect plot(:K, [:FinalBiomass], group = :Z, \n legend = :bottomright,\n ylabel = \"Biomass\", \n\txlabel = \"Karrying Kapacity\",\n seriestype = [:scatter, :line])\n \np3 = @df df_collect plot(:K, [:FinalRichness], group = :Z, \n legend = :bottomright,\n ylabel = \"Richness\", \n\txlabel = \"Karrying Kapacity\",\n seriestype = [:scatter, :line])\n\n# create a layout of 3 graphs stacked on top of each other.\nplot(p1, p2, p3, layout=(3,1), legend = false)\n```\n:::\n\n\n### Interpretation!\n\n#### Challenge - can you get the number of extinctions into the data frame?\n\n### Experiment 2: The Functional Response\n\nThe 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.\n\nA 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.\n\nA 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**.\n\nFinally, 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.\n\n::: {.callout-tip icon=false}\n\nRemember that the logistic growth equation, with a carying capacity specified, is also a source of _density dependent negative feedback_\n:::\n\n::: {.callout-tip icon=false}\n\nThe 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.\n:::\n\n::: {#5eeb3011 .cell execution_count=8}\n``` {.julia .cell-code}\nf_t1(n) = 0.5*n\nf_t2(n) = 0.5*n/(0.2+0.01*n)\nf_t3(n) = 0.5*n^2/(10 + 0.01*n^2)\n\nplot(f_t1, 0, 100, label = \"Type I\")\nplot!(f_t2, 0, 100, label = \"Type II\")\nplot!(f_t3, 0, 100, label = \"Type III\")\n```\n\n::: {.cell-output .cell-output-display execution_count=7}\n```{=html}\n\n\n\n \n \n \n\n\n\n \n \n \n\n\n\n \n \n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n```\n:::\n:::\n\n\n#### How does the BEFW make a functional response?\n\nThere 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.\n\nThe 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.\n\nThe 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.\n\n#### Creating Type II vs. Type III with the Bioenergetic response\n\nLet'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.\n\n::: {#55111e82 .cell execution_count=9}\n``` {.julia .cell-code}\nRandom.seed!(12352)\n\n# fixed parameters\nS = 20\nC = 0.15\n\n# set the hill exponent to move from Type II to Type III)\nh_levels = [1.0, 1.1, 1.25, 2.0]\n\n# set collecting data frame \n# we will look at how Richness, Biomass and Stability are affected by the hill exponent\ndf_collect_h = DataFrame(h = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])\n\n# create look across values of h\nfor h in h_levels \n println(\"***> This is iteration with h = $h\\n\")\n \n # make the network\n # Note that we specify the Environment, but there is no K or T set (using defaults)\n # Note the new BioenergeticResponse function\n fw_h = Foodweb(:niche; S = S, C = C)\n \n # set body sizes and parameters \n B0 = rand(S)\n params = default_model(fw_h)\n\n # here we now update the exponent of the hill function\n params.hill_exponent = h\n\n # specify number of time steps\n t = 300\n\n # simulate\n sim_niche = simulate(params, B0, t)\n\n # collect data \n fin_rich = richness(sim_niche)\n fin_bio = total_biomass(sim_niche)\n s_div = shannon_diversity(sim_niche)\n\n push!(df_collect_h, [h, fin_rich, fin_bio, s_div])\nend\n\ndf_collect_h\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n***> This is iteration with h = 1.0\n\n***> This is iteration with h = 1.1\n\n***> This is iteration with h = 1.25\n\n***> This is iteration with h = 2.0\n\n```\n:::\n\n::: {.cell-output .cell-output-display execution_count=8}\n```{=html}\n
4×4 DataFrame
RowhFinalRichnessFinalBiomassShannonDiversity
AnyAnyAnyAny
11.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]
21.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]
31.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]
42.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]
\n```\n:::\n:::\n\n\nNow, we can visualise these data\n\n::: {#be0fe852 .cell execution_count=10}\n``` {.julia .cell-code}\n# Visualize the results\np1_h = @df df_collect_h plot(:h, [:FinalStability],\n legend = :bottomright,\n ylabel = \"Stability\",\n xlabel = \"Functional response\",\n seriestype = [:scatter, :line])\n\np2_h = @df df_collect_h plot(:h, [:FinalBiomass],\n legend = :bottomright,\n ylabel = \"Biomass\",\n xlabel = \"Functional response\",\n seriestype = [:scatter, :line])\n\np3_h = @df df_collect_h plot(:h, [:FinalRichness],\n legend = :bottomright,\n ylabel = \"Richness\",\n xlabel = \"Functional response\",\n seriestype = [:scatter, :line])\n\nplot(p1_h, p2_h, p3_h, layout=(3,1), legend = false, size = (1000, 1000))\n```\n:::\n\n\n#### INTERPRETATION?\n\nWhat can you see happening as we move away from the destabilising Type II functional response?\n\nCan 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.\n\n::: {#f8aca1ed .cell execution_count=11}\n``` {.julia .cell-code}\n # make the network\n fw_h = Foodweb(:niche; S = S, C = C)\n\n # set body sizes and parameters \n B0 = rand(S)\n params = default_model(fw_h, LogisticGrowth(; K = k))\n\n # update the exponent of the hill function\n params.hill_exponent = h\n```\n:::\n\n\n### Experiment 3: What is Z\n\nOne 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. \n\nIn 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.\n\nIn 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):\n\n$M_C = Z^(T-1)$\n\n[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()`\n\n::: {#b1b7b1ce .cell execution_count=12}\n``` {.julia .cell-code}\nRandom.seed!(12352)\n\n# fixed parameters\nS = 20\nC = 0.15\n\n# set the PPRM\nz_levels= [0.1, 1, 10, 100]\n\n# set collecting data frame \n# we will look at how Richness, Biomass and Stability are affected by the hill exponent\ndf_collect_z = DataFrame(z = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])\n\n# create look across values of h\nfor z in z_levels \n println(\"***> This is iteration with z = $z\\n\")\n \n # make the network\n # Note that we specify the Environment, but there is no K or T set (using defaults)\n # Note Z is specified when building the FoodWeb() network\n fw_z = Foodweb(:niche; S = S, C = C)\n \n # set body sizes and parameters \n B0 = rand(S)\n params = default_model(fw_z, BodyMass(; Z = z))\n\n # specify number of time steps\n t = 300\n\n # simulate\n out_z = simulate(params, B0, t)\n\n # collect data \n fin_rich = richness(out_z)\n fin_bio = total_biomass(out_z)\n s_div = shannon_diversity(out_z)\n\n push!(df_collect_z, [z, fin_rich, fin_bio, s_div])\nend\n\ndf_collect_z\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n***> This is iteration with z = 0.1\n\n***> This is iteration with z = 1.0\n\n***> This is iteration with z = 10.0\n\n***> This is iteration with z = 100.0\n\n```\n:::\n\n::: {.cell-output .cell-output-display execution_count=9}\n```{=html}\n
4×4 DataFrame
RowzFinalRichnessFinalBiomassShannonDiversity
AnyAnyAnyAny
10.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]
21.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]
310.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]
4100.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]
\n```\n:::\n:::\n\n\nAs 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.\n\n::: {#ab493d43 .cell execution_count=13}\n``` {.julia .cell-code}\n# Visualize the results\np1_z = @df df_collect_z plot(:z, [:FinalStability],\n legend = :bottomright,\n ylabel = \"Stability\",\n xlabel = \"Z value (PPMR)\",\n seriestype = [:scatter, :line])\n\np2_z = @df df_collect_z plot(:z, [:FinalBiomass],\n legend = :bottomright,\n ylabel = \"Biomass\",\n xlabel = \"Z value (PPMR)\",\n seriestype = [:scatter, :line])\n\np3_z = @df df_collect_z plot(:z, [:FinalRichness],\n legend = :bottomright,\n ylabel = \"Richness\",\n xlabel = \"Z value (PPMR)\",\n seriestype = [:scatter, :line])\n\nplot(p1_z, p2_z, p3_z, layout=(3,1), legend = false, size = (1000, 1000))\n```\n:::\n\n\n#### Challenge\n\nCould 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.\n\n### Experiment 4: Manipulating Competition among producers\n\nOur 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.\n\nWe 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.\n\nThe 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. \n\nWe 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$. \n\nWhat 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$.\n\n::: {#12de3283 .cell execution_count=14}\n``` {.julia .cell-code}\nS = 20 # define the number of species\nC = 0.2 # define the connectance (complexity) of the network\nZ = 100 # Predator Prey Mass Ratio\n\n# here we set the \ninterspecific_vals = 0.8:0.05:1.2 # a set of (9) values between 0.8 and 1.2 in steps of 0.05\n# collect(0.8:0.05:1.2) # see them if you want to\n\n\n# set collecting data frame \n# we will look at how Richness, Biomass and Stability are affected by the hill exponent\ndf_collect_comp = DataFrame(InterComp = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])\n\nfor alpha_ij in interspecific_vals\n println(\"***> This is iteration with alpha_ij = $alpha_ij\\n\")\n \n # this will make always the same network and body mass\n Random.seed!(123)\n foodweb = Foodweb(:niche; S = S, C = C, Z = Z)\n\n # enhance detail in the network\n # specify K\n LogisticGrowth(foodweb; K = 10)\n # we fix intraspecific = 1 and vary interspecific\n ProducerCompetiton(foodweb; αii = 1.0, αij = 1)\n # set the hill exponent to 1 (type II Functional Response)\n BioenergeticResponse(foodweb, h = 1)\n\n # define parameters with extras\n params_comp = default_model(foodweb, BodyMass(; Z = z))\n\n # set bodymass\n B0 = rand(S)\n\n # simulate\n # note verbose = false ; remove this to see extinction detail for each sim\n out_comp = simulate(params_comp, B0, verbose = false)\n\n # generate stats and add to collector\n # collect data \n fin_rich = richness(out_comp)\n fin_bio = biomass(out_comp).total\n stab = community_cv(out_comp)\n\n push!(df_collect_comp, [alpha_ij, fin_rich, fin_bio, stab])\nend\n\ndf_collect_comp\n```\n:::\n\n\nLet'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.\n\n::: {#7d288b7a .cell execution_count=15}\n``` {.julia .cell-code}\np1 = @df df_collect_comp plot(:InterComp, [:FinalRichness],\n ylabel = \"Richness\",\n xlabel = \"alpha_ij\")\np2 = @df df_collect_comp plot(:InterComp, [:FinalBiomass],\n ylabel = \"Biomass\",\n xlabel = \"alpha_ij\")\np3 = @df df_collect_comp plot(:InterComp, [:FinalStability],\n ylabel = \"Stability\",\n xlabel = \"alpha_ij\")\n\nplot(p1, p2, p3, layout = (3,1), legend = false)\n```\n:::\n\n\n#### Challenge - Competition\n\nPerhaps consider expanding the code above to assess one of these?\n\nIs this pattern sensitive to specie richness?\nIs this pattern sensitive to the functional response?\nIs this pattern sensitive to the PPMR?\nIs this pattern sensitive to values of K?\n\n## Experiment 5: Multiple networks (replicates)\n\nTo 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...\n\n", + "engine": "julia", + "markdown": "---\ntitle: \"Tutorial 10: Complex Experiments with the END\"\ndate: last-modified\nauthor: \"Danet and Becks, based on originals by Delmas and Griffiths\"\nformat:\n html:\n embed-resources: true\ntitle-block-banner: true\nengine: julia\n---\n\n\n\n\n::: {.callout-caution}\n## Warning\n\nStability was replaced by Shannon diversity - need to review context as well as if we want that or rather re-integrate stability\n:::\n\nThe 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.\n\nIn 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.\n\nYou 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`).\n\n\n\n\n::: {#2 .cell execution_count=1}\n``` {.julia .cell-code}\nusing Random, Plots, Distributions, DataFrames, StatsPlots\nusing EcologicalNetworksDynamics\n```\n:::\n\n\n\n\n\n\n### Experiment 1: Carrying Capacity and the Predator Prey Mass Ratio\n\nNow 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. \n\nZ 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. \n\nThe 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.\n\n\n\n\n\n::: {#4 .cell execution_count=1}\n``` {.julia .cell-code}\n#Fixed Parameters\nS = 20\nC = 0.15\n\n# Variable Parameters\nZ_levels = [0.1, 1, 10, 100]\nK_levels = [0.1, 1, 10, 100]\n\n# run this to get same results as in the document\nRandom.seed!(123)\n```\n\n::: {.cell-output .cell-output-display execution_count=1}\n```\nTaskLocalRNG()\n```\n:::\n:::\n\n\n\n\n\n\nNow, lets set up the collecting data frame.\n\n\n\n\n::: {#6 .cell execution_count=1}\n``` {.julia .cell-code}\ndf_collect = DataFrame(Z = [], K = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])\n```\n\n::: {.cell-output .cell-output-display execution_count=1}\n```{=html}\n
0×5 DataFrame
RowZKFinalRichnessFinalBiomassShannonDiversity
AnyAnyAnyAnyAny
\n```\n:::\n:::\n\n\n\n\n\n\nNow, 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).\n\nThe 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.\n\nNote too our use of `println` and the values of `Z` and `K` to produce an informative _break_ between each combination.\n\n::: {.callout-note icon=false}\n\nCan 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?\n:::\n\n\n\n\n\n\n::: {#8 .cell execution_count=1}\n``` {.julia .cell-code}\nfor z in Z_levels\n for k in K_levels\n\n println(\" ***> This is iteration with Z = $z and K = $k\\n\")\n\n # Define the food web\n fw = Foodweb(:niche; S = S, C = C)\n # specify the K value of the producer growth function\n\n B0 = rand(S)\n # specify model to simulate logistic growth as well as BM ratio\n params = default_model(fw, BodyMass(; Z = z), LogisticGrowth(; K = k))\n \n # number of timestamps\n t = 300\n\n out = simulate(params, B0, t)\n\n # calculate metrics\n fin_rich = richness(out)\n fin_biomass = total_biomass(out)\n s_div = shannon_diversity(out)\n\n push!(df_collect, [z, k, fin_rich, fin_biomass, s_div])\n end\nend\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n ***> This is iteration with Z = 0.1 and K = 0.1\n\n ***> This is iteration with Z = 0.1 and K = 1.0\n\n ***> This is iteration with Z = 0.1 and K = 10.0\n\n ***> This is iteration with Z = 0.1 and K = 100.0\n\n ***> This is iteration with Z = 1.0 and K = 0.1\n\n ***> This is iteration with Z = 1.0 and K = 1.0\n\n ***> This is iteration with Z = 1.0 and K = 10.0\n\n ***> This is iteration with Z = 1.0 and K = 100.0\n\n ***> This is iteration with Z = 10.0 and K = 0.1\n\n ***> This is iteration with Z = 10.0 and K = 1.0\n\n ***> This is iteration with Z = 10.0 and K = 10.0\n\n ***> This is iteration with Z = 10.0 and K = 100.0\n\n ***> This is iteration with Z = 100.0 and K = 0.1\n\n ***> This is iteration with Z = 100.0 and K = 1.0\n\n ***> This is iteration with Z = 100.0 and K = 10.0\n\n ***> This is iteration with Z = 100.0 and K = 100.0\n\n```\n:::\n:::\n\n\n\n\n\n\nWonderful. Now we are in a position to learn about two new plotting methods. First, let's look at the data frame we've created.\n\n\n\n\n::: {#10 .cell execution_count=1}\n``` {.julia .cell-code}\ndf_collect\n```\n\n::: {.cell-output .cell-output-display execution_count=1}\n```{=html}\n
16×5 DataFrame
RowZKFinalRichnessFinalBiomassShannonDiversity
AnyAnyAnyAnyAny
10.10.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]
20.11.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]
30.110.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]
40.1100.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]
51.00.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]
61.01.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]
71.010.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]
81.0100.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]
910.00.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]
1010.01.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]
1110.010.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]
1210.0100.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]
13100.00.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]
14100.01.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]
15100.010.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]
16100.0100.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 … 465.906, 464.433, 463.337, 461.4, 458.545, 455.022, 451.776, 450.914, 456.694, 471.903][18.3089, 18.2607, 18.15, 17.9446, 17.542, 16.8141, 15.6206, 13.9258, 12.0175, 9.86971 … 7.1311, 7.20354, 7.11395, 6.98256, 6.85529, 6.73983, 6.6778, 6.78037, 7.20202, 7.95913]
\n```\n:::\n:::\n\n\n\n\n\n\n#### Visualising the experiment\n\nOne 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?\n\nLet's first look at a single plot of stability\n\n\n\n::: {#12 .cell execution_count=0}\n``` {.julia .cell-code}\n@df df_collect plot(:K, [:FinalStability], group = :Z, \n ylabel = \"Stabilty\", \n\txlabel = \"Karrying Kapacity\",\n seriestype = [:scatter, :line],\n legend = false)\n```\n:::\n\n\n\n\n\n\nNow some new ploting tricks... 3 plots in a layout.\n\n\n\n\n::: {#14 .cell execution_count=0}\n``` {.julia .cell-code}\np1 = @df df_collect plot(:K, [:FinalStability], group = :Z, \n legend = :bottomright,\n ylabel = \"Stabilty\", \n\txlabel = \"Karrying Kapacity\",\n seriestype = [:scatter, :line])\n\np2 = @df df_collect plot(:K, [:FinalBiomass], group = :Z, \n legend = :bottomright,\n ylabel = \"Biomass\", \n\txlabel = \"Karrying Kapacity\",\n seriestype = [:scatter, :line])\n \np3 = @df df_collect plot(:K, [:FinalRichness], group = :Z, \n legend = :bottomright,\n ylabel = \"Richness\", \n\txlabel = \"Karrying Kapacity\",\n seriestype = [:scatter, :line])\n\n# create a layout of 3 graphs stacked on top of each other.\nplot(p1, p2, p3, layout=(3,1), legend = false)\n```\n:::\n\n\n\n\n\n\n### Interpretation!\n\n#### Challenge - can you get the number of extinctions into the data frame?\n\n### Experiment 2: The Functional Response\n\nThe 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.\n\nA 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.\n\nA 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**.\n\nFinally, 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.\n\n::: {.callout-tip icon=false}\n\nRemember that the logistic growth equation, with a carying capacity specified, is also a source of _density dependent negative feedback_\n:::\n\n::: {.callout-tip icon=false}\n\nThe 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.\n:::\n\n\n\n\n::: {#16 .cell execution_count=1}\n``` {.julia .cell-code}\n\nf_t1(n) = 0.5*n\nf_t2(n) = 0.5*n/(0.2+0.01*n)\nf_t3(n) = 0.5*n^2/(10 + 0.01*n^2)\n\nplot(f_t1, 0, 100, label = \"Type I\")\nplot!(f_t2, 0, 100, label = \"Type II\")\nplot!(f_t3, 0, 100, label = \"Type III\")\n```\n\n::: {.cell-output .cell-output-display execution_count=1}\n```{=html}\n\n\n\n \n \n \n\n\n\n \n \n \n\n\n\n \n \n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n```\n:::\n:::\n\n\n\n\n\n\n#### How does the BEFW make a functional response?\n\nThere 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.\n\nThe 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.\n\nThe 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.\n\n#### Creating Type II vs. Type III with the Bioenergetic response\n\nLet'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.\n\n\n\n\n::: {#18 .cell execution_count=1}\n``` {.julia .cell-code}\n\nRandom.seed!(12352)\n\n# fixed parameters\nS = 20\nC = 0.15\n\n# set the hill exponent to move from Type II to Type III)\nh_levels = [1.0, 1.1, 1.25, 2.0]\n\n# set collecting data frame \n# we will look at how Richness, Biomass and Stability are affected by the hill exponent\ndf_collect_h = DataFrame(h = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])\n\n# create look across values of h\nfor h in h_levels \n println(\"***> This is iteration with h = $h\\n\")\n \n # make the network\n # Note that we specify the Environment, but there is no K or T set (using defaults)\n # Note the new BioenergeticResponse function\n fw_h = Foodweb(:niche; S = S, C = C)\n \n # set body sizes and parameters \n B0 = rand(S)\n params = default_model(fw_h)\n\n # here we now update the exponent of the hill function\n params.hill_exponent = h\n\n # specify number of time steps\n t = 300\n\n # simulate\n sim_niche = simulate(params, B0, t)\n\n # collect data \n fin_rich = richness(sim_niche)\n fin_bio = total_biomass(sim_niche)\n s_div = shannon_diversity(sim_niche)\n\n push!(df_collect_h, [h, fin_rich, fin_bio, s_div])\nend\n\ndf_collect_h\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n***> This is iteration with h = 1.0\n\n***> This is iteration with h = 1.1\n\n***> This is iteration with h = 1.25\n\n***> This is iteration with h = 2.0\n\n```\n:::\n\n::: {.cell-output .cell-output-display execution_count=1}\n```{=html}\n
4×4 DataFrame
RowhFinalRichnessFinalBiomassShannonDiversity
AnyAnyAnyAny
11.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]
21.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]
31.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]
42.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]
\n```\n:::\n:::\n\n\n\n\n\n\nNow, we can visualise these data\n\n\n\n\n::: {#20 .cell execution_count=0}\n``` {.julia .cell-code}\n# Visualize the results\np1_h = @df df_collect_h plot(:h, [:FinalStability],\n legend = :bottomright,\n ylabel = \"Stability\",\n xlabel = \"Functional response\",\n seriestype = [:scatter, :line])\n\np2_h = @df df_collect_h plot(:h, [:FinalBiomass],\n legend = :bottomright,\n ylabel = \"Biomass\",\n xlabel = \"Functional response\",\n seriestype = [:scatter, :line])\n\np3_h = @df df_collect_h plot(:h, [:FinalRichness],\n legend = :bottomright,\n ylabel = \"Richness\",\n xlabel = \"Functional response\",\n seriestype = [:scatter, :line])\n\nplot(p1_h, p2_h, p3_h, layout=(3,1), legend = false, size = (1000, 1000))\n```\n:::\n\n\n\n\n\n\n#### INTERPRETATION?\n\nWhat can you see happening as we move away from the destabilising Type II functional response?\n\nCan 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.\n\n\n\n\n::: {#22 .cell execution_count=0}\n``` {.julia .cell-code}\n # make the network\n fw_h = Foodweb(:niche; S = S, C = C)\n\n # set body sizes and parameters \n B0 = rand(S)\n params = default_model(fw_h, LogisticGrowth(; K = k))\n\n # update the exponent of the hill function\n params.hill_exponent = h\n```\n:::\n\n\n\n\n\n\n### Experiment 3: What is Z\n\nOne 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. \n\nIn 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.\n\nIn 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):\n\n$M_C = Z^(T-1)$\n\n[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()`\n\n\n\n\n::: {#24 .cell execution_count=1}\n``` {.julia .cell-code}\nRandom.seed!(12352)\n\n# fixed parameters\nS = 20\nC = 0.15\n\n# set the PPRM\nz_levels= [0.1, 1, 10, 100]\n\n# set collecting data frame \n# we will look at how Richness, Biomass and Stability are affected by the hill exponent\ndf_collect_z = DataFrame(z = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])\n\n# create look across values of h\nfor z in z_levels \n println(\"***> This is iteration with z = $z\\n\")\n \n # make the network\n # Note that we specify the Environment, but there is no K or T set (using defaults)\n # Note Z is specified when building the FoodWeb() network\n fw_z = Foodweb(:niche; S = S, C = C)\n \n # set body sizes and parameters \n B0 = rand(S)\n params = default_model(fw_z, BodyMass(; Z = z))\n\n # specify number of time steps\n t = 300\n\n # simulate\n out_z = simulate(params, B0, t)\n\n # collect data \n fin_rich = richness(out_z)\n fin_bio = total_biomass(out_z)\n s_div = shannon_diversity(out_z)\n\n push!(df_collect_z, [z, fin_rich, fin_bio, s_div])\nend\n\ndf_collect_z\n```\n\n::: {.cell-output .cell-output-stdout}\n```\n***> This is iteration with z = 0.1\n\n***> This is iteration with z = 1.0\n\n***> This is iteration with z = 10.0\n\n***> This is iteration with z = 100.0\n\n```\n:::\n\n::: {.cell-output .cell-output-display execution_count=1}\n```{=html}\n
4×4 DataFrame
RowzFinalRichnessFinalBiomassShannonDiversity
AnyAnyAnyAny
10.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]
21.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]
310.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]
4100.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]
\n```\n:::\n:::\n\n\n\n\n\n\nAs 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.\n\n\n\n\n::: {#26 .cell execution_count=0}\n``` {.julia .cell-code}\n# Visualize the results\np1_z = @df df_collect_z plot(:z, [:FinalStability],\n legend = :bottomright,\n ylabel = \"Stability\",\n xlabel = \"Z value (PPMR)\",\n seriestype = [:scatter, :line])\n\np2_z = @df df_collect_z plot(:z, [:FinalBiomass],\n legend = :bottomright,\n ylabel = \"Biomass\",\n xlabel = \"Z value (PPMR)\",\n seriestype = [:scatter, :line])\n\np3_z = @df df_collect_z plot(:z, [:FinalRichness],\n legend = :bottomright,\n ylabel = \"Richness\",\n xlabel = \"Z value (PPMR)\",\n seriestype = [:scatter, :line])\n\nplot(p1_z, p2_z, p3_z, layout=(3,1), legend = false, size = (1000, 1000))\n```\n:::\n\n\n\n\n\n\n#### Challenge\n\nCould 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.\n\n### Experiment 4: Manipulating Competition among producers\n\nOur 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.\n\nWe 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.\n\nThe 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. \n\nWe 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$. \n\nWhat 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$.\n\n\n\n\n\n::: {#28 .cell execution_count=0}\n``` {.julia .cell-code}\nS = 20 # define the number of species\nC = 0.2 # define the connectance (complexity) of the network\nZ = 100 # Predator Prey Mass Ratio\nt = 300 # time steps\n\n# here we set the \ninterspecific_vals = 0.8:0.05:1.2 # a set of (9) values between 0.8 and 1.2 in steps of 0.05\n# collect(0.8:0.05:1.2) # see them if you want to\n\n\n# set collecting data frame \n# we will look at how Richness, Biomass and Shannon Diversity are affected by interspecific competition\ndf_collect_comp = DataFrame(InterComp = [], FinalRichness = [], FinalBiomass = [], ShannonDiversity = [])\n\nfor alpha_ij in interspecific_vals\n println(\"***> This is iteration with alpha_ij = $alpha_ij\\n\")\n \n # this will make always the same network and body mass\n Random.seed!(123)\n foodweb = Foodweb(:niche; S = S, C = C, Z = Z)\n\n # enhance detail in the network\n #We fix intraspecific = 1 (diag) and vary interspecific (off)\n p = ProducersCompetition(foodweb; diag = 1.0, off = alpha_ij)\n #Specify K and competition in LogisticGrowth\n r = LogisticGrowth(K = 10, producers_competition = p)\n\n #Set the hill exponent to 1 (type II Functional Response)\n b = BioenergeticResponse(h = 1)\n\n #Define parameters with extras\n params_comp = default_model(foodweb, BodyMass(; Z = Z), b, r)\n\n # set initial biomasses\n B0 = rand(S)\n\n # simulate\n # note verbose = false ; remove this to see extinction detail for each sim\n out_comp = simulate(params_comp, B0, t, verbose = false)\n\n # generate stats and add to collector\n # collect data \n fin_rich = richness(out_comp)\n fin_bio = total_biomass(out_comp)\n s_div = shannon_diversity(out_comp)\n\n push!(df_collect_comp, [alpha_ij, fin_rich, fin_bio, s_div])\nend\n\ndf_collect_comp\n```\n:::\n\n\n\n\n\n\nLet'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.\n\n\n\n\n\n::: {#30 .cell execution_count=0}\n``` {.julia .cell-code}\np1 = @df df_collect_comp plot(:InterComp, [:FinalRichness],\n ylabel = \"Richness\",\n xlabel = \"alpha_ij\")\np2 = @df df_collect_comp plot(:InterComp, [:FinalBiomass],\n ylabel = \"Biomass\",\n xlabel = \"alpha_ij\")\np3 = @df df_collect_comp plot(:InterComp, [:ShannonDiversity],\n ylabel = \"Shannon Diversity\",\n xlabel = \"alpha_ij\")\n\nplot(p1, p2, p3, layout = (3,1), legend = false)\n```\n:::\n\n\n\n\n\n\n#### Challenge - Competition\n\nPerhaps consider expanding the code above to assess one of these?\n\nIs this pattern sensitive to species richness?\nIs this pattern sensitive to the functional response?\nIs this pattern sensitive to the PPMR?\nIs this pattern sensitive to values of K?\n\n## Experiment 5: Multiple networks (replicates)\n\nTo 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...\n\n", "supporting": [ - "10_complex_experiments_end_files" + "10_complex_experiments_end_files/figure-html" ], "filters": [], "includes": {