In keepsimpler/perturbEcol: Simulation of ecological dynamic systems under environmental perturbations

library(perturbEcol)
library(deSolve)
library(rootSolve)
library(plyr)
library(igraph)

library(doMC)  # 
registerDoMC(1)  # register Multi Cores
getDoParWorkers()  # get available Cores

extinct_threshold <- .Machine$double.eps * 100  # threshold of species become extinct
s = 10  # number of species

# Initialize an object of camModel
cam <- new("camModel")
cam@times <- 1:200  # time steps
cam@perturb <- perturb_growthrate  # perturbation function
cam@perturbNum <- 500  # number of iterated steps
cam@solver <- sim_ode_press # 

First, we evaluate the special case of all interactions are mutualistic.

k = 3
hybrid_graph = gen_hybrid_network(s, k, type = 'er', pc = 0., pa = 0., pm = 1.0)
coeff <- c(alpha.mu = 0.2, alpha.sd = 0.1, beta0.mu = 1., beta0.sd = 0.1, beta1.mu = 0.0, beta1.sd = 0.0, antago.mu = 0., antago.sd = 0.0, gamma.mu = 0.1, gamma.sd = 0., h.mu = 0.1, h.sd = 0., g.mu = 0., g.sd = 0., delta = 0.)
cam@params <- params_lv2_cam(hybrid_graph, coeff)
cam@init <- init_lv2_cam(cam@params)
cam <- sim(cam, extinct_threshold = extinct_threshold, r.delta.mu = 0.005, r.delta.sd = 0.00)
plot(cam, main = 'hybrid, full graph')

gamma.mus = seq(from = 0.1, to = 0.28, by = 0.05)
gamma.mus = rep(gamma.mus, times = 10)
ret <- lapply(gamma.mus, function(gamma.mu) {
  coeff['gamma.mu'] = gamma.mu
  cam@params <- params_lv2_cam(hybrid_graph, coeff)
  cam@init <- init_lv2_cam(cam@params)
  cam <- sim(cam, extinct_threshold = extinct_threshold, r.delta.mu = 0.005, r.delta.sd = 0.00)
  list(fragility = cam@out$fragility, coeff = coeff)
})

ret2 <- sapply(ret, function(one) {
  c(gamma.mu = one$coeff[['gamma.mu']], fragility.variance = one$fragility$variance, fragility.entropy = one$fragility$entropy, resistance = one$fragility$resistance)
})
plot(ret2['gamma.mu', ], ret2[c('fragility.variance'), ])
plot(ret2['gamma.mu', ], ret2[c('fragility.entropy'), ])
plot(ret2['gamma.mu', ], ret2[c('resistance'), ])

Next, we add antagonistic interactions

k = 3
hybrid_graph = gen_hybrid_network(s, k, type = 'er', pc = 0., pa = 0.2, pm = 0.8)
coeff <- c(alpha.mu = 0.2, alpha.sd = 0.1, beta0.mu = 1., beta0.sd = 0.1, beta1.mu = 0.0, beta1.sd = 0.0, antago.mu = 0.1, antago.sd = 0.0, gamma.mu = 0.1, gamma.sd = 0., h.mu = 0.1, h.sd = 0., g.mu = 1., g.sd = 0., delta = 0.)
cam@params <- params_lv2_cam(hybrid_graph, coeff)
cam@init <- init_lv2_cam(cam@params)
cam <- sim(cam, extinct_threshold = extinct_threshold, r.delta.mu = 0.01, r.delta.sd = 0.00)
plot(cam, main = 'hybrid, full graph')

gamma.mus = seq(from = 0.1, to = 0.5, by = 0.05)
gamma.mus = rep(gamma.mus, times = 10)
ret <- llply(gamma.mus, .parallel = T, function(gamma.mu) {
  coeff['gamma.mu'] = gamma.mu
  coeff['antago.mu'] = gamma.mu
  cam@params <- params_lv2_cam(hybrid_graph, coeff)
  cam@init <- init_lv2_cam(cam@params)
  cam <- sim(cam, extinct_threshold = extinct_threshold, r.delta.mu = 0.01, r.delta.sd = 0.00)
  list(fragility = cam@out$fragility, coeff = coeff)
})

