lv_interaction {gauseR}

Lotka-Volterra Interactions

Description

Calculates dn/dt for n species in a Lokta-Volterra system, following the form: dni/dt = ni * (ri + aii * ni + sum_j(aij * nj)) Note that aii coefficients can be positive or negative, although positive coefficients risk having the system run to infinite population sizes, which will crash the function.

Usage

lv_interaction(time, n, parms)

Arguments

time

The time steps corresponding to each observation - exists to interface with ode function, but should be left blank.

n

A vector of species abundances

parms

A vector of parameters - the first n elements should be the growth rates r1, r2, ... rn for all n species. The remaining terms should be the elements of the interaction matrix A, listed in the order a11, a12, ... a1n, a21, a22, ... a2n, ... an1, an2, ... ann.

Value

vector of growth rates for each species

Examples

# load data from competition experiment
data(gause_1934_science_f02_03)

# subset data to include just mixtures
mixturedata<-gause_1934_science_f02_03[gause_1934_science_f02_03$Treatment=="Mixture",]

# get time-lagged observations for each species
Pc_lagged<-get_lag(x = mixturedata$Volume_Species1, time = mixturedata$Day)
Pa_lagged<-get_lag(x = mixturedata$Volume_Species2, time = mixturedata$Day)

# calculate per-capita growth rates
Pc_dNNdt<-percap_growth(x = Pc_lagged$x, laggedx = Pc_lagged$laggedx, dt = Pc_lagged$dt)
Pa_dNNdt<-percap_growth(x = Pa_lagged$x, laggedx = Pa_lagged$laggedx, dt = Pa_lagged$dt)

# fit linear models to dNNdt, based on average
# abundances between current and lagged time steps
Pc_mod_dat<-data.frame(Pc_dNNdt=Pc_dNNdt, Pc=Pc_lagged$laggedx, Pa=Pa_lagged$laggedx)
mod_comp_Pc<-lm(Pc_dNNdt~Pc+Pa, data=Pc_mod_dat)

Pa_mod_dat<-data.frame(Pa_dNNdt=Pa_dNNdt, Pa=Pa_lagged$laggedx, Pc=Pc_lagged$laggedx)
mod_comp_Pa<-lm(Pa_dNNdt~Pa+Pc, data=Pa_mod_dat)

# model summaries
summary(mod_comp_Pc)
summary(mod_comp_Pa)

# extract parameters
# note - linear regressions give us dynamics in the form:
# dni/nidt ~ (Intercept) + (n1_slope) * n1 + (n2_slope) n2
# and thus:
# dni/dt = n1*((Intercept) + (n1_slope) * n1 + (n2_slope) n2)

# growth rates
r1 <- unname(coef(mod_comp_Pc)["(Intercept)"])
r2 <- unname(coef(mod_comp_Pa)["(Intercept)"])

# self-limitation
a11 <- unname(coef(mod_comp_Pc)["Pc"])
a22 <- unname(coef(mod_comp_Pa)["Pa"])

# effect of Pa on Pc
a12 <- unname(coef(mod_comp_Pc)["Pa"])
# effect of Pc on Pa
a21 <- unname(coef(mod_comp_Pa)["Pc"])

# run ODE:
# make paramter vector:
parms <- c(r1, r2, a11, a12, a21, a22)
initialN <- c(1, 1)
out <- deSolve::ode(y=initialN, times=1:25, func=lv_interaction, parms=parms)
matplot(out[,1], out[,-1], type="l",
   xlab="time", ylab="N", col=c("black","red"), lty=c(1,3), lwd=2, ylim=c(0, 150))
legend("topleft", c("Pc", "Pa"), col=c(1,2), lwd=2, lty=c(1,3))

# now, plot in points from data
points(mixturedata$Day, mixturedata$Volume_Species1, col=1)
points(mixturedata$Day, mixturedata$Volume_Species2, col=2)