miércoles, 12 de enero de 2011

Sistemas Lotka-volterra. Aplicación en R (paquete deSolve)

Sistemas de dinámica poblacional Lotka-volterra


Utilizando el paquete “deSolve” la simulación se reduce a:
#Lotka-Volterra single species discrete population growth.
time <- 20
Ns <- dlogistic(t=20)
plot(0:time, Ns)

#A function to simulate discrete 2 species Lotka-Volterra competition.
alphs <- matrix(c( .01, .005,
.008, .01), ncol=2, byrow=TRUE)
t <- 20
N <- matrix(NA, nrow=t+1, ncol=2)
N[1,] <- c(10,10)
for(i in 1:t) N[i+1,] <- dlvcomp2(N[i,], alphs)
matplot(0:t, N, type='l', col=1, ylim=c(0,110))
abline(h=1/alphs[1,1], lty=3)
text(0, 1/alphs[1,1], "K", adj=c(0,0))
legend("right", c(expression("Sp.1 "*(alpha[21]==0.008)),
expression("Sp.2 "*(alpha[12]==0.005))), lty=1:2, bty='n')

#A Lotka-Volterra model of intraguild predation, after Holt and Polis (1997). For use with ode in the deSolve package.
library(deSolve)
params <- c(bpb= 0.032, abp=10^-8, bpn=10^-5, anp=10^-4, mp=1,
bnb=0.04, abn=10^-8, mn=1,
r=1, abb=10^-9.5)
t=seq(0, 60, by=.1)
N.init <- c(B = 10^9, N = 10^4, P = 10^3)
igp.out <- ode(N.init, t, igp, params)
matplot(t, log10(igp.out[,-1]+1), type="l",ylab="log(Abundance)")
legend('right', c("B", "N", "P"), lty=1:3, col=1:3, lwd=2,bty="n")

#System of ordinary differential equations for two species Lotka-Volterra competition. For use with ode in the deSolve package.
function (t, n, parms)
{
with(as.list(parms), { dn1dt <- r1 * n[1] * (1 - a11 * n[1] - a12 * n[2]); dn2dt <- r2 * n[2] * (1 - a22 * n[2] - a21 * n[1]); list(c(dn1dt, dn2dt)) })
}
library(deSolve)
parms <- c(r1 = 1, r2 = 0.1, a11 = 0.2, a21 = 0.1, a22 = 0.02, a12 = 0.01)
initialN <- c(1, 1)
out <- ode(y = initialN, times = 1:100, func = lvcomp2, parms = parms)
matplot(out[, 1], out[, -1], type = "l")

#System of ordinary differential equations for three species Lotka-Volterra competition. For use with ode in the deSolve package.
function (t, n, parms)
{
with(as.list(parms), { dn1dt <- r1 * n[1] * (1 - a11 * n[1] - a12 * n[2] - a13 *; n[3]); dn2dt <- r2 * n[2] * (1 - a22 * n[2] - a21 * n[1] - a23 * n[3]); dn3dt <- r3 * n[3] * (1 - a33 * n[3] - a31 * n[1] - a32 * n[2]); list(c(dn1dt, dn2dt, dn3dt)) })
}
library(deSolve)
parms <- c(r1 = 0.1, r2 = 0.2, r3 = 0.3,
a11 = 0.1, a12 = 0.01, a13 = 0.01,
a21 = 0.01, a22 = 0.15, a23 = 0.01,
a31 = 0.01, a32 = 0.01, a33 = 0.2)
initialN <- c(1, 1, 1)
out <- ode(y = initialN, times = 1:100, func = lvcomp3, parms = parms)
matplot(out[, 1], out[, -1], type = "l")

#A general Lotka-Volterra competition model, for any number of species. For use with ode in the deSolve package. This function uses a vector and matrix within the list of parameters.
S <- 10
alpha <- .01
r <- runif(S)*2
a <- matrix(rnorm(S^2, m=alpha, sd=alpha/10), nrow=S, ncol=S)
parms <- list(r,a)
t=seq(0,40, by=.1)
N0 <- runif(S)/(S*alpha)
library(deSolve)
lvout <- ode(N0, t, lvcompg, parms)
matplot(t, lvout[,-1], type="l", ylab="N", log='y')

#Analysis of Jacobian Lotka-Volterra Food Web Matrices: Used primarily to repeat simulations and analyses of Pimm and Lawton (1977), given a Jacobian matrix. Analyses include eigenanalysis, but also measuring average interaction strength (May 1972), average intraspecific negative density dependence, and the strength of the omnivory interaction, if present.
Aq = matrix(c(
-1, -10, 0, 0,
0.1, 0, -10, 0,
0, 0.1, 0, -10,
0, 0, 0.1, 0),
nrow=4, byrow=TRUE)
pimmlawton(Aq, N=1)
out <- pimmlawton(Aq, N=2000)
out <- subset(out, -1/DomEig <150)
hist(-1/out$DomEig, main="Frequencies of Return Time")

Otros modelos interesantes de interacción inter-específica
#The Lotka-Volterra predator-prey model, for use with ode in the deSolve package.
params1 <- c(b=.5, a=.01, s=.2, e=.1)
Time <- seq(0,100, by=.1) # Set time here
LV.out <- ode(c(H0=25,P0=5), Time, predpreyLV, params1)
matplot(Time, (LV.out[,2:3]), type="l", ylab="Population Size")