--- layout: page title: Competitive Lotka–Volterra Model ---

Lotka–Volterra Equations

\[ \frac{dx_1}{dt}=r_1x_1(1-\frac{x_1 + \alpha_{12}x_2}{K_1})\] \[ \frac{dx_2}{dt}=r_2x_2(1-\frac{x_2 + \alpha_{21}x_1}{K_2})\]

Case 1: Stable Equilibrium (with random noises)

\[K_1 < \frac{K_2}{\alpha_{21}}\] \[K_2 < \frac{K_1}{\alpha_{12}}\]

Set \(K_1=100\), \(K_2=100\), \(\alpha_{12}=0.6\), \(\alpha_{21}=0.6\),

\(r_1 = 0.1\), \(r_2 = 0.1\),

and initial population size: \(x_1=5\), \(x_2=10\).

Time series

plot_DE(cp_lotka(r1=0.1,r2=0.1, a12=0.6,a21=0.6, K1=100,K2=100, x1_int=5, x2_int=10))
## [1] 62.5 62.5

Phase portrait

p <- c(r1=0.1,r2=0.1, a12=0.6,a21=0.6, K1=100,K2=100) # p is a named vector of parameters
s <- c(x1=5,x2=10) # s is the state
plane(tstep=0.5,portrait=T, xmax=150,ymax=150)

Case 2: Unstable Equilibrium (with random noises)

\[K_1 > \frac{K_2}{\alpha_{21}}\] \[K_2 > \frac{K_1}{\alpha_{12}}\] Set \(K_1=100\), \(K_2=100\), \(\alpha_{12}=1.2\), \(\alpha_{21}=1.2\),

and same initial population size for the two species: \(x_1=10\), \(x_2=10\).

Time series: species coexist

Set \(r_1 = 0.1\), \(r_2 = 0.1\)

plot_DE(cp_lotka(r1=0.1,r2=0.1, a12=1.2,a21=1.2, K1=100,K2=100, x1_int=10, x2_int=10))
## [1] 45.45455 45.45455

Time series: species 2 lose

Set \(r_1 = 0.1001\), \(r_2 = 0.1\)

plot_DE(cp_lotka(r1=0.1001,r2=0.1, a12=1.2,a21=1.2, K1=100,K2=100, x1_int=10, x2_int=10))
## [1] 45.45455 45.45455

Phase portrait

p <- c(r1=0.1,r2=0.1, a12=1.2,a21=1.2, K1=100,K2=100) # p is a named vector of parameters
s <- c(x1=10,x2=10) # s is the state
plane(tstep=0.5,portrait=T, xmax=150,ymax=150)

Case 3: Species 2 must lose (with random noises)

\[K_1 > \frac{K_2}{\alpha_{21}}\] \[K_2 < \frac{K_1}{\alpha_{12}}\]

Population size of species 2 has weaker negative impact on the growth rate of species 1 than vice versa.

Set \(K_1=100\), \(K_2=100\), \(\alpha_{12}=0.6\), \(\alpha_{21}=1.2\),

\(r_1 = 0.1\), \(r_2 = 0.15\),

and initial population size: \(x_1=10\), \(x_2=15\).

Time series

plot_DE(cp_lotka(r1=0.1,r2=0.15, a12=0.6,a21=1.2, K1=100,K2=100, x1_int=10, x2_int=15))
## [1] 142.85714 -71.42857

Phase portrait

p <- c(r1=0.1,r2=0.15, a12=0.6,a21=1.2, K1=100,K2=100) # p is a named vector of parameters
s <- c(x1=10,x2=15) # s is the state
plane(tstep=0.5,portrait=T, xmax=150,ymax=150)

Case 4: Species 1 must lose (with random noises)

\[K_1 < \frac{K_2}{\alpha_{21}}\] \[K_2 > \frac{K_1}{\alpha_{12}}\]

Population size of species 1 has weaker negative impact on the growth rate of species 2 than vice versa.

Set \(K_1=100\), \(K_2=100\), \(\alpha_{12}=1.2\), \(\alpha_{21}=0.6\),

\(r_1 = 0.3\), \(r_2 = 0.1\),

and initial population size: \(x_1=30\), \(x_2=10\).

Time series

plot_DE(cp_lotka(r1=0.3,r2=0.1, a12=1.2,a21=0.6, K1=100,K2=100, x1_int=30, x2_int=10))
## [1] -71.42857 142.85714

Phase portrait

p <- c(r1=0.3,r2=0.1, a12=1.2,a21=0.6, K1=100,K2=100) # p is a named vector of parameters
s <- c(x1=30,x2=10) # s is the state
plane(tstep=0.5,portrait=T, xmax=150,ymax=150)

Code

Function to generate data

cp_lotka <- function(r1,r2,a12,a21,K1,K2,x1_int,x2_int){
    library(deSolve)
    
    Pars <- c(r1,r2,a12,a21,K1,K2)
    State <- c(x1=x1_int, x2=x2_int)
    
    
    LotVmod <- function (Time, State, Pars) {
        with(as.list(c(State, Pars)), {
            dx1 = (1+rnorm(1,0,0.01))*r1*x1*(1-(x1+a12*x2)/(K1+rnorm(1,0,0.05)))   # add random noise: rnorm(1,0,0.01) to replace 0
            dx2 = (1+rnorm(1,0,0.01))*r2*x2*(1-(x2+a21*x1)/(K2+rnorm(1,0,0.05)))
            return(list(c(dx1, dx2)))
        })
    }
    ## equilibrium point
    x1 <- (K1-a12*K2)/(1-a12*a21)
    x2 <- (K2-a21*K1)/(1-a12*a21)
    print(c(x1,x2))

    Time <- seq(0, 1000, by = 1)
    out <- as.data.frame(ode(func = LotVmod, y = State, parms = Pars, times = Time))
}

Function to plot 2 species time series

plot_DE <- function(out){
    p <- ggplot(out[,-1], aes(x=1:nrow(out)))+
            geom_line(aes(y=x1,color="Species 1"))+
            geom_line(aes(y=x2,color="Species 2"))
    p + labs(x = "Generation",y="N",colour = "Species")
}

Function for phase portrait

source('PhasePortraitFc.R')

model <- function (Time, State, Pars) {
    with(as.list(c(State, Pars)), {
        dx1 = (1+rnorm(1,0,0.01))*r1*x1*(1-(x1+a12*x2)/(K1+rnorm(1,0,0.05)))   # add random noise: rnorm(1,0,0.01) to replace 0
        dx2 = (1+rnorm(1,0,0.01))*r2*x2*(1-(x2+a21*x1)/(K2+rnorm(1,0,0.05)))
        return(list(c(dx1, dx2)))
    })
}

PhasePortraitFc.R: credit to Rob de Boer