31:T456,viewof e1_plot = { var simline= vl.markLine({strokeWidth: 1, interpolate: 'step-after', }) .data(timecourse) .encode( vl.x().fieldQ("t") .axis({titleFont:"Helvetica Neue"}), vl.y().fieldQ("N1") .title("N₁") .scale({domain: [0,K]}) .axis({titleFont:"Helvetica Neue"}), vl.color().value(accent_color) ) .width(600); var delta_arrows = vl.markPoint({shape:"arrow"}) .data(delta_N1_per_N1) .transform( vl.calculate("(abs(datum.deltaN1))").as("absdN1"), vl.calculate("datum.deltaN1 > 0 ? 180 : 0").as("sdN1")) .encode( vl.size().fieldQ("absdN1").scale({range: [0, 280]}).legend(false), //vl.x().fieldQ("deltaN1"), vl.angle().fieldQ("sdN1").scale({domain: [0, 360], range: [180, -180]}), vl.fill().value(accent_color), vl.opacity().fieldQ("absdN1"), vl.color().value(false), vl.y().fieldQ("N1") .title(false) .axis(false), ) .config({view: {stroke: null}}) return vl.hconcat(simline,delta_arrows).config({ concat: { spacing: 0 }, view: {stroke: null}}).render() }32:T4a8,viewof e2_plot = { var e2_simline = vl.markLine({strokeWidth: 1, interpolate: 'step-after', }) .data(e2_timecourse) .encode( vl.x().fieldQ("t") .axis({titleFont:"Helvetica Neue"}), vl.y().fieldQ("N1") .title("N₁") .scale({domain: [_.minBy(e2_timecourse, "N1").N1,_.max([K,_.maxBy(e2_timecourse, "N1").N1])]}) .axis({titleFont:"Helvetica Neue"}), vl.color().value(accent_color), ) .width(600) var delta_arrows = vl.markPoint({shape:"arrow"}) .data(e2_delta_N1_per_N1) .transform( vl.calculate("(abs(datum.deltaN1))").as("absdN1"), vl.calculate("datum.deltaN1 > 0 ? 180 : 0").as("sdN1")) .encode( vl.size().fieldQ("absdN1").scale({range: [0, 280]}).legend(false), //vl.x().fieldQ("deltaN1"), vl.angle().fieldQ("sdN1").scale({domain: [0, 360], range: [180, -180]}), vl.fill().value(accent_color), vl.opacity().fieldQ("absdN1"), vl.color().value(false), vl.y().fieldQ("N1") .title(false) .axis(false), ) .config({view: {stroke: null}}) return vl.hconcat(e2_simline,delta_arrows).config({ concat: { spacing: 0 }, view: {stroke: null}}).render() }33:[] 2f:{"title":"Species coexistance and variablity.","cells":[{"id":0,"value":"md`\n# Species coexistance and variablity. `","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"md `A notebook to keep track of different mechanisms that facilitate species coexistance and are stchastisity-driven. _WIP_`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"md`## Long-term average growth rate`","pinned":false,"hidden":false,"mode":"ojs"},{"id":20,"value":"md `This is based on one species having a better arithmetic-mean reproductive fitness but more variability. Because the payoff of reproductive fitness is non-linear and concave, the variability decreases it's long-term average growth rate when Resident. Thus the lower-mean-fitness species can invade.` ","pinned":false,"hidden":false,"mode":"ojs"},{"id":32,"value":"md `Let's go though the example given [here](https://doi.org/10.1111/j.1461-0248.2011.01643.x). \nThis is based on the Lottery model developed by [Chesson & Warner](https://doi.org/10.1086/283778).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":39,"value":"md `> «Consider two species, denoted 1 and 2, competing for vacant space following the random death of adults. The recruitment probability of a species is simply its relative abundance in the pool of seeds reaching the vacant site.»`","pinned":false,"hidden":false,"mode":"ojs"},{"id":51,"value":"md `The update rule for species 1 is: `","pinned":false,"hidden":false,"mode":"ojs"},{"id":53,"value":"tex `N_{1,t+1} = (1-d_1) N_{1,t} + (d_1 N_{1,t} + d_2 N_{2,t}) \\frac{f_1(t) N_{1,t}}{f_1(t) N_{1,t} + f_2(t) N_{2,t}}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":68,"value":"texmd ` $d_i$ is the death rate and $f_i$ the fecundity of species $i$. Because the sum of $N_1$ and $N_2$ is conserved, we don't need $N_{2,t+1}$.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":404,"value":"texmd `Let's go through the update rule $N_{1,t+1}$ quickly to make sure we understand what happens here. \nThe first term $(1-d_1) N_{1,t}$ is the fraction of species 1 that survives from _t_ to _t+1_. It's the survival probalility of species 1 $(1-d_1)$ times $N_{1,t}$ the abundance of species 1.\nThe last term $(d_1 N_{1,t} + d_2 N_{2,t}) \\frac{f_1(t) N_{1,t}}{f_1(t) N_{1,t} + f_2(t) N_{2,t}}$ is the fraction of newly “born” members of species 1 replacing the dead of both species. $(d_1 N_{1,t} + d_2 N_{2,t})$ is the fraction of the combined population dying from _t_ to _t+1_. The term $f_i(t) N_{i,t}$ is the rate at which species $i$ can produce offspring. The fraction $\\frac{f_1(t) N_{1,t}}{f_1(t) N_{1,t} + f_2(t) N_{2,t}}$ is fraction of vacant space species 1 can succesfully ocupy in competition with species 2.` ","pinned":false,"hidden":false,"mode":"ojs"},{"id":63,"value":"texmd `Let's assume species 2 has the advantage of a higher average fecundity ($\\bar{f_2} > \\bar{f_1}$) and only $f_2(t)$ fluctuates in response to the enviroment. For the sake of simplicity we assume the most simple model possible for fluctuations: Two “Enviroments” drawn with equal probability. The death rates are equal $d = d_1 = d_2$. Note that this is the mean-field approximation (valid if the number of individuals $\\rightarrow \\infty\n$) and does not include variablity from the random drawing of deaths and recruitment. Since the $N_i$ are not natural numbers, we set the maximum capacity $K = N_1 + N_2 = 1$.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":27,"value":"md`Okay, simulate this!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":96,"value":"viewof d = slider({\n min: 0,\n max: 1,\n step: 0.05,\n value: 0.1,\n description:\n \"d - the death rate of the species\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":103,"value":"viewof f2divf1 = slider({\n min: 1,\n max: 1.1,\n step: 0.01,\n value: 1.02,\n description:\n \"f₂/f₁ - the relative advantage in fecundity of the species 2\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":376,"value":"viewof deltaf2 = slider({\n min: 0,\n max: 2,\n step: 0.1,\n value: 0.4,\n description:\n \"Δf₂ - the relative stochastic variation of f₂\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":252,"value":"viewof N1_0 = slider({\n min: 0,\n max: 1,\n step: 0.01,\n value: 0.01,\n description: \"N₁(t=0) - the inital value of the abundance N₁ of species 1\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":349,"value":"viewof run_button = button({value: \"Run\", description: \"Simulate a new trajectory!\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":178,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":301,"value":"md `Okay, this seems to work! For a certain parameter regime, both species are able to recover when close to extinction. Let's look why this is possible. `","pinned":false,"hidden":false,"mode":"ojs"},{"id":503,"value":"texmd `To be able to recover, the growth $\\Delta N_{1} = N_{i,t+1} - N_{i,t}$ has to be positive on average for small $N_{i,t} \\ll K$. Let's look at the function $\\overline{\\Delta N_{1}(N_{1,t})}$. _This will change if you modify the simulation parameters above. The arrows in the plot above show the the direction and strength of $\\overline{\\Delta N_{1}(N_{1,t})}$ as well._`","pinned":false,"hidden":false,"mode":"ojs"},{"id":453,"value":"viewof deltline = {\n \n var max_deltaN1 = _.max([_.maxBy(delta_N1_per_N1, 'deltaN1').deltaN1, \n -_.minBy(delta_N1_per_N1, 'deltaN1').deltaN1]);\n \n console.log(max_deltaN1)\n var l = vl.markLine({interpolate: 'monotone'})\n .data(delta_N1_per_N1)\n .encode(\n vl.x().fieldQ(\"N1\")\n .title(\"N₁\")\n .axis({titleFont:\"Helvetica Neue\"}),\n vl.y().fieldQ(\"deltaN1\")\n .title(\"ΔN₁\")\n .axis({titleFont:\"Helvetica Neue\"})\n .scale({domain: [-1.05*max_deltaN1,1.05*max_deltaN1]}),\n vl.color().value(accent_color)\n )\n .width(350).height(150).render()\n return l\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":560,"value":"texmd `We see: If the simulation produces stable trajectories, $\\overline{\\Delta N_{1}(N_{1,t})}$ is $>0$ for $N_{1,t} \\approx 0$ and $<0$ for $N_{1,t} \\approx 1$!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":619,"value":"md `### The Storage Effect`","pinned":false,"hidden":false,"mode":"ojs"},{"id":626,"value":"md `The example above can also be framed as the storage effect. Even if one species is favored during one generation it's total loss is bounded because the adult population has a fixed death rate. The old, (by chance) more successful generation is _stored_ and can fuel future growth. Overlapping generations are a nessecary requirement.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":681,"value":"md `### Do the species need to be different?`","pinned":false,"hidden":false,"mode":"ojs"},{"id":687,"value":"md `Let's have two identical species both with a random component.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":703,"value":"viewof e2_d = slider({\n min: 0,\n max: 1,\n step: 0.05,\n value: 0.1,\n description:\n \"d - the death rate of the species\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":708,"value":"viewof e2_delta = slider({\n min: 0.1,\n max: 1.9,\n step: 0.1,\n value: 0.3,\n description:\n \"Δf - the relative stochastic variation of both f\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":733,"value":"viewof e2_f2divf1 = slider({\n min: 1,\n max: 1.05,\n step: 0.005,\n value: 1.01,\n precision: 4,\n description:\n \"f₂/f₁ - the relative advantage in fecundity of the species 2\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":715,"value":"viewof e2_N1_0 = slider({\n min: 0,\n max: 1,\n step: 0.01,\n value: 0.01,\n description:\n \"N₁(t=0) - the inital value of the abundance N₁ of species 1\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":725,"value":"viewof e2_run_button = button({value: \"Run\", description: \"Simulate a new trajectory!\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":719,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":884,"value":"texmd `Works too! The arrows show the direction and strenght of $\\Delta N_1$.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":70,"value":"md `## Appendix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":545,"value":"md `### simulation`","pinned":false,"hidden":false,"mode":"ojs"},{"id":122,"value":"N1_step = (N1) => (1-d)*N1+(d*K)*(f1*N1)/(f1*N1+f2()*(K-N1))","pinned":true,"hidden":false,"mode":"ojs"},{"id":148,"value":"f2 = () => {\n var f2 = f1*f2divf1 \n var delta = 1/2*f2*deltaf2\n return f2 + (Math.random() < 0.5 ? -delta : delta)\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":288,"value":"timecourse = {\n run_button;\n var timecourse = Array(50000);\n timecourse[0] = {t:0, N1:N1_0};\n for (var i = 1; i < timecourse.length; i++) {\n timecourse[i] = {t:i, N1:N1_step(timecourse[i-1].N1)};\n };\n return timecourse\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":334,"value":"md `additional (usually costant) variables`","pinned":false,"hidden":false,"mode":"ojs"},{"id":138,"value":"K = 1","pinned":true,"hidden":false,"mode":"ojs"},{"id":100,"value":"f1 = 1.0","pinned":true,"hidden":false,"mode":"ojs"},{"id":542,"value":"md `### analysis`","pinned":false,"hidden":false,"mode":"ojs"},{"id":511,"value":"delta_N1_per_N1 = {\n var data = Array(16);\n for (var i = 0; i < 16; i++) {\n var N1 = i/15\n var f2_high = f1*f2divf1 + 1/2*f1*f2divf1*deltaf2\n var f2_low = f1*f2divf1 - 1/2*f1*f2divf1*deltaf2\n var N1plus_low = (1-d)*N1+(d*K)*(f1*N1)/(f1*N1+f2_low*(K-N1))\n var N1plus_high = (1-d)*N1+(d*K)*(f1*N1)/(f1*N1+f2_high*(K-N1))\n var av_N1plus = (N1plus_low+N1plus_high)/2\n data[i] = {N1:N1, deltaN1:av_N1plus - N1};\n };\n return data\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":813,"value":"md `### simulation 2`","pinned":false,"hidden":false,"mode":"ojs"},{"id":696,"value":"e2_N1_step = (N1) => {\n var f1 = e2_f1()\n return (1-e2_d)*N1+(e2_d*K)*(f1*N1)/(f1*N1+e2_f2()*(K-N1))}","pinned":true,"hidden":false,"mode":"ojs"},{"id":701,"value":"e2_f2 = () => {\n var f2 = f1*e2_f2divf1 \n var delta = 1/2*e2_delta\n return f2 + (Math.random() < 0.5 ? -delta : delta)\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":706,"value":"e2_f1 = () => {\n var delta = 1/2*e2_delta\n return f1 + (Math.random() < 0.5 ? -delta : delta)\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":714,"value":"e2_timecourse = {\n e2_run_button;\n //var epsilon = 10e-9\n var timecourse = Array(50000);\n timecourse[0] = {t:0, N1:e2_N1_0};\n for (var i = 1; i < timecourse.length; i++) {\n var next_N2 = e2_N1_step(timecourse[i-1].N1)\n /*if (next_N2 < epsilon || timecourse[i-1].N1 < 0){\n next_N2 = 0}\n else if (next_N2 > (1-epsilon) || timecourse[i-1].N1 > (1-epsilon)) {\n next_N2 = 1}*/\n timecourse[i] = {t:i, N1:next_N2};\n };\n return timecourse\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":824,"value":"e2_delta_N1_per_N1 = {\n var data = Array(16);\n for (var i = 0; i < 16; i++) {\n var N1 = i/15\n var f2_high = f1*e2_f2divf1 + 1/2*e2_f2divf1*e2_delta\n var f2_low = f1*e2_f2divf1 - 1/2*e2_f2divf1*e2_delta\n \n var f1_high = f1 + 1/2*e2_delta\n var f1_low = f1 - 1/2*e2_delta\n \n var N1plus = Array(4)\n N1plus[0] = (1-d)*N1+(d*K)*(f1_low*N1)/(f1_low*N1+f2_low*(K-N1))\n N1plus[1] = (1-d)*N1+(d*K)*(f1_low*N1)/(f1_low*N1+f2_high*(K-N1))\n N1plus[3] = (1-d)*N1+(d*K)*(f1_high*N1)/(f1_high*N1+f2_low*(K-N1))\n N1plus[4] = (1-d)*N1+(d*K)*(f1_high*N1)/(f1_high*N1+f2_high*(K-N1))\n var av_N1plus = _.sum(N1plus)/4\n data[i] = {N1:N1, deltaN1:av_N1plus - N1};\n };\n return data\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":934,"value":"md `### style`","pinned":false,"hidden":false,"mode":"ojs"},{"id":940,"value":"accent_color = \"#2980b9\"//\"dodgerblue\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":399,"value":"md `### imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":960,"value":"import {debounce} from '@mbostock/debouncing-input'","pinned":true,"hidden":false,"mode":"ojs"},{"id":75,"value":"import {texmd} from \"@airbornemint/texmd\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":109,"value":"import {slider, button} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":213,"value":"import { vl } from \"@vega/vega-lite-api\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":172,"value":"_ = require(\"lodash\");","pinned":false,"hidden":false,"mode":"ojs"},{"id":155,"value":"PD = require('https://bundle.run/probability-distributions@0.9.1')","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}