I. INTRODUCCIÓN

El modelo presa-depredador es bien conocido y utilizado para describir comportamientos poblacionales de sistemas muy variados desde virus, bacterias, algas, parásitos, insectos y pequeños y grandes mamíferos, así como sinergias adaptativas e inclusive vínculos con y entre vegetales (^{Beauchamp, 2007}; ^{Blanchard et al., 2017}; ^{Dattner et al., 2017}; ^{Barbosa y Castellanos, 2005}; ^{Johannessen, 2014}). La representación analítica de la relación presa-depredador, desarrollada en paralelo por Alfred Lotka y Vito Volterra da lugar a un sistema de ecuaciones diferenciales acopladas, conocidas como ecuaciones de Lotka-Volterra (LV) (Lotka, 1910), ha sido ampliamente estudiada y aplicada (^{Skvortsov et al., 2018}; ^{Bunin, 2017}; ^{Portalier et al., 2019}). Las ecuaciones de LV describen la evolución temporal de las poblaciones de presas y depredadores mediante parámetros externos que establecen el nacimiento y muerte de ambas especies y funciones que regulan la interacción. En su versión más simple, sea N(t) la población de la especie de presa y P(t) la población de la especie depredadora:

Siendo α, β, γ y δ constantes que regulan la evolución de cada especie y su interacción. En ausencia de depredadores, las presas crecen exponencialmente siendo α su tasa de crecimiento y pudiendo interpretarla como la tasa de reproducción de la población. En ausencia de presas, los depredadores decrecen exponencialmente siendo γ su tasa de decrecimiento y pudiendo interpretarse como la tasa de mortalidad. β y δ pueden interpretarse como la susceptibilidad de las presas a ser atrapadas y la eficacia de los depredadores para atrapar presas, respectivamente tal como se realizó en el artículo original de ^{Lotka-Volterra (Lotka, 1910}). Las poblaciones de presas y depredadores, de acuerdo con este planteo, evolucionan cíclicamente. Las presas, en presencia de pocos depredadores, comienzan a crecer, pero la abundancia de presas produce, con un cierto retraso, el crecimiento de los depredadores (Lotka, 1910). Cuando los depredadores son muchos la población de presas comienza a disminuir y esto provoca la disminución de éstos. La frecuencia de las oscilaciones poblacionales y su amplitud dependen de las constantes y de las condiciones iniciales (^{Jiang et al., 2017}).

La posibilidad de describir el comportamiento de cada individuo mediante un sistema multiagente (SMA) avanza sobre el problema presa-depredador desde un enfoque que permite, además del seguimiento temporal de las poblaciones, obtener otros atributos poblacionales de interés, tales como la distribución espacial, la estructura poblacional o la distribución de biomasa y su densidad temporal o espacial (^{Seitbekova y Bakibayev, 2018}). Podemos notar, además, que una diferencia fundamental del modelado por SMA es que éste no es determinista, es decir que, dadas las mismas condiciones, el sistema puede evolucionar de manera distinta, permitiendo que las trayectorias en el espacio de las fases (N, P) se intersequen, dando lugar a ramificaciones (^{Strogatz, 2018}). De esta manera no se presentan necesariamente las órbitas cerradas, típicas de estos sistemas, durante las oscilaciones poblacionales (LV, ^{Lotka, 1910}). El modelo que se presenta permitiría estudiar las condiciones bajo las cuales se producen estas ramificaciones.

El crecimiento y desarrollo de un ecosistema también es un tema de interés actual, entendiendo al crecimiento como el aumento de la biomasa total mientras que la maduración y desarrollo se vinculan a las interacciones que dan lugar a los flujos de biomasa en la estructura trófica (^{Hatton et al., 2015}; ^{Ulanowicz, 1986}). Estos aspectos, desde la descripción de SMA y sus definiciones individuales, podrían describirse y explicarse como comportamientos emergentes más allá de la descripción puramente fenomenológica (^{Gilbert y Trena, 2000}). En este sentido la información y datos empíricos son escasos y con incertidumbres importantes. A modo de descripción de la fenomenología se muestran en la Figura 1, datos empíricos de la abundancia poblacional de liebres y linces durante un período de 90 años (1845 - 1935) en la tundra de Alaska con un período de aproximadamente 9 años. Esta información si bien ha sido cuestionada por surgir del número de cueros de ambas especies, sigue siendo, actualmente, muy utilizada como una serie temporal típica a los fines ilustrativos (^{Odum y Barrett, 1953}).

