Modele Prada - Pradator

Modelele de tip prada-pradator modeleaza interactiunea dintre doua specii, una din ele reprezentand pradatorul, cel care beneficiaza de pe urma interactiunii, iar cealalta specie fiind prada, cea care pierde de pe urma interactiunii. Acest tip de modele au un rol important in ecologie si prezinta dinamica populatiilor din ecosisteme.

Modelul clasic Lotka-Volterra

Cel mai simplu model pentru interactiunea prada-pradator este cel propus de Lotka si Volterra si este definiti de 2 ecuatii diferentiale ordinare:

unde:

• x reprezinta numarul de indivizi ce reprezinta prada
• y reprezinta numarul de indivizi ce reprezinta pradatorii
• t reprezinta timpul
• α, β, δ, γ reprezinta parametrii sistemului ce descriu interactiunea dintre cele 2 specii (numere reale pozitive)
• si reprezinta ratele de modificare (crestere sau scadere) instatanee ale celor 2 populatii

Pentru a intelege ce semnifica cele 2 ecuatii, le vom lua pe fiecare pe rand pentru a le analiza.

Ecuatia evolutiei populatiei prada

• αx reprezinta functia de crestere a populatiei de prada in absenta pradatorilor. Astfel Δx = αx defineste o functie exponentiala de crestere, unde rata de crestere a populatiei este direct proportionala cu numarul de indivizi care exista in momentul respectiv.

  Astfel populatia de prada se presupune ca va creste nelimitat, intr-o maniera exponentiala in lipsa pradatorilor. De altfel, aceasta presupunere este naiva, iar vom incerca imbunatatirea modelului mai tarziu, folosind o functie diferita de crestere.

• -βxy reprezinta functia raspuns la activitatea de vanatoare a pradatorilor asupra populatiei de prada. Este o functie ce depinde atat de numarul de pradatori, cat si de numarul de prazi care exista in momentul respectiv. Are semn negativ deoarece populatia de prada are de pierdut in urma interactiunii cu pradatorii.

  Astfel, activitatea de vanatoare a pradatorilor limiteaza si induce un trend descrescator in populatia prazii.

Ecuatia evolutiei populatiei pradator

• -γy reprezinta functia de scadere a populatiei pradatorilor in absenta prazii. Astfel Δy = -γy defineste o functie exponentiala de scadere, unde rata de scadere a populatiei este direct proportionala cu numarul de pradatori existenti in acel moment.

  Aceasta scadere in randul pradatorilor simbolizeaza faptul ca, in lipsa unei prazi pe care o pot vana, pradatorii vor muri intr-un ritm exponential deoarece nu vor avea ce manca.

• δxy reprezinta functia raspuns la activitatea de vanatoare a pradatorilor asupra populatiei lor. Functia depinde atat de numarul de pradatori, cat si de numarul de prazi care exista in momentul respectiv. Semnul pozitiv al functiei simbolizeaza o crestere a numarului de pradatorii in urma interactiunii cu prada.

  Astfel, doar in urma actiunii de vanatoare, populatia pradatorilor poate sa creasca

Implementarea modelului

Definim in R modelul Lotka-Volterra

lv_model <- function(t, state, params) {
  with(as.list(c(state, params)), {
    d_prey <- alpha * prey - beta * prey * predator
    d_predator <- delta * prey * predator - gamma * predator
    return(list(c(prey = d_prey, predator = d_predator)))
  })
}

Dam valori parametrilor ce caracterizeaza modelul si un numar initial al populatiei prada, respectiv pradator si putem rezolva pe un interval de timp sistemul ordinal de ecuatii diferentiale utilizand functia ode din bilbioteca deSolve.

library(deSolve)

parameters <- c(alpha = 0.2, beta = 0.01, delta = 0.005, gamma = 0.1)
init <- c(prey = 10, predator = 5)
times <- seq(0, 150, by=1)

lv_results <- ode(init, times, lv_model, parameters)

Dinamica populatiilor in timp se poate observa in graficul rezultat.

param_names = paste(names(parameters), parameters, sep = " = ", collapse = "; ")
title = "Model Lotka-Volterra, dinamica in timp a populatiilor"
plot(lv_results[,1], lv_results[,2], main=title, sub=param_names,col="red",type="l", lwd=2, xlab="Timp",ylab="Populatie")
points(lv_results[,1], lv_results[,3], col="blue",type="l",lwd=2, lty=2)

legend("topright",lwd=c(2,2),lty=c(1,2),col=c("red","blue"),legend=c("prada","pradator"),bty="n")

param_names = paste(names(parameters), parameters, sep = " = ", collapse = "; ")
title = "Model Lotka-Volterra, diagrama de faza a populatiilor"
plot(lv_results[,2], lv_results[,3], main=title, sub=param_names, col="#f1a340", type="l", lwd=2, xlab="Prada", ylab="Pradator")
points(lv_results[1,2], lv_results[1,3], col="#f1a340", pch=20, cex=3)

Diagrama de faza a populatiilor ne arata ca exista o relatie periodica intre populatiile de prada si pradator. De altfel, daca se iau Δx si Δy ca 0 in cele 2 ecuatii diferentiale, se pot calcula starile de echilibru ale sistemului. Astfel exista 2 stari de echilibru:

• (0, 0), ceea ce reprezinta starea in care nu exista nici un individ, fie prada sau pradator, si s-ar intelege de ce e o stare de echilibru
• (α/β, γ/δ), este punctul in care dezavantajul interactiunii cu pradatorul compenseaza exact cresterea naturala a populatiei prazii, respectiv avantajul obtinut din interactiunea cu prada compenseaza exact scaderea naturala a populatiei pradatorilor. Pentru modelul nostru parametrizat cu α=0.2, β=0.01, γ=0.1, δ=0.005 acest punct este (20, 20)

Punctele de echilibru se afla la intersectia celor 2 linii unde ecuatiile diferentiale din sistem se anuleaza (nullcline). Putem observa ca pentru acelasi model, indiferent de populatiile initiale, cele 2 populatii vor oscila periodic in jurul punctului de echilibru (α/β, γ/δ).

library(phaseR)

param_names = paste(names(parameters), parameters, sep = " = ", collapse = "; ")
title = "Model Lotka-Volterra, diagrama de faza a populatiilor"
xlim = c(0, 120)
ylim = c(0, 80)

plot(20, 20, main=title, sub=param_names, xlim=xlim, ylim=ylim, pch=20, xlab="Prada", ylab="Pradator")

invisible(nullclines(lv_model, xlim=xlim, ylim=ylim, parameters = parameters, col=c("red", "blue"), lty=5, system="two.dim", state.names=c("prey", "predator")))

invisible(flowField(lv_model, xlim=xlim, ylim=ylim, parameters = parameters, system="two.dim", state.names=c("prey", "predator")))


inits <- rbind(c(3,3), c(15, 18), c(30, 40), c(55, 45))
for(i in 1:nrow(inits)) {
  init <- c(prey = inits[i, 1], predator = inits[i, 2])
  lvr <- ode(init, times, lv_model, parameters)
  
  lines(lvr[,2], lvr[,3], main=title, sub=param_names, col="#f1a340", type="l", lwd=2, xlab="Prada", ylab="Pradator")
  points(lvr[1,2], lvr[1,3], col="#f1a340", pch=20, cex=3)
}

Imbunatatiri ale modelului clasic

Din multe privinte modelul clasic prada-pradator Lotvka-Volterra este naiv si prezinta lucrurile intr-o maniera mult simplificata. De altfel, el functioneaza pe baza urmatoarelor presupuneri:

• populatia de prada se inmulteste la nesfarsit si este limitata doar de interactiunea cu pradatorii
• pradatorii pot manca la infinit din populatia prada, iar infometarea lor nu e satisfacuta niciodata
• populatia de pradatori supravietuiesc hranindu-se doar cu o specie de prada
• nu exista competitie intre pradatori pentru prada
• prada nu se poate ascunde din calea pradatorilor

1. Introducem capacitate maxima pentru populatia de prada

Folosind o functie de crestere logistica, in loc de cea exponentiala, putem limita populatia de prada la o capacitate maxima. Acesta este un scenariu mult mai realist, pentru ca in lumea reala, chiar si in lipsa pradatorilor, exista constrangeri legate de resursele existente.

Modelul logistic de crestere unde K reprezinta capacitatea maxima.

Ecuatiile modelului vor arata astfel acum:

Definim modelul in R:

lv_improved_1 <- function(t, state, params) {
  with(as.list(c(state, params)), {
    d_prey <- alpha * prey * (1 - prey / kappa) - beta * prey * predator
    d_predator <- delta * prey * predator - gamma * predator
    return(list(c(prey = d_prey, predator = d_predator)))
  })
}

Simularea modelului:

parameters_1 <- c(alpha = 0.6, beta = 0.01, delta = 0.005, gamma = 0.1, kappa = 250)
init_1 <- c(prey = 10, predator = 5)
times_1 <- seq(0, 150, by=1)

lv_improved_1_results <- ode(init_1, times_1, lv_improved_1, parameters_1)

param_names = paste(names(parameters_1), parameters_1, sep = " = ", collapse = "; ")
title = "Model prada-pradator imbunatatit 1, dinamica in timp a populatiilor"
plot(lv_improved_1_results[,1], lv_improved_1_results[,2], main=title, sub=param_names,col="red",type="l", lwd=2, xlab="Timp",ylab="Populatie")
points(lv_improved_1_results[,1], lv_improved_1_results[,3], col="blue",type="l",lwd=2, lty=2)

legend("topright",lwd=c(2,2),lty=c(1,2),col=c("red","blue"),legend=c("prada","pradator"),bty="n")

In graficul ce reprezinta dinamica in timp a populatiilor se observa cum oscilatiile se reduc si tind spre o stare de echilibru.

param_names = paste(names(parameters_1), parameters_1, sep = " = ", collapse = "; ")
title = "Model prada-pradator imbunatatit 1, diagrama de faza a populatiilor"
xlim = c(0, 160)
ylim = c(0, 130)

plot(20, 55.2, main=title, sub=param_names, xlim=xlim, ylim=ylim, pch=20, xlab="Prada", ylab="Pradator")

invisible(nullclines(lv_improved_1, xlim=xlim, ylim=ylim, parameters = parameters_1, col=c("red", "blue"), lty=5, system="two.dim", state.names=c("prey", "predator")))

invisible(flowField(lv_improved_1, xlim=xlim, ylim=ylim, parameters = parameters_1, system="two.dim", state.names=c("prey", "predator")))

colors <- c("#f1a340",  "#998ec3")
inits_1 <- rbind(c(3,3), c(95, 18))
for(i in 1:nrow(inits_1)) {
  init <- c(prey = inits_1[i, 1], predator = inits_1[i, 2])
  lvr <- ode(init, times_1, lv_improved_1, parameters_1)
  
  lines(lvr[,2], lvr[,3], main=title, sub=param_names, col=colors[i], type="l", lwd=2, xlab="Prada", ylab="Pradator")
  points(lvr[1,2], lvr[1,3], col=colors[i], pch=20, cex=3)
}

Diagrama de faza a populatiilor reflecta faptul ca oscilatiile populatiilor nu sunt periodice, iar acestea tind spre unul din punctele de echilibru unde atat pradatorii, cat si prazile se afla in numar egal.

Pentru modelul acesta se pot calcula 3 stari de echilibru: - (0,0) - (γ/δ, α(δK-γ)/(βδK)) care la nivelul diagramei de mai sus reprezinta punctul (20, 55.2) - (K, 0) poate cel mai interesat care da de inteles ca exista o stare in care pradatorii au disparut si mai exista doar populatie prada

2. Schimbam functia de raspuns la interactiunea prada-pradator

In modelul actual se presupune ca numarul de prazi consumate de catre pradatori creste intr-o maniera lineara pe masura ce populatia de prazi creste. Presupunerea acesata este naiva pentru ca inseamna ca pradatorii nu se satura niciodata. Acest tip de functie de raspuns are tipul I.Pe langa aceasta, voi prezenta si functia de raspuns de tipul II si III.

type1F <- function(x, beta) {
  x*beta
}

type2F <- function(x, beta, h) {
  (beta*x)/(1+beta*h*x)
}

type3F <- function(x, beta, h) {
  (beta*x*x)/(1+beta*h*x*x)
}

Tipul I: unde:

• β caracterizeaza rata de atac a pradatorului

curve(type1F(x, 0.1), 0, 1000, main="Functie raspuns la interactiunea prada-pradator de tip I", col="red", type="l", lwd=2, xlab="Populatie prada",ylab="Indivizi prada mancati / pradator")

Tipul II: unde:

• β caracterizeaza rata de atac a pradatorului
• h reprezinta timpul necesar ca pradatorul sa isi atace victima

Acest raspuns tine cont de satietatea pradatorului. Numarul de prazi consumate de pradator creste proportional cu populatia prazilor, dar pana la un punct.

curve(type2F(x, 0.1, 0.1), 0, 1000, main="Functie raspuns la interactiunea prada-pradator de tip II", col="red", type="l", lwd=2, xlab="Populatie prada",ylab="Indivizi prada mancati / pradator")

Tipul III: unde:

• β caracterizeaza rata de atac a pradatorului
• h reprezinta timpul necesar ca pradatorul sa atace o prada

In plus fata de raspunsul de tip II, pradatorul isi formeaza o imagine a prazii pe masura prezentei acesteia. Astfel, daca populatia prada este redusa, pradatorul nu o va cunoaste si nu o va cauta, dar pe masura ce apar tot mai multe prazii, se va familiariza cu aceasta.

curve(type3F(x, 0.0001, 0.1), 0, 1000, main="Functie raspuns la interactiunea prada-pradator de tip III", col="red", type="l", lwd=2, xlab="Populatie prada",ylab="Indivizi prada mancati / pradator")

Folosind functia raspuns de tip III, modelul imbunatatit va fi:

lv_improved_2 <- function(t, state, params) {
  with(as.list(c(state, params)), {
    d_prey <- alpha * prey * (1 - prey / kappa) - type3F(prey, beta, h) * predator
    d_predator <- type3F(prey, beta, h) * predator - gamma * predator
    return(list(c(prey = d_prey, predator = d_predator)))
  })
}

parameters_2 <- c(alpha = 0.1, beta = 0.0012, delta = 0.05, gamma = 0.9, kappa = 250, h = 0.1)
init_2 <- c(prey = 10, predator = 5)
times_2 <- seq(0, 150, by=1)

lv_improved_2_results <- ode(init_2, times_2, lv_improved_2, parameters_2)

param_names = paste(names(parameters_2), parameters_2, sep = " = ", collapse = "; ")
title = "Model prada-pradator imbunatatit 2, dinamica in timp a populatiilor"
plot(lv_improved_2_results[,1], lv_improved_2_results[,2], main=title, sub=param_names,col="red",type="l", lwd=2, xlab="Timp",ylab="Populatie", ylim=c(0, 50))
points(lv_improved_2_results[,1], lv_improved_2_results[,3], col="blue",type="l",lwd=2, lty=2)

legend("topright",lwd=c(2,2),lty=c(1,2),col=c("red","blue"),legend=c("prada","pradator"),bty="n")

param_names = paste(names(parameters_2), parameters_2, sep = " = ", collapse = "; ")
title = "Model prada-pradator imbunatatit 2, diagrama de faza a populatiilor"
plot(lv_improved_2_results[,2], lv_improved_2_results[,3], main=title, sub=param_names, col="#f1a340", type="l", lwd=2, xlab="Prada", ylab="Pradator")
points(lv_improved_2_results[1,2], lv_improved_2_results[1,3], col="#f1a340", pch=20, cex=3)

Bibliografie

• The prey-predator model: https://rpubs.com/Jeet1994/Prey-predator-model
• Mathematical biology, Lecture Notes for MATH 4333, Jeffrey, R. Chasnov, Chapter 1.4. “The Lotka-Volterra predator-prey model”: https://www.math.ust.hk/~machas/mathematical-biology.pdf
• Applied Population Ecology, Lecture Notes, Kevin Shoemaker, Species Interactions: prey-predator: https://kevintshoemaker.github.io/NRES-470/LECTURE17.html
• University of California Irvine OpenCourseWare, Introduction to Mathematical Modeling in Biology: Predator Prey Model: http://ocw.uci.edu/lectures/math_113b_lec_14_introduction_to_mathematical_modeling_in_biology_predator_prey_model.html