Justin Frank wrote:>
> I'm fairly new to R, and I'm trying to write out a population model
that
> satisfies the following;
>
> the system consists of s species, i= 1, 2,...,s
> network of interactions between species is specified by a (s x s) real
matrix,
> C[i,j]
>
> x[i] being the relative population of the "ith" species (0 =<
x[i] =< 1,
> sum(x[i]=1)
>
> the evolution rule being considered is as follows;
>
> xprime[i] = f[i] if x[i] > 0 or f[i] >= 0
>
> xprime[i] = 0 if x[i] = 0 and f[i] < 0
>
> where f[i] = sum(C[i,j]*x[j]) - x[i]*sum(C[k,j]*x[j])
>
>
> I have a bit of attempted code written out, but are there any tricks or
tips
> that would condense or make this mess look nicer?
>
> -Justin
Hi Justin,
you may consider using package simecol (http://www.simecol.de). You find
several examples on the help pages and also some PDF files (so called
vignettes). If you have many species, you should consider using matrix
formulations. The example below demonstrates a multi-species
Lotka-Volterra model with an interaction matrix A.
While the example is not yet your model it may serve as a starting
point. If you have any further questions, please let me know.
Thomas Petzoldt
library(simecol)
##########################################################
## Basic Multi Species Predator Prey Model
## 1) parameters for pairwise single species interaction
##########################################################
LVPP <- new("odeModel",
main = function(t, n, parms) {
with(parms, {
dn <- r * n + n * (A %*% n)
return(list(c(dn)))
})
},
parms = list(
r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
## only pairwise interactions:
A = matrix(c(0.0, 0.0, -0.2, 0.0, # prey 1
0.0, 0.0, 0.0, -0.1, # prey 2
0.2, 0.0, 0.0, 0.0, # predator 1; eats prey 1
0.0, 0.1, 0.0, 0.0), # predator 2; eats prey 2
nrow = 4, ncol = 4, byrow=TRUE)
),
times = seq(from=0, to=500, by = 0.1),
init = c(prey1=1, prey2=1, pred1=2, pred2=2),
solver = "lsoda"
)
plot(sim(LVPP))
##########################################################
## 2) multi species interactions
##########################################################
## make two clones of LVPP
LVPPweak <- LVPPstrong <- LVPP
## a helper function
## this copies the negative of the upper triangular to the lower
makeSym <- function(A, f = -1) {
ind <- lower.tri(A)
A[ind] <- f * t(A)[ind]
A
}
## weak coupling
A <- matrix(c(0.0, 0.0, -0.1, -0.001, # prey 1
NA, 0.0, -0.002, -0.2, # prey 2
NA, NA, 0.0, 0.0, # predator 1
NA, NA, NA, 0.0), # predator 2
nrow = 4, ncol = 4, byrow=TRUE)
parms(LVPPweak)$A <- makeSym(A)
## stronger coupling
A <- matrix(c(0.0, 0.0, -0.1, -0.05, # prey 1
NA, 0.0, -0.02, -0.2, # prey 2
NA, NA, 0.0, 0.0, # predator 1
NA, NA, NA, 0.0), # predator 2
nrow = 4, ncol = 4, byrow=TRUE)
parms(LVPPstrong)$A <- makeSym(A)
LVPPweak <- sim(LVPPweak)
LVPPstrong <- sim(LVPPstrong)
plot(LVPPweak)
plot(LVPPstrong)
o <- out(LVPPweak)
par(mfrow=c(2,2))
plot(o$prey1, o$pred1)
plot(o$prey2, o$pred2)
plot(o$prey1, o$pred2)
plot(o$prey2, o$pred1)
o <- out(LVPPstrong)
par(mfrow=c(2,2))
plot(o$prey1, o$pred1)
plot(o$prey2, o$pred2)
plot(o$prey1, o$pred2)
plot(o$prey2, o$pred1)