Victor committed Oct 11, 2020 1 2 3 4 5 # A first model of meta population In this script, we model the evolution of a population where agents are simply defined by their position on some landscape. We implement the simplest possible birth and death function. ## The landscape Victor committed Oct 12, 2020 6 7 8 Let's start by a linear landscape. We define a discrete segment of length 9, with reflecting boundary conditions. In fact, reflecting boundary conditions are implemented for any finite space. !!! warning "1D Reflections" Victor committed Oct 12, 2020 9 As of v1, only reflections on one dimensional space are implemented. We have a two dimensional reflection method that will be released in the future. Victor committed Oct 12, 2020 10 11 12 13 14 15 16 17 There are two ways of implementing a linear landscape. The first one uses a DiscreteSegment while the second relies on LightGraphs.jl library. Both are *almost* equivalent. ### DiscreteSegment julia using ABMEv nodes = 10 mysegment = DiscreteSegment(1,nodes) Victor committed Oct 12, 2020 18 wholespace = (mysegment,) Victor committed Oct 12, 2020 19 20 21 22 23 24 25 26  ### grid julia using ABMEv, LightGraphs nodes = 10 g = grid([nodes,1]) mysegmentgraph = GraphSpace(g) Victor committed Oct 12, 2020 27 wholespace = (mysegmentgraph,) Victor committed Oct 12, 2020 28  Victor committed Oct 12, 2020 29 30 !!! warning "Space tuple" Notice that the whole space should be a *tuple of spaces*. Even where there is only one sub vector space as here, you need to have brackets and comma around the unit vector space. Victor committed Oct 12, 2020 31 32 33 Here is how you can visualise the landscape. julia Victor committed Oct 12, 2020 34 35 using GraphPlot gplot(g, collect(1:nodes), collect(1:nodes)) Victor committed Oct 12, 2020 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55  ## Defining competition processes We propose that any individual have a constant birth rate, and competes with all the individuals present in the same patch. Let i \in \N,x_i \in \{1,2,\dots,9\}. The competition pressure experience by individual i is such that math d(x_i) = \sum_j \delta(x_i-x_j)  where \delta is the dirac function. In this way, we recover a logistic growth function for subpopulation within a patch. julia K0 = 1000 # We will have in total 1000 individuals b(X,t) = 1 / nodes d(X,Y,t) = (X[1] ≈ Y[1]) / K0  At equilibrium, population size in each deme will reach K0 / nodes. Victor committed Oct 12, 2020 56 57 58 !!! warning "Time dependency" Even though time is not used, you have to write birth and death functions with time dependency. Victor committed Oct 12, 2020 59 60 61 62 63 64 65 66 67 ## Dispersal We assume that anytime an offspring is born, it is given a chance to move (\mu = 1). julia mu = [1.] D = (1.5,)  ## Running the world We initialise the world with initial population of size K_0 / 9 located on patch 5. We keep track of individuals' ancestors by setting ancestors=true. Because we wish to use Gillepsie algorithm, we need rates=true as agents' internal birth and death rates are updated at every time step. !!! note "Warning" Victor committed Oct 12, 2020 68 rates treatment is something we might implement in the library internals. Victor committed Oct 12, 2020 69 70 julia Victor committed Oct 12, 2020 71 using UnPack# useful macro @pack! Victor committed Oct 12, 2020 72 NMax = 2000 Victor committed Oct 12, 2020 73 tend = 300. Victor committed Oct 12, 2020 74 75 76 77 78 p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax myagents = [Agent(myspace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes] w0 = World(myagents,myspace,p,0.) @time sim = run!(w0,Gillepsie(),tend)  Victor committed Oct 12, 2020 79 80 This is the simplest run you can do. Now time to more interesting things Victor committed Oct 12, 2020 81 ## Analysis Victor committed Oct 12, 2020 82 83 ### Size of the world Victor committed Oct 12, 2020 84 Let's verify that the population's growth is logistic. We will plot the population size over time. Victor committed Oct 12, 2020 85 To do so, one need to define dt_saving < tend to save every dt_saving time steps of the world. Victor committed Oct 12, 2020 86 87 julia Victor committed Oct 12, 2020 88 89 90 91 myagents = [Agent(wholespace,(5,),ancestors=true,rates=true) for i in 1:K0/nodes] w0 = World(myagents,wholespace,p,0.) # we need to reinitialise the world @time sim = run!(w0,Gillepsie(),tend,dt_saving=2.) wsize = [length(w) for w in sim[:]] Victor committed Oct 12, 2020 92 using Plots Victor committed Oct 12, 2020 93 94 95 96 97 Plots.plot(get_tspan(sim),wsize, label = "", ylabel = "Metapopulation size", xlabel ="time", grid = false) Victor committed Oct 12, 2020 98  Victor committed Oct 12, 2020 99 100 ![](tutorials/delta_comp_wsize.png) Victor committed Oct 12, 2020 101 Victor committed Oct 12, 2020 102 103 !!! notes "Callbacks" Note that one could also use a callback function to obtain time series of size of the world computed at simulation time. See [Simulation page](../manual/simulation.md). Victor committed Oct 12, 2020 104 Victor committed Oct 12, 2020 105 106 107 108 109 110 111 112 113 114 115 116 ### Position through time One might be tempted to plot the agents position for some time steps. julia Plots.plot(sim, label = "", ylabel = "Geographical position", grid = false, markersize = 10)  ![](tutorials/delta_comp_pos.png)