Assessment - BIOLM0029 - 1719614
Create three species Lotka-Volterra model
Create cluster, register and load packages to it:
library(parallel)
library(doSNOW)
cl <- makeCluster(detectCores() -1) #Make the cluster which uses the computer's number of cores -1.
registerDoSNOW(cl) #Register the cluster.
clusterEvalQ(cl,{ #Load the relevant packages to the cluster.
library(parallel)
library(base)
library(tidyverse)
library(data.table)
library(foreach)
library(deSolve)
}) Define three species model, based on the equations:
## Define the three species LV model:
lv3 <- function(t, start, parms) {
##Setting the start values for each of the species:
## basal prey species
R <- start["R"]
## primary consumers
N <- start["N"]
## top predators
P <- start["P"]
##Allowing R to look within parms for the values of r, a, e, etc
with(as.list(parms), {
## dynamics of the basal prey species #R grows at rate r.
dR <- r*R - a*N*R #a - predation rate of N on R.
## dynamics of the primary consumers #f - conversion efficiency of R to new N.
dN <- f*a*N*R - q*N - e*N*P #q - exponential decline rate of N in absence of R.
##dynamics of the top predator #e - predation of P on N.
dP <- z*e*N*P - b*P #z - conversion efficiency of N to new P.
#b - exponential decline rate of P in absence of N.
##return a list of the abundances of each species
##deSolve will then group these into a data frame
list(c(dR, dN, dP)) #return all 3 species
})
}
#By passing these two equations to the ODE solver (Ordinary Differential Equation solver), we tell the simulation to set all abundances to zero once one species reaches an abundance of zero.
eventFun <- function(t,y,p){
y <- c(0, 0, 0)
return (y)
}
rootfun <- function (t,y,p) {
if(min(y)<1){
y <- c(0, 0, 0)
}
return(y)
}Create ODE (Ordinary Differential Equation) solver function for simulations:
#A wrapper function containing the ODE solver and lines for reshaping the data outputs, as well as a section to inlcude a list of the simulation results and the parameters used.
lv3_sim <- function(start, parms, time){
##run the simulation
library(tidyverse)
sim <- as_tibble(lsoda(y = start,
times = time,
func = lv3,
parms = parms,
##pass the event function
events = list(func = eventFun,
##yes we want to use a root
##function
root = TRUE,
##tell lsoda to stop the
##simulation if the event
##is triggered
terminalroot = 1),
##and the root function
rootfun = rootfun))
##reshape the data to long format:
longer_sim <- sim %>% pivot_longer(-time,
values_to = "abundance",
names_to="species")
##number of years the time series ran for:
longer_sim$sim_length<-max(longer_sim$time)
##add a column to allow us to easily split the results up later
longer_sim$parms <- paste(names(parms), parms, sep="=", collapse = ",")
##make a list of the simulation results and the
##parameters which made them
res <- list("sim" = longer_sim,
"parameters" = parms)
return(res)
}Specify parameters used simulations:
##Make an object of the model parameters to be passed to the ODE solver:
parms <- c(r = 0.1, #growth rate of resources
a = 0.01, #feeding rate of primary consumer
#on basal species
f = 0.01, #efficiency of feeding on resources
q = 0.1, #intrinsic death rate of the
#primary consumer
e = 0.01, #feeding rate of top predator
#on primary consumer
z = 0.1, #efficiency of predating on
#primary consumer
b = 0.1) #intrinsic death rate of top predator
#Define the starting population parameters:
start <- c(R = 1125, # basal species is most numerous
N = 225, # still many primary consumers
P = 20) #fewest number of top predators
##Set the length of the simulation (100 time steps), and the resolution (set size, in this case 0.1 time steps)
time <- seq(0, 100, 1)Run and plot initial simulation using the specified parameters:
dd <- lv3_sim(start, parms, time) #Passes the specified parameters to the ODE solver, saves result in dd#Plots the returned result from dd
initial <- ggplot(dd$sim, aes(x=time, y=abundance, col=species)) +
geom_line() +
scale_y_log10() +
theme_minimal() +
ggtitle("Three interacting species with defined starting values and parameters")
initial## Don't know how to automatically pick scale for object of type deSolve/matrix. Defaulting to continuous.
Three species interacting through time for the static specified parameter and starting values:
## r a f q e z b
## 0.10 0.01 0.01 0.10 0.01 0.10 0.10Simulate multiple model outputs
We can run multiple simulations in parallel. We create an object containing all possible combinations of the parameters, shifted to range between 0.01 and 0.2. Length.out specifies how many values in the sequence between 0.01 and 0.2, i.e. how many simulation replicates we run of different possible parameter combinations.
#make a list where the sequence 0.01 to 0.2 (with a specified length of the
#sequence) is replicated once for each of the parameters in parms
parms_list <- rep(list(seq(0.01, 0.2, length.out = 4)),
length(parms))
##use expand.grid() to make a data.frame of every possible
##combination of the parameters
all_var <- expand.grid(parms_list) #expanding this creates a data.table of the different possible combinations of the parameters, and is the width of the parameters from parms, whilst multiplied by the length.out value to the power of the number of parms, i.e. all possible combinations
##look at the data
head(all_var)## Var1 Var2 Var3 Var4 Var5 Var6 Var7
## 1 0.01000000 0.01000000 0.01 0.01 0.01 0.01 0.01
## 2 0.07333333 0.01000000 0.01 0.01 0.01 0.01 0.01
## 3 0.13666667 0.01000000 0.01 0.01 0.01 0.01 0.01
## 4 0.20000000 0.01000000 0.01 0.01 0.01 0.01 0.01
## 5 0.01000000 0.07333333 0.01 0.01 0.01 0.01 0.01
## 6 0.07333333 0.07333333 0.01 0.01 0.01 0.01 0.01## the dimensions of the data
dim(all_var)## [1] 16384 7##convert that into a list, where each object is 1 row of the data frame
all_var_list <- as.list(as.data.frame(t(all_var)))
head(all_var_list)## $V1
## [1] 0.01 0.01 0.01 0.01 0.01 0.01 0.01
##
## $V2
## [1] 0.07333333 0.01000000 0.01000000 0.01000000 0.01000000 0.01000000 0.01000000
##
## $V3
## [1] 0.1366667 0.0100000 0.0100000 0.0100000 0.0100000 0.0100000 0.0100000
##
## $V4
## [1] 0.20 0.01 0.01 0.01 0.01 0.01 0.01
##
## $V5
## [1] 0.01000000 0.07333333 0.01000000 0.01000000 0.01000000 0.01000000 0.01000000
##
## $V6
## [1] 0.07333333 0.07333333 0.01000000 0.01000000 0.01000000 0.01000000 0.01000000
The length.out determines the number of simulations run: e.g. length.out=1 -> 7 simulations, length.out=3-> 2187, length.out=5 -> 78125 etc. Therefore, with our parameters we see that by running more simulations, some reach longer times without any species becoming extinct. With length.out=4, we get 16,384 simulations to analyse.
Create progress bar, run parallel simulations on cluster:
##make a progress bar, we need to set the maximum to the
##maximum length of the sequence we are iterating across
pb <- txtProgressBar(max = max(length(all_var_list)), style = 3)##
|
| | 0%##make a function to pass the progress bar to the clusters
progress <- function(n) setTxtProgressBar(pb, n)
##and set the options for snow to include a progress bar
opts <- list(progress = progress)
#foreach already loaded into cluster
make_parallel <- foreach(x = all_var_list, #run the simulations of all the possible combinations for the parameters
.options.snow = opts) %dopar% { #use a progress bar
#add names to parameters
names(x) <- names(parms)
dd <- lv3_sim(start = start, parms = x, time = time) #use the paramaters from all_var_list, and pass each of them through the simulation
return(dd)
}##
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%Extract and plot results of all simulations:
#data.table and parallel loaded into cluster already
extract_res <- rbindlist(lapply(make_parallel, function(h){
##extract and return the simulation values
return(h$sim)
}))
##plot the results from the parallel simulation
plot1 <- ggplot(extract_res, aes(x = time,
y = abundance,
group = parms)) +
geom_line(col = "black", alpha=0.3) +
facet_wrap(~species) +
##log the y axis
scale_y_log10() +
##dark theme
theme_light() +
ggtitle("All simulations")
plot1## Don't know how to automatically pick scale for object of type deSolve/matrix. Defaulting to continuous.
By looking at simulations of all parameter combinations, we see the majority fail with a species becoming extinct early on. We can filter simulations to see how many surpass time>=20.
plot2 <- ggplot(filter(extract_res, sim_length >= 20),
aes(x = time,
y = abundance,
group = parms)) +
geom_line(col = "black", alpha=0.3) +
facet_wrap(~species) +
##log the y axis
scale_y_log10() +
##dark theme
theme_light() +
ggtitle("Simulations surpassing time >= 20")
plot2## Don't know how to automatically pick scale for object of type deSolve/matrix. Defaulting to continuous.
b <- filter(extract_res, sim_length >=20) #filter for simulations which last longer than 20
sum_b <- b%>% group_by(parms) %>% summarise() #group by parms, how many unique sim paths are there## `summarise()` ungrouping output (override with `.groups` argument)sum_b #print the unique simulations## # A tibble: 58 x 1
## parms
## <chr>
## 1 r=0.01,a=0.01,f=0.01,q=0.0733333333333333,e=0.01,z=0.2,b=0.136666666666667
## 2 r=0.01,a=0.01,f=0.01,q=0.136666666666667,e=0.01,z=0.136666666666667,b=0.1366~
## 3 r=0.01,a=0.01,f=0.01,q=0.136666666666667,e=0.01,z=0.2,b=0.2
## 4 r=0.01,a=0.01,f=0.01,q=0.2,e=0.01,z=0.01,b=0.136666666666667
## 5 r=0.0733333333333333,a=0.01,f=0.01,q=0.01,e=0.01,z=0.2,b=0.136666666666667
## 6 r=0.0733333333333333,a=0.01,f=0.01,q=0.0733333333333333,e=0.01,z=0.136666666~
## 7 r=0.0733333333333333,a=0.01,f=0.01,q=0.0733333333333333,e=0.01,z=0.2,b=0.136~
## 8 r=0.0733333333333333,a=0.01,f=0.01,q=0.0733333333333333,e=0.01,z=0.2,b=0.2
## 9 r=0.0733333333333333,a=0.01,f=0.01,q=0.136666666666667,e=0.01,z=0.01,b=0.073~
## 10 r=0.0733333333333333,a=0.01,f=0.01,q=0.136666666666667,e=0.01,z=0.0733333333~
## # ... with 48 more rows58 simulations surpass time>=20.
There is still lots of data, so we filter for time>=30:
plot3 <- ggplot(filter(extract_res, sim_length >= 30),
aes(x = time,
y = abundance,
group = parms)) +
geom_line(col = "black", alpha=0.3) +
facet_wrap(~species) +
##log the y axis
scale_y_log10() +
##dark theme
theme_light() +
ggtitle("Simulations surpassing time >= 30")
plot3## Don't know how to automatically pick scale for object of type deSolve/matrix. Defaulting to continuous.
c <- filter(extract_res, sim_length >= 30) #filter for simulations which last longer than 30
sum_c <- c %>% group_by(parms) %>% summarise() #group by parms, how many unique sim paths are there## `summarise()` ungrouping output (override with `.groups` argument)sum_c## # A tibble: 8 x 1
## parms
## <chr>
## 1 r=0.0733333333333333,a=0.01,f=0.01,q=0.01,e=0.01,z=0.2,b=0.136666666666667
## 2 r=0.0733333333333333,a=0.01,f=0.01,q=0.0733333333333333,e=0.01,z=0.1366666666~
## 3 r=0.136666666666667,a=0.01,f=0.01,q=0.01,e=0.01,z=0.2,b=0.136666666666667
## 4 r=0.136666666666667,a=0.01,f=0.01,q=0.0733333333333333,e=0.01,z=0.13666666666~
## 5 r=0.2,a=0.01,f=0.01,q=0.01,e=0.01,z=0.136666666666667,b=0.136666666666667
## 6 r=0.2,a=0.01,f=0.01,q=0.01,e=0.01,z=0.2,b=0.136666666666667
## 7 r=0.2,a=0.01,f=0.01,q=0.0733333333333333,e=0.01,z=0.01,b=0.0733333333333333
## 8 r=0.2,a=0.01,f=0.01,q=0.0733333333333333,e=0.01,z=0.136666666666667,b=0.13666~Where the starting population values and range of specified parameters influence the species interactions, we see very few simulations surpassing time>=30.
We can further reduce simulations in order to compare the parameters:
plot4 <- ggplot(filter(extract_res, sim_length >= 33),
aes(x = time,
y = abundance,
group = parms)) +
geom_line(col = "black", alpha=0.3) +
facet_wrap(~species) +
##log the y axis
scale_y_log10() +
##dark theme
theme_light() +
ggtitle("Simulations surpassing time >= 33")
plot4## Don't know how to automatically pick scale for object of type deSolve/matrix. Defaulting to continuous.
d <- filter(extract_res, sim_length >= 33) #filter for simulations which last longer than 33
sum_d <- d %>% group_by(parms) %>% summarise() #group by parms, how many unique sim paths are there## `summarise()` ungrouping output (override with `.groups` argument)sum_d ## # A tibble: 4 x 1
## parms
## <chr>
## 1 r=0.0733333333333333,a=0.01,f=0.01,q=0.01,e=0.01,z=0.2,b=0.136666666666667
## 2 r=0.136666666666667,a=0.01,f=0.01,q=0.01,e=0.01,z=0.2,b=0.136666666666667
## 3 r=0.2,a=0.01,f=0.01,q=0.01,e=0.01,z=0.2,b=0.136666666666667
## 4 r=0.2,a=0.01,f=0.01,q=0.0733333333333333,e=0.01,z=0.01,b=0.0733333333333333Although we see a pattern in these parameters, we can check it using another set of values.
Other starting values
Having previously run simulations with another set of starting values (where simulations surpass time>=15) using the same parameter and length.out values, we can compare graph outputs.
After filtering for the longest simulations, there are three which surpass time>=25. Their parameters are:
By looking at the parameter values for these longest running simulations, with different starting values, we see that almost all rely on r (basal species (R) growth rate) being at least 0.073, alongside the values z (the conversion efficiency of N to P) and b (top predator (P) exponential decline rate) being similarly high, at 0.136 or higher. From this we can predict the model’s stability depends on these three values being relatively high for longer running simulations.
The other long-running simulation has a low z (0.01) and relatively low b (0.073) value, alongside a higher q value (primary consumer (N) exponential decline rate) at 0.073. Therefore, where N is reduced by this higher q, it’s population is stabilized by the high r value and low r & z values, keeping the three species relatively stable.
The starting values and number of replications have obvious contributions to the modelling of these simulations, however to achieve relatively stable long-running simulations, the most impactful parameters are r, z and b.