# 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. ```r 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. ## Returns vector of growth rates for each species ## Examples ```r # 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