Resumen: En este trabajo se presenta un método numérico para observar la distribución y el comportamiento entre la interacción de presas y depredadores bajo un modelo difusivo bidimensional, con crecimiento logístico para las presas y funcional de depredación tipo Holling II. Al realizar algunas perturbaciones en los parámetros del modelo, determinar condiciones de contornos apropiadas y establecer intervalos de tiempo para la convergencia del método, se presentan diversos patrones sobre las soluciones simuladas en el modelo. En vista que el modelo matemático sin difusión presenta ciclo límite, un equilibrio que puede ser localmente un nodo o una espiral estable, las soluciones numéricas del modelo difusivo reflejan dichos comportamientos.
Palabras clave: Método de diferencias finitas, operador laplaciano, bifurcación de Hopf, coeficiente de difusión.
Abstract: This paper presents a numerical method to observe the distribution and behavior between prey and predator interaction under a two-dimensional diffusive model, with logistic growth for prey and Holling II type predation functional. By performing some perturbations in the model parameters, determining appropriate boundary conditions and establishing time intervals for the convergence of the method, several patterns are presented on the simulated solutions in the model. Since the mathematical model without diffusion presents a limit cycle, an equilibrium that can be locally a node or a stable spiral, the numerical solutions of the diffusive model reflect these behaviors.
Keywords: Finite difference method, Laplacian operator, Hopf bifurcation, diffusion coefficient.
Artículos de revisión
Formación de patrones en un modelo difusivo bidimensional depredador - presa tipo Holling II
Pattern formation in a two-dimensional diffusive predator model - Holling II type prey
Recepción: 08 Marzo 2021
Aprobación: 29 Septiembre 2021
Autor de correspondencia: cc.cortes@cnb.csic.es
La formulación matemática que describe un fenómeno biológico determinista puede ser descrito por un sistema dinámico, el cual permite prever cualquier estado futuro si el operador de evolución y su estado inicial son conocidos. Si el tiempo para ese operador se define en la recta real, se habla de un sistema dinámico continuo.
En particular, las ecuaciones diferenciales de la forma (Ecuación 1):
Con la variable de estados y un vector de parámetros es un sistema dinámico y del cual pueden ser usados para describir y analizar todas las posibles dinámica de un fenómeno.
Los retratos de fases para el modelo (1) del cual cambia drásticamente cuando se altera el parámetro , se denomina bifurcación. Los trabajos publicados por Hirsch y Smale (1974), Sotomayor (1979), Arrowsmith (1992), Berkooz, et al., (1996), Kuznetsov (1998), Blanchard y Devaney (1999), Chicone (2006) y, Renardy y Rogers (2006) y Cortés (2017) describen la forma matemática para analizar cualitativamente el modelo (1) respecto a sus equilibrios, además de todos los posibles casos de bifurcaciones, tanto local como global.
Los modelos biológicos son planteados para analizar la dinámica en la interacción de dos especies bajo un sistema de ecuaciones diferenciales ordinarias, conocidos generalmente como modelo depredador – presa. Un ejemplo para este tipo de modelos es expresado por (Ecuación 2):
En los trabajos publicados por Freedman (1980), Edelstein (1988), Kuang y Freedman (1988), Wolkowicz (1988), Murray (2002), Zhou (2019), Zhu, et al., (2003), Du y Shi (2006), Liu (2010), Hasik (2010), Zu y Mimura (2010), Tian y Xu (2011), Tang, et al., (2014), Lin, et al., (2018), Mukherjee y Maji (2020) y, Cortés (2021a) se plantean y se analizan cualitativamente algunos modelos biológicos de la forma (Ecuación 1) que describen el comportamiento de las presas y depredadores bajo ciertas suposiciones.
Sin embargo, al incluir algunos elementos matemáticos, como la difusión unidimensional de la especie, el modelo (1) se replantea como un sistema de ecuaciones diferenciales parciales de la forma (Ecuación 3):
Con D el vector de difusión, con entradas no negativas, del cual permite analizar de diversas maneras el comportamiento y la dinámica de dicha población de una manera más acorde con el ambiente. En particular, al agregar un término difusivo unidimensional a las presas y/o depredadores al modelo (2) se tiene que (Ecuación 4):
con y las constantes de difusión para las presas y depredadores, respectivamente.
El estudio analítico para modelos biológicos de la forma (Ecuación 3) están centrados en mostrar soluciones tipo ondas viajeras de la forma , con la velocidad de onda. En los trabajos publicados por Grindrod (1991), Huang (2003), Li y Wu (2008), Hong y Weng (2013), Wu, et al., (2014), Ghazaryan, et al., (2015), Ai, et al., (2017) y, Cortés y Ramírez (2021b), y muestran mecanismos matemáticos para demostrar la existencia de ondas viajeras como soluciones al modelo de la forma (3).
Sin embargo, el análisis matemático hace cada vez más complejo si se al incluir un término difusivo bidimensional al modelo (1), esto es, (Ecuación 5):
Con un operador Laplaciano de W. En particular, si el modelo (2) es modificado por (Ecuación 6):
Su estudio se centra en determinar todos los posibles patrones que presentan dichas soluciones de y en el plano xy para un tiempo definido. Los trabajos publicados por Cohen y Grossberg (1983), Murray y Oster (1984), Doelman (1995), Golubitsky, et al., (1999), Cruywagen, et al., (2000), Hoyle y Hoyle (2006), Mancini (2007), Cross y Greenside (2009), Cárdenas, et al., (2013), Rao y Kang (2016) y, Wu, et al., (2018) y, muestran diversas maneras para encontrar patrones a modelos biológicos de la forma (Ecuación 5).
Para encontrar una solución explicita de un modelo biológico, y en general para cualquier modelo matemático compuesto por ecuaciones diferenciales parciales, no siempre es fácil. De esa forma, se buscan soluciones aproximadas mediante la aplicación de métodos computacionales apropiados para resolver problemas matemáticos difíciles de resolver o que no se pueden resolver analíticamente. En este caso, existen diversos métodos numéricos para aproximar soluciones en modelos de ecuaciones diferenciales y que podrían ser consultados en los trabajos de Mathews, et al., (2000), Iserles (2009) y Sauser (2013).
Por tanto, en este trabajo se determinan los posibles patrones, al ser uso de un método numérico, para un modelo depredador - presa tipo Holling II (6) y término difusivo bidimensional. Para ello, al tomar los trabajos publicados por Garvie (2005ab, 2007), en la segunda sesión se presenta un resumen del método numérico por diferencias finitas, ejecutado por Garvie, para resolver numéricamente un modelo de ecuaciones diferenciales parciales con término difusivo bidimensional.
Por otro lado, y a partir del trabajo publicado por Cortés (2021b), en la tercera sesión se describe un modelo depredador-presa tipo Holling II con termino difusivo bidimensional y, en la cuarta sesión, se presenta un análisis cualitativo al modelo cuando los parámetros de difusión son nulos.
Por último, al usar los lineamientos que se presentaron en las secciones anteriores y al variar algunos parámetros y condiciones iniciales, en la quinta sesión se determinan una serie de patrones que muestran el comportamiento del modelo difusivo bidimensional depredador-presa con respuesta funcional tipo Holling II y se realiza una comparación con la dinámica del modelo sin difusión.
Para generar soluciones numéricas al sistema (Ecuación 7):
Con y conocido como un operador Laplaciano de z, se emplea un método numérico de diferencia finita, el cual sustituye las derivadas parciales de una ecuación diferencial por aproximaciones discretas y graficar las soluciones numéricas del modelo (7) sobre una malla con una programación adecuada de sus nodos
Según Garvie (2007), para discretizar el modelo (2) en el intervalo de tiempo y niveles de tiempo , para todo , con , considere una malla cuadrática de puntos dada por , con , y una subdivisión uniforme de puntos en la cuadrícula , con incremento de espacio .
Como lo observado en la Figura 1 y al considerar , los nodos rellenos representan los valores conocidos de la solución aproximada, a partir de las condiciones iniciales y de fronteras, y los nodos ficticios son los puntos de la malla que se deben programar mediante el método, esto es, incógnitas para resolver, donde representa la solución aproximada de del modelo (7) en el punto .
Los nodos en la malla están numerados naturalmente, de izquierda a derecha, y comienza con la fila inferior, desde k=0,1,..., J. La relación entre la numeración de los nodos y la indexación (i,j) está dada por , donde k= i +j (J+1) para i, j = 0, ..., J.
Una forma discreta del modelo (7), con para n = 1,..., N e i, j = 0,...J, es dado por (Ecuación 8):
donde ^f y ^g son funciones discretas, en términos del tiempo anterior tn-1 , modificadas de f y g, esto es (Ecuación 9),
y las aproximaciones de y del operador Laplaciano discreto en dos dimensiones son representadas por (Ecuación 10):
A partir de las condiciones iniciales del modelo (Ecuación 7), dadas por , y las condiciones de contorno a lo largo de los cuatro bordes del cuadrado, es decir, para , se tiene que (Ecuación 11),
y de forma análogamente para las aproximaciones de v.
Como lo observado en la Figura 1, las condiciones de frontera de las esquinas para u son dadas por (Ecuación 12):
donde se obtiene de igual forma para las aproximaciones de v.
Para calcular las 2(J+1) ecuaciones lineales que representan los nodos interiores, se usa la siguiente matriz de bloque, con la numeración natural de los nodos (Ecuación 13),
con y para k=0,...,J.
La matriz de coeficientes constantes B1 y B2 son calculadas por (Ecuación 14):
donde la matriz L, de dimensión (J+1)2 x (J+1)2 y bloques (J+1) x (J+1), es dada por (Ecuación 15):
con
Sean y la población de presas y depredadores, respectivamente, cuya distribución bidimensional x, y siguen un modelo de la siguiente forma (Ecuación 16):
donde
son los coeficientes de difusión de las presas y depredadores, respectivamente,
describe el crecimiento logístico de las presas en ausencia de los depredadores, con la tasa de crecimiento y K >0 la capacidad de carga,
es la disminución en las presas al interactuar con los depredadores, donde es la tasa de encuentro entre las dos especies y e>0 la constante de saturación de las presas,
describe el aumento en la población de los depredadores al interactuar con las presas, donde m>0 es la tasa de encuentro entre las dos especies,
es el decrecimiento de los depredadores en ausencia de las presas con d>0 la tasa de muerte de los depredadores.
Nótese que el modelo (16) es equivalente al modelo (7) con y
En el caso que los coeficientes de difusión en el modelo (16) sean nulos, es decir D1=D2=0, el modelo (16) toma la forma de un sistema de ecuaciones diferenciales ordinarias, esto es (Ecuación 17),
el cual describe la cantidad de presas y depredadores a lo largo del tiempo.
El modelo (17) tiene dos puntos de equilibrio sobre los ejes xy dados por: E1= (0,0) y E2=K,0. Además, si m <de y d(eK+1), el modelo (11) presenta un punto de equilibrio interior dado por
Para determinar la estabilidad local de cada equilibrio sobre el modelo (17) se debe analizar el signo de la traza y el determinante de la matriz jacobiana evaluada en cada equilibrio. En efecto, la matriz Jacobiana del modelo (17) evaluada en el punto de equilibrio E1 es dada por (Ecuación 18).
con Det A(E1) = -rd<0. Por consiguiente, el punto E1 es localmente un punto silla.
De igual forma, la matriz Jacobiana del modelo (17) asociada al punto E2 esta dada por (Ecuación 19):
con valores propios: y . Si entonces y así el punto E2es localmente un nodo estable, y para se tiene , esto es, el punto E2 es localmente un punto silla.
Proposición 1. El modelo (12) tiene dos puntos de equilibrio E1=(0,0) y E2= (K,0) y un punto de equilibrio interior, si m>de y mK>d(eK+1) dado por
Las condiciones de estabilidad local para los puntos de equilibrios están dadas de la siguiente manera:
El punto E1= (0,0) es locamente un punto silla.
El punto E2=(K,0) es locamente
una silla si ,
un nodo estable si .
El punto E* es localmente
una espiral inestable si y ,
una espiral estable si y,
un nodo estable si y ,
un nodo inestable si y .
Además, si el punto de equilibrio E* interior es localmente inestable, entonces el modelo (6) presenta un único ciclo límite globalmente estable en el primer cuadrante.
Por otro lado, el modelo (17) presenta una bifurcación de Hopf al colisionar el equilibrio E*con el ciclo limite estable, esto es, si Tr A=0 , y cuya prueba puede ser vista en Cortés, et al., (2021).
Lema 1. Si entonces el modelo (17) presenta una bifurcación de Hopf.
De igual forma, se deduce que la colisión del equilibrio E* con el punto de silla E2 genera un cambio en la estabilidad de E2. En este caso, se forma una bifurcación transcrítica si la componente sobre el eje x del equlibrio E* concuerda con la componente sobre el eje x de E2 cómo se observa en el siguiente resultado.
Lema 2. Si entonces el modelo (17) presenta una bifurcación transcrítica.
En la Figura 2 se muestran las curvas de bifurcación del modelo (17) con sus retratos de fase. En la región 1 se observa que el punto E2 es un nodo estable que se transforma en un punto silla a medida que K cruza la curva y posteriormente se forma un punto de equilibrio interior estable dado en la región 2.
En la región 3, al traspasar los parámetros la curva , el nodo estable se transforma en un foco estable como lo observado en la región 2, y a medida que traspasa la curva , el modelo (17) presenta un ciclo límite estable como se muestra en la región 4.
Al aplicar el método numérico, dado en la sección uno, al modelo (16), se procede a realizar simulaciones para detectar posibles patrones en la forma como se distribuyen las especies sobre una malla. El código de programación fue ejecutado en Matlab 2016.
Si D1=1, D2=10, r=1, K=1, b=2.5, e=2.5, m=5 y d=0.6, se observa en la Figura 3 que tanto los depredadores como las presas se distribuyen en forma de espiral a lo largo de la malla [0,500]x[0,500], y se van desvaneciendo a medida que aumenta el valor de T, formando diversas manchas. Este fenómeno concuerda con la dinámica del modelo (17) presentado en la Figura 2.
Al usar los mismos parámetros que generan la Figura 3 con otras condiciones iniciales u(x,y,0) y v(x,y,0), se observa en la Figura 4 unos cambios drásticos en la distribución de las presas y depredadores para T pequeño, el cual, a diferencia de la Figura 3, los depredadores como las presas se distribuyen en forma de nodos a lo largo de la malla [0,500]x[0,500].
Para los parámetros y condiciones dadas en la Figura 5 solo presenta un único patrón. En dicha Figura se puede observar que la distribución de ambas especies presenta ciclos, lo cual concuerda con la dinámica del modelo (11) mostrado en la Figura 2.
El método numérico de diferencias finitas resultó ser efectivo para el desarrollo de la simulación, de tal forma que se logró simular las distintas formas de distribución y patrones entre presas y depredadores, cuyo comportamiento es descrito por el modelo difusivo (16). Sin embargo, el método no permitió detectar los posibles comportamientos posteriores al ciclo formado. Para analizar dichos comportamientos cuando se aumenta el valor de T, se podría aplicar otros métodos numéricos, como el método de elementos finitos.
La formación de patrones dados en el modelo (16) concuerda con las presentadas en la dinámica del modelo sin difusión (17), por lo que se podría pensar que dicho modelo difusivo presenta una bifurcación de Hopf.
cc.cortes@cnb.csic.es