Figura 1: Abundancia poblacional (cantidad de individuos de cada especie) de liebres (en rojo) y de linces (en azul) desde 1845 a 1935 en la tundra de Alaska. 

II. METODOLOGÍA

El modelo se implementó en la plataforma Matlab. Este modelo es genérico, sin especificar un vínculo presa-depredador particular entre especies. Por lo tanto, el comportamiento que se obtiene es característico de este tipo de relación. El modelo se desarrolla espacialmente en una región plana denominada hábitat caracterizada por un radio R_h. En la simulación los agentes no pueden sobrepasar los límites del hábitat. Inicialmente se encuentran en el hábitat, una cierta cantidad de agentes de dos especies: presas (N ₀) y depredadores (P₀), cada una con características que definirán sus interacciones, tamaño, velocidad y dirección de movimiento y masa. Los agentes se modelan como círculos de radios fijos que determina el espacio que ese agente puede sensar a su paso y se mueven en línea recta en direcciones sorteadas al azar. En su camino se cuentan las interacciones con otros agentes. El modelo, para contar las interacciones entre agentes, calcula la superficie de superposición entre las áreas de los mismos en los encuentros ocurridos durante el viaje o recorrido realizado. Las distintas interacciones (Iij) se acumulan de la siguiente manera: superficie total de interacciones con otros agentes de la misma especie (I_1,1: presa-presa. I_2,2: depredador-depredador) y superficie total de interacciones con agentes de especie diferente (I_{1,2 =} I_2,1: presa-depredador).

El movimiento de los agentes se divide en intervalos temporales de duración igual a ∆T al cabo del cual se evalúa:

las interacciones que definen si hay generación de nuevos agentes y/o desaparición de algunos de los existentes

el conteo de los agentes de distintas especies

la posición de los agentes. Si han sobrepasado los límites del hábitat, se sortea una nueva dirección de viaje hacia el interior del hábitat

la estadística de la distribución espacial (^{Papoulis, 1991}; ^{Cáceres, 2002}; ^{Spiegel y Stephens, 2009}) de los agentes: momentos de primer, segundo, tercer y cuarto orden, los cuales dan información sobre la posición media, dispersión, asimetría alrededor de la media y curtosis, respectivamente.

Al cabo de cada ∆T se evalúa la cantidad de agentes nuevos de cada especie y la cantidad de agentes que desaparecen. El modelo no contempla procesos tales como inmigración o emigración. Los agentes que representan a las presas tienen un crecimiento determinado por el número de presas existentes y un parámetro r₁, equivalente al α del modelo de LV (Ec. 1): ∆N(T) = r₁ N(T) ∆T. Para determinar la posición inicial donde aparecen esos nuevos agentes presas, se evalúa la trayectoria de su generador. Los agentes generadores serán aquellos que durante la última iteración tuvieron mayor interacción con otros agentes de la misma especie. El punto medio del desplazamiento durante el último ∆T del generador, determina la posición inicial de su progenie. Este tipo de cómputo, que tiene implícito un promedio en las especies y sus interacciones, implica que ∆T represente un tiempo comparable al de los ciclos de reproducción.

Para el caso de los agentes depredadores la situación es diferente al planteo de LV. En el presente modelo, para calcular el crecimiento del número de agentes depredador, utilizamos un parámetro a₂ multiplicado por la superficie total de todas las interacciones entre presas y depredadores que ocurrieron durante ese intervalo (I_1,2 ). Este resultado se ve afectado no sólo por la cantidad de agentes sino también por sus distribuciones espaciales. Los progenitores de los nuevos agentes depredador serán los que hayan acumulado mayor interacción con agentes presa (contribución de cada agente depredador a I_1,2). La posición se define de manera análoga al caso de los nuevos agentes presa ya descripto.

