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')
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]:
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)
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.
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]:
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)
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)
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]:
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()
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}$")
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")
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")