R/lv_interaction.R

In gauseR: Lotka-Volterra Models for Gause's 'Struggle for Existence'

#' Lotka-Volterra Interactions
#'
#' 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.
#' @param time The time steps corresponding to each observation - exists to interface with ode function, but should be left blank.
#' @param n A vector of species abundances
#' @param 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.
#' @concept Lokta-Volterra
#' @concept Gause
#' @concept interaction
#' @return vector of growth rates for each species
#' @export
#' @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)

lv_interaction <- function(time, n, parms) {
  nsp<-length(n)
  
  r<-parms[1:nsp]
  A<-matrix(parms[-(1:nsp)], nrow=nsp, ncol=nsp)
  
  dndt <- n*(r+n%*%A)
  
  list(c(dndt))
}