Skip to content

coupled ODE population model

3 messages · Justin Frank, Thomas Petzoldt, Spencer Graves

Original
Justin Frank
# 
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
Thomas Petzoldt
#
Justin Frank wrote:
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)
Spencer Graves
# 
Hi, Justin: 


      In addition to the example below, Thomas Petzoldt also has two 
vignettes included with his "simecol" package.  These might help you 
more easily use his package. 


      Have you looked for what you want using the "RSiteSearch.function" 
in the "RSiteSearch" package?  Consider, for example, the following: 


library(RSiteSearch)
popec <- RSiteSearch.function('population ecology')
hits(popec)  #  100 help pages
summary(popec) #   29 packages
HTML(popec) # displays a table in a web browser


      This finds 100 different help pages in 29 different packages 
including this search term.  The last line opens a web browser with a 
table of the results sorted to put the package with the most hits first 
and providing links to HTML versions of those 100 help pages.  The 
RSiteSearch package also includes the capability to combine different 
searches using union and intersection. 


      Hope this helps. 
      Spencer
Thomas Petzoldt wrote: