Modelo Rosenzweig-MacArthur

Un modelo razonable para las poblaciones de presas y depredadores a través del tiempo, considerando su interacción, podría tener esta forma:

$\dot{x}=$ tasa de nacimiento - tasa de muertes no debidas a $ y $ - tasa de muertes por depredación

$\dot{y}=$ tasa de reproducción - tasa de mortalidad

Donde $ x $ representa la densidad de presas (número de presas por unidad de área), y $ y $ la densidad de depredadores (número de depredadores por unidad de área).

Un modelo cualquiera es una elección de cómo computar estas tasas, así lo abstrae el modelo de Rosenzweig-McArthur:

$$\begin{cases}\dot{x}=f(x)-\phi(x,y)\\\dot{y}=-ey+k\phi(x,y)\end{cases}$$

Donde $f(x)$ es la tasa de cambio de $x$ cuando no hay depredación, $\phi(x, y)$ es la tasa de depredación, $ k $ la eficiencia con la que los depredadores asimilan a las presas y $ e $ la tasa de muertes del depredador ($ e, k>0 $).


In [1]:
# Para hacer experimentos numéricos importamos numpy
import numpy as np

# y biblioteca para plotear
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline

# cómputo simbólico con sympy
from sympy import *
init_printing()

# definimos qué variables usar como símbolos
x, y, e, phi, f, k, r, K = symbols('x y e phi f k r K')

Tasa de cambio de la población de presas

La tasa de cambio de $ x $ (nacimientos y muertes de presas) podría modelarse como crecimiento Maltusiano $f(x)=\frac{dx}{dt}=rx\left(1-\frac{x}{K}\right)$, donde la tasa de reproducción es proporcional tanto al tamaño de la población actual como a la cantidad de recursos disponible ($ r,K>0 $).


In [2]:
malthusian = r*x*(1-x/K)
malthusian


Out[2]:
$$r x \left(1 - \frac{x}{K}\right)$$

In [3]:
F = lambdify((x, r, K), malthusian)
pop = np.zeros(100)

# población inicial
pop[0] = 1
reproduccion = 0.1
capacidad = 200

# sumamos población en cada t
for t in range(99):
    pop[t+1] = pop[t] + F(pop[t], 
                          reproduccion,
                          capacidad)
    
plt.xlabel("tiempo")
plt.ylabel("población de presas sin depredación")
fig = plt.plot(pop)


Tasa de cambio de la población de depredadores

Si se asume que la tasa de muertes del depredador es lineal y su densidad no es tal que se interfieran entre sí para cazar, entonces sólo la cantidad de presas limita la población de depredadores.

$$\begin{cases}\dot{x}=f(x)-y\phi(x)\\\dot{y}=-ey+ky\phi(x)\end{cases}$$

$\phi(x)$ es la cantidad de presas atrapadas por depredador por unidad de tiempo, en función de la densidad de presas.

Tasa de depredación según Holling

Lo esencial es que la tasa de depredación esté en función sólo de la densidad de presas. Este artículo discute la creación de tal función considerando el promedio de presas atrapadas a lo largo de un lapso extendido. El argumento $a$ modela el tiempo que le toma al depredador encontrar una presa, mientras que $b$ modela cuánto tiempo le toma atraparla, comérsela, digerirla, dormir la siesta, y empezar a acechar nuevamente.

$$\phi(x)=\frac{ax}{1+bx}$$

In [4]:
a, b = symbols('a b')
phi_holling = a*x / (1 + b*x)
phi_holling


Out[4]:
$$\frac{a x}{b x + 1}$$

In [5]:
# convertimos los símbolos de phi_holling en función ejecutable
phi_holling = lambdify((x, a, b), phi_holling)

Con un tiempo para encontrar presa de 20 unidades y un tiempo para consumirlas de 9 encontramos una curva parecida a la del artículo de Hal L Smith.


In [6]:
t_encontrar = 20
t_consumir = 9

presas = np.linspace(0, 3.5, 100)

# computamos depredaciones para un rango de densidades de presas
depredaciones = [phi_holling(n, 
                             t_encontrar,
                             t_consumir) 
                 for n in presas]

plt.xlabel("densidad de presas")
plt.ylabel("tasa de depredación")
fig = plt.plot(presas, depredaciones)


Interacción de ambas poblaciones

La sección 2 del mismo artículo, Hal L. Smith discute una forma más conveniente del mismo modelo en la que los parámetros de tiempo $a$ y $b$ se reemplazan por el argumento adimensional $m$.

$$\begin{cases}\dot{x}=x\left(1-\frac{x}{k}\right)-\frac{mxy}{1+x}\\\dot{y}=-ey+\frac{mxy}{1+x}\end{cases}$$

In [7]:
x, y, k, m, e = symbols('x y k m e')

dx = x*(1-x/k) - (m*x*y)/(1+x)
dy = (-e*y)+(m*x*y)/(1+x)

Equilibrios


In [8]:
# ¿en qué lugares las tasas de cambio son iguales a cero?
f0 = Eq(dx, 0)
g0 = Eq(dy, 0)
equilibrios = solve((f0, g0), x, y)

# dos equilibrios dependen de los valores de e, m y k
eq3x = lambdify((e,m), equilibrios[2][0])
eq3y = lambdify((e,m,k), equilibrios[2][1])

equilibrios


Out[8]:
$$\left [ \left ( 0, \quad 0\right ), \quad \left ( k, \quad 0\right ), \quad \left ( - \frac{e}{e - m}, \quad \frac{- e k - e + k m}{k \left(e - m\right)^{2}}\right )\right ]$$

In [9]:
def presas(m, k):
    """
    devuelve una función de (x, y)
    con los argumentos (m, k) interpolados
    """
    def f(x, y):
        return x*(1-x/k)-((m*x*y)/(1+x))
    return f

def depredadores(e, m):
    """
    devuelve una función de (x, y)
    con los argumentos (e, m) interpolados
    """
    def g(x, y):
        return (-e*y)+((m*x*y)/(1+x))
    return g

In [10]:
m = 3
k = 3
e = 1
f = presas(m, k)
g = depredadores(e, m)

In [11]:
x, y = np.meshgrid(np.linspace(0, 3.5, 22),
                   np.linspace(-0.5, 1, 22))

u = f(x, y)
v = g(x, y)

# campo de vectores
plt.quiver(x, y, u, v)

# equilibrios
plt.plot(0, 0, 'red', marker='o', markersize=10)
plt.plot(k, 0, 'red', marker='o', markersize=10)
plt.plot(eq3x(e,m), eq3y(e,m,k), 'red', marker='o', markersize=10)

plt.xlabel("$\dot{x}$")
plt.ylabel("$\dot{y}$")

plt.show()


Trayectorias


In [12]:
# funciones para graficar las trayectorias
def paso(x, y, dt, f, g):
    """
    Devuelve una nueva coordenada a partir de:
     - la actual (x, y)
     - delta t
     - funciones de tasa de cambio
    """
    return (x + dt * f(x, y),
            y + dt * g(x, y))

def trayectoria(x0, y0, f, g, dt=0.01, pasos=100):
    """
    Calcula una trayectoria dado un punto inicial, las funciones
    de tasa de cambio, delta t y un número de iteraciones.
    """
    x = x0
    y = y0
    t = list()
    for n in range(pasos):
        t.append((x, y))
        x, y = paso(x, y, dt, f, g)
    return t

In [13]:
m = 3
k = 3
e = 1.5
f = presas(m, k)
g = depredadores(e, m)

# Calcular 800 pasos tamaño 0.1 de las trayectorias
dt = 0.1
pasos = 800

# Condiciones iniciales dentro del ciclo límite, en color índigo
x0 = 1
y0 = 0.4
x_dentro, y_dentro = zip(*[coord for coord in
                          trayectoria(x0, y0,
                                      f, g,
                                      dt, pasos)])

# Condiciones iniciales fuera del ciclo límite, en verde
x0 = 1
y0 = 1
x_fuera, y_fuera = zip(*[coord for coord in
                         trayectoria(x0, y0,
                                     f, g,
                                     dt, pasos)])

In [14]:
# dibuja campo de vectores
x, y = np.meshgrid(np.linspace(0, 2.2, 22),
                   np.linspace(0, 1, 22))
u = f(x, y)
v = g(x, y)
plt.quiver(x, y, u, v)

# dentro del ciclo límite
plt.plot(x_dentro, y_dentro, color='indigo')
plt.arrow(x_dentro[-2], y_dentro[-2],
          x_dentro[-1] - x_dentro[-2], y_dentro[-1] - y_dentro[-2],
          color="indigo", head_width=0.1, head_length=0.1)

# fuera del ciclo límite
plt.plot(x_fuera, y_fuera, color='green')
plt.arrow(x_fuera[-2], y_fuera[-2],
          x_fuera[-1] - x_fuera[-2], y_fuera[-1] - y_fuera[-2],
          color="green", head_width=0.1, head_length=0.1)

label = plt.xlabel("$\dot{x}$")
label = plt.ylabel("$\dot{y}$")


Condicion inicial dentro del ciclo límite


In [15]:
t = np.arange(pasos)

plt.xlabel("tiempo")
plt.ylabel("$\dot{x},\dot{y}$")
fig = plt.plot(t, x_dentro, color="darkmagenta")
fig = plt.plot(t, y_dentro, color="darkorchid")


Condicion inicial fuera del ciclo límite


In [16]:
t = np.arange(pasos)

plt.xlabel("tiempo")
plt.ylabel("$\dot{x},\dot{y}$")
fig = plt.plot(t, x_fuera, color="lawngreen")
fig = plt.plot(t, y_fuera, color="lime")