De manera similar, para el cálculo del número de agentes presa que desaparecen al cabo de un intervalo ∆T se utiliza la acumulación de interacciones de los agentes presa con depredadores (I_1,2) y se multiplica por un parámetro a₁. Los agentes presa que desaparecen serán los que hayan acumulado mayor interacción con depredadores (contribución de cada agente presa a I_1,2).

Los diferentes agentes también tienen una masa asociada. En el caso de las presas esa masa es la misma para todos los agentes de la especie (m₀). Para mantenerse con vida, los agentes presa tienen una fuente inagotable de alimento disponible y consumen una masa de alimento proporcional a su masa individual, con factor f_v. Ese flujo de masa entre el reservorio y la especie presa también considera la masa necesaria para la generación de los nuevos agentes.

La interacción entre presas y depredadores genera un flujo de masa. El flujo de masa global de las presas hacia los depredadores es el correspondiente a la masa de las presas muertas. Este flujo se distribuye entre lo necesario para la aparición de los nuevos agentes depredador y para mantener con vida a los ya existentes. Los depredadores requieren, para mantenerse, un flujo de masa mínimo proporcional a su propia masa a través de un factor f_v, mientras que los agentes nuevos nacen cada uno con una masa predeterminada (M₀).

Una vez calculado el flujo de masa de presas hacia depredadores, el cual está determinado por el parámetro a₁ y las interacciones entre los agentes (I_1,2), se calcula la diferencia entre ese flujo y los requerimientos para la generación de los nuevos agentes y el mantenimiento de los existentes. Si la diferencia es positiva, el flujo restante se distribuye por igual entre la totalidad de los agentes aumentando su masa. Si, por el contrario, la diferencia resulta negativa, los agentes disminuyen su masa en lugar de aumentar. Cuando la masa de algún agente depredador cae por debajo de un cierto umbral = μ_cM₀, ese agente desaparece. No hay en el modelo muerte de depredadores con un factor de proporcionalidad al número de depredadores existentes, es decir no existe el término equivalente al −γP de LV (Ec. 2). Los agentes depredadores mueren por no encontrar presas suficientes para su mantenimiento. La masa de los agentes depredadores que mueren constituye un flujo hacia el reservorio ilimitado de donde las presas se alimentan, cerrando el ciclo.

Se puede notar que el modelo propuesto puede adaptarse a distintas situaciones mediante la elección de distintas combinaciones de sus parámetros. La tabla 1 muestra un conjunto de parámetros y datos iniciales con los que será ejemplificado la versatilidad del mismo. 1/3

Tabla 1 Parámetros de la simulación. Unidades arbitrarias. 

Estos parámetros están expresados en unidades arbitrarias, ya que son un ejemplo abstracto, donde la masa de los depredadores es 5 veces la masa de las presas, m₀= 1 mientras M₀= 5, pero sin especificar cuánto valen esas masas en kilogramos. Lo mismo ocurre con los tiempos, que según la especie que en su momento se esté modelando, se asignarán parámetros de crecimiento y decrecimiento adecuados para ser representados en tiempos comparables a los de los ciclos de reproducción.

III. RESULTADOS Y ANÁLISIS

A continuación, se muestran los resultados obtenidos con la simulación en base a un conjunto de parámetros y datos iniciales con los que podemos ejemplificar la información que se obtiene y la versatilidad del modelo. Los mismos se consignan en la Tabla 1. El comportamiento del modelo multiagente se muestra en la figura 2. En 2a) se presenta la evolución temporal del número de agentes presa y depredador. El comportamiento obtenido es semejante a un comportamiento periódico con un período aproximado de 50 ciclos de reproducción y valores máximos de número de agentes presas que se repite en aproximadamente 300 ciclos. Podemos notar también en esta figura, retrasos entre el número de agentes presas y agentes depredador y oscilaciones acopladas. En la figura 2b) se representa la evolución temporal en el espacio de las fases junto a la correspondiente predicha por LV para una situación similar con α=0,2; β=0,27; γ=0,9 y δ=0,05 de la ec (1), en la trayectoria cerrada de color verde. En esta imagen se observa la diferencia entre un comportamiento determinista, dada por la órbita cerrada en verde (^{Strogatz, 2018}), y nuestro modelo estocástico. En la figura 2c) se muestra la evolución de las biomasas de presas y depredadores donde se observa una evolución con quiebres menos pronunciados en la biomasa que en el número de depredadores. En la figura 2d) se representa esta evolución en el espacio de fases para las biomasas totales de presas y depredadores.

Figura 2: a) Evolución temporal del número de agentes presa y depredador b) Evolución del sistema en el espacio de las fases para el número de agentes presa y predador. En verde LV para α=0,2; β=0,27; γ=0,9 y δ=0,05 c) Evolución temporal de la biomasa de los agentes presa y depredador d) Evolución del sistema en el espacio de las fases para la biomasa de agentes presa y depredador 

Por otra parte, se obtuvo la distribución espacial de los agentes y sus densidades a través de este modelo. Podemos observar, a partir de la Figura 3, heterogeneidad, asimetrías y desplazamientos de la posición media, las cuales fueron realizadas a cuatros unidades de tiempos diferentes (60; 100; 160 y 240).

Figura 3: Ejemplos de distribuciones espaciales de agentes presa (en rojo) y depredador (en azul) a diferentes tiempos de evolución: 60; 100; 160 y 240 unidades de tiempo (∆T). Los ejes muestran las coordenadas de los agentes en el hábitat. 

Para caracterizar cuantitativamente estas distribuciones, el código implementado calcula para cada unidad de tiempo, los cuatro primeros momentos (^{Papoulis, 1991}; ^{Cáceres, 2002}; ^{Spiegel y Stephens, 2009}) para cada eje de las distribuciones espaciales de los agentes de ambas especies. Estas distribuciones espaciales, nos permitió discriminar las diferentes situaciones observadas en la Figura 3, tales como vaciamientos de áreas, agrupamientos en los extremos del hábitat y aumentos simétricos o asimétricos de la densidad. La figura 4 muestra la evolución de estos parámetros para la ordenada (eje vertical): ȳ , σ, asimetría y curtosis de la distribución de agentes presa y depredador respectivamente. Se puede notar apartamientos de los cuatro primeros momentos de la distribución espacial respecto de los correspondientes a una distribución espacial homogénea. Estos apartamientos son mayores para la distribución de presas que de depredadores (ver Fig. 3). Por ejemplo, para el tiempo correspondiente a 160, se observa que la distribución de presas presenta un notorio desplazamiento hacia valores positivos en la ordenada (coordenada y), asociado con el máximo que muestra la curva en la figura 4a) mostrando a ȳ = 50 muy apartado del centro del hábitat disponible. En Fig. 4b), la desviación cuadrática media, para este mismo tiempo, cae muy por debajo del valor R_h (valor esperado para una distribución homogénea), mostrando una concentración de la población en ese momento. Las figuras 4c) y 4d) muestran la asimetría alrededor de la media y el tipo de distribución leptocúrtica (Spiegel y Stephens, 2009) que presenta un extremo muy por encima del valor 3 para la curtosis a dicho tiempo.

IV. DISCUSIÓN, PERSPECTIVAS Y CONCLUSIONES

En el presente trabajo se desarrolló un modelo no determinista a través de la herramienta de sistemas multiagentes. Los resultados obtenidos se muestran acordes de manera cualitativa con comportamientos empíricos, tales como los observados en Parque de la Isla Nacional Royale (^{Peterson y Vucetich, 2022}). Es decir, en este estudio de campo se reportan oscilaciones acopladas con desfasaje temporal en la abundancia de alces y lobos. Siguiendo en esta línea, notamos que la escasez de datos experimentales en la bibliografía especializada torna dificultoso la contrastación con nuestro modelo, por lo que se decidió contrastar el modelo multiagente con los resultados del modelo analítico de LV. Este modelo analítico, aunque muy restrictivo respecto a sus suposiciones y el cual se cumple en escasas situaciones reales, puede considerarse un patrón a reproducir cuando estas situaciones se satisfacen. En este trabajo, fue modelado la situación simple, que se ajusta a las condiciones en las cuales LV es confiable. Por lo tanto, la semejanza cuali-cuantitativa de los resultados del modelo desarrollado en la predicción del comportamiento de la cantidad de agentes de ambas especies genera una confianza inicial en el modelo. El valor de los resultados en esta instancia no es tanto sus valores cuantitativos sino su semejanza en el tipo de comportamiento global predicho, que establece oscilaciones y acoplamiento. Es decir, nuestro modelo reproduce características distintivas de LV tales como ciclos, retrasos entre especies y oscilaciones acoplada.

El modelo es versátil y permite el estudio de distribuciones espacio temporales de los agentes. En particular, se puede notar que la distribución espacial de los depredadores no presenta las notorias alteraciones anteriormente descritas en la metodología para las presas. Su distribución espacial aproximadamente homogénea se visualiza en el valor medio de la coordenada cercano a cero, y las curvas de desviación cuadrática media, asimetría y curtosis sin mayores desviaciones (no presenta los extremos encontrados para las presas). Este resultado es razonable dado que los depredadores se mueven al azar por el hábitat sin ser perturbados por otra especie y su velocidad es grande por lo cual en cada ciclo barren gran parte del espacio disponible para sensar presas y nacer o morir, homogeneizando espacialmente a la población de depredadores. Una característica importante de este modelo multiagente es la posibilidad de obtener diferentes evoluciones aun pasando por el mismo punto en el espacio de las fases (para el número de agentes o biomasa de presas/depredador) dado que este modelo sigue la evolución de cada agente y no sólo su número como en el modelo analítico de LV. La distribución espacial de los agentes y las direcciones de viaje cambian las interacciones y la evolución puede por lo tanto ser diferente. Respecto al estudio de la masa de los agentes, podemos observar una evolución con quiebres menos pronunciados en la biomasa de los agentes depredadores, lo que obedece a que los depredadores que desaparecen del hábitat son aquellos con masa menor a la masa crítica, por lo que el cambio de mayor relevancia relativa es en el número de depredadores que en su biomasa.

Por lo tanto, el punto de partida presentado en este modelado multiagente, permitirá estudiar aspectos vinculados a la circulación de la biomasa y, por lo tanto, de la energía en el ecosistema (^{Granda Velepucha et al, 2016}), puede también utilizarse para evaluar la relación de masa presa/predador (^{Uchmanski, 1987}) o presencia de contaminantes en el sistema. En este último caso considerando al contaminante como una especie con tasa de reproducción nula si el contaminante es conservativo o con una tasa de reproducción negativa en caso de que éste presente reacciones secundarias de degradación en el ecosistema (^{Fjeld et al, 2007}).

El modelo propuesto es pasible de complejizarse en múltiples formas, tales como la descripción de la estructura del hábitat (externo a los agentes: topografía, alimento, refugio, zonas de apareamiento, residencia, alimentación) o la estructura de las poblaciones (interna: sexos, edades, habilidades). También existe el camino de agregar otras poblaciones de especies diferentes (productores primarios, otras presas o predadores que compitan por el alimento y/o el espacio). La incorporación de la distribución espacial de los agentes en las interacciones y memorias individuales (costosa en cálculo y almacenamiento) o colectiva, permitiría asignar cierta inteligencia a los agentes de manera que autónomamente puedan sensar o recordar dónde conviene alimentarse (donde encuentren mucha concentración de alimento) o dónde es riesgoso ubicarse (presas, donde interactúen con mucha concentración de predadores). Por último, evaluar probabilidades de ramificaciones mediante el cálculo de probabilidad frecuentista con múltiples corridas a bajas poblaciones, o comportamientos emergentes con foco y detalle en las distribuciones espaciales.