[R] coupled ODE population model

[Thomas Petzoldt]

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)