ret2 <- sapply(ret, function(one) {
  c(gamma.mu = one$coeff[['gamma.mu']], fragility.variance = one$fragility$variance, fragility.entropy = one$fragility$entropy, resistance = one$fragility$resistance)
})
plot(ret2['gamma.mu', ], ret2[c('fragility.variance'), ])
plot(ret2['gamma.mu', ], ret2[c('fragility.entropy'), ])
plot(ret2['gamma.mu', ], ret2[c('resistance'), ])

Next, we add antagonistic interactions

k = 3
hybrid_graph = gen_hybrid_network(s, k, type = 'er', pc = 0., pa = 0.5, pm = 0.5)


coeff <- c(alpha.mu = 0.2, alpha.sd = 0.1, beta0.mu = 1., beta0.sd = 0.1, beta1.mu = 0.0, beta1.sd = 0.0, antago.mu = 0.1, antago.sd = 0.0, gamma.mu = 0.1, gamma.sd = 0., h.mu = 0.1, h.sd = 0., g.mu = 1., g.sd = 0., delta = 0.)
cam@params <- params_lv2_cam(hybrid_graph, coeff)
cam@init <- init_lv2_cam(cam@params)
cam <- sim(cam, extinct_threshold = extinct_threshold, r.delta.mu = 0.01, r.delta.sd = 0.00)
plot(cam, main = 'hybrid, full graph')

gamma.mus = seq(from = 0.1, to = 0.5, by = 0.05)
gamma.mus = rep(gamma.mus, times = 10)
ret <- lapply(gamma.mus, function(gamma.mu) {
  coeff['gamma.mu'] = gamma.mu
  coeff['antago.mu'] = gamma.mu
  cam@params <- params_lv2_cam(hybrid_graph, coeff)
  cam@init <- init_lv2_cam(cam@params)
  cam <- sim(cam, extinct_threshold = extinct_threshold, r.delta.mu = 0.01, r.delta.sd = 0.00)
  list(fragility = cam@out$fragility, coeff = coeff)
})

ret2 <- sapply(ret, function(one) {
  c(gamma.mu = one$coeff[['gamma.mu']], fragility.variance = one$fragility$variance, fragility.entropy = one$fragility$entropy, resistance = one$fragility$resistance)
})
plot(ret2['gamma.mu', ], ret2[c('fragility.variance'), ])
plot(ret2['gamma.mu', ], ret2[c('fragility.entropy'), ])
plot(ret2['gamma.mu', ], ret2[c('resistance'), ])

Next, we add antagonistic interactions

k = 3
hybrid_graph = gen_hybrid_network(s, k, type = 'er', pc = 0., pa = 1., pm = 0.)
coeff <- c(alpha.mu = 0.2, alpha.sd = 0.1, beta0.mu = 1., beta0.sd = 0.1, beta1.mu = 0.0, beta1.sd = 0.0, antago.mu = 0.1, antago.sd = 0.0, gamma.mu = 0.1, gamma.sd = 0., h.mu = 0.1, h.sd = 0., g.mu = 1., g.sd = 0., delta = 0.)
cam@params <- params_lv2_cam(hybrid_graph, coeff)
cam@init <- init_lv2_cam(cam@params)
cam <- sim(cam, extinct_threshold = extinct_threshold, r.delta.mu = 0.01, r.delta.sd = 0.00)
plot(cam, main = 'hybrid, full graph')

gamma.mus = seq(from = 0.1, to = 0.5, by = 0.05)
gamma.mus = rep(gamma.mus, times = 10)
ret <- lapply(gamma.mus, function(gamma.mu) {
  coeff['gamma.mu'] = gamma.mu
  coeff['antago.mu'] = gamma.mu
  cam@params <- params_lv2_cam(hybrid_graph, coeff)
  cam@init <- init_lv2_cam(cam@params)
  cam <- sim(cam, extinct_threshold = extinct_threshold, r.delta.mu = 0.01, r.delta.sd = 0.00)
  list(fragility = cam@out$fragility, coeff = coeff)
})

ret2 <- sapply(ret, function(one) {
  c(gamma.mu = one$coeff[['gamma.mu']], fragility.variance = one$fragility$variance, fragility.entropy = one$fragility$entropy, resistance = one$fragility$resistance)
})
plot(ret2['gamma.mu', ], ret2[c('fragility.variance'), ])
plot(ret2['gamma.mu', ], ret2[c('fragility.entropy'), ])
plot(ret2['gamma.mu', ], ret2[c('resistance'), ])

keepsimpler/perturbEcol documentation built on May 20, 2019, 8:45 a.m.