© Revista de Matemática: Teoría y Aplicaciones 2022 29(1) : 105–138 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v29i1. 48408

influencia del efecto allee débil en las presas en un modelo de depredación del

tipo leslie-gower con respuesta funcional sigmoidea

influence of the weak allee effect on prey in a leslie-gower type predation model with sigmoid functional response

Sebastián Valenzuela-Figueroa^{∗}

Eduardo González-Olivares^{†}  Alejandro Rojas-Palma^{‡}

Received: 10/Sep/2021; Revised: 6/Dec/2021; Accepted: 10/Dec/2021

∗Universidad Austral de Chile, Centro de Docencia Superior en Ciencias Básicas, Puerto Montt, Chile. E-Mail: sebastian.valenzuela@uach.cl

†Pontificia Universidad Católica de Valparaíso, Instituto de Matemáticas, Valparaíso, Chile. E-Mail: ejgonzal@ucv.cl

‡Universidad Católica del Maule, Departamento de Matemática, Física y Estadística, Facultad de Ciencias Básicas, Talca, Chile. E-Mail: amrojas@ucm.cl

106 s. valenzuela-figueroa et al.

Resumen

Utilizando un sistema topológicamente equivalente al original, dependiente sólo de cuatro parámetros, en este trabajo analizamos un modelo de depredación del tipo Leslie-Gower, considerando que el consumo de los depredadores es modelado por una respuesta funcional sigmoidea. Además, asumimos que las presas están afectadas por un efecto Allee y que los depredadores son generalistas. Mostramos que el sistema de ecuaciones diferenciales ordinarias que describe el modelo puede tener hasta cuatro puntos de equilibrio positivos. Dadas las dificultades para obtener explícitamente las coordenadas de estos puntos analizamos parcialmente el sistema considerando que la pobla- ción de presas está afectada por un efecto Allee débil. Entre los resultados más importantes obtenidos, se demuestra la existencia de una curva separatriz, dividiendo el comportamiento de las soluciones o trayectorias del sistema en el plano de fase. Dos soluciones muy cercanas pero en un lado diferente de esa separatriz, tendrían ω − limites diferentes y distantes. Esto implica que teniendo un mismo tamaño poblacional de las presas, para distintos tamaños poblacionales de depredadores, pero muy cercanos, ambas poblaciones podrían coexistir o las presas podrían ir a la extinción.

Palabras clave: modelo depredador-presa; bifurcación; ciclo límite; curva separatriz; estabilidad; respuesta funcional.

Abstract

Using a topologically equivalent system to the original, dependent only on four parameters, in this work we analyze a Leslie-Gower type predation model, considering that predator consumption is modeled by a sigmoid functional response. Furthermore, we assume that the prey are affected by an Allee effect and that the predators are generalists. We show that the system of ordinary differential equations (ODE) that the model describes can have up to four positive equilibrium points. Given the difficulties in obtaining explicitly the coordinates of these points, we partially analyze the system considering that the prey population is affected by a weak Allee effect. Among the most important results obtained, the existence of a separator curve is demonstrated, dividing the behavior of the solutions or trajectories of the system in the phase plane. Two very close solutions, but on a different side of that separatrix, would have different and distant ω − limit sets. This implies that having the same population size of the prey, for different population sizes of predators, but very close, both populations could coexist or the prey could go to extinction. We estimate that the analytical results obtained have an adequate ecological interpretation, under the underlying assumptions in the modeling of the predation interaction with ODEs.

Keywords: predator-prey model; bifurcation; limit cycle; separatrix curve; stability; functional response.

Mathematics Subject Classification: 92D25; 34C23; 58F14; 58F21.

1 Introducción

Es bien sabido que la principal característica de los modelos depredador-presa del tipo Leslie [28], o modelo logistico [7, 45], o modelos del tipo Leslie-Gower [1, 2, 21, 35], es que la capacidad de carga del medio ambiente de los depredadores K_{y} es una función del tamaño de la población de presas x = x (t), existente en el medio ambiente donde se realiza la interacción, es decir, la capacidad de carga de los consumidores, depende de los recursos alimenticios disponibles [45].

En el modelo original propuesto por el ecólogo británico Patrick Holt Leslie en 1948 [3], como alternativa al cuestionado modelo de Lotka-Volterra [3, 4] formulado en 1926, se asume que la capacidad de carga ambiental de los depredadores es proporcional a la abundancia de presas, es decir, K_{y} = K (x) = nx.

En este trabajo analizamos un modelo de depredación del tipo Leslie-Gower [28, 29], considerando tres aspectos importantes:

i. La acción de los depredadores o respuesta funcional es una función sigmoidea racional [20, 22, 41],

ii. los depredadores cuentan con un alimento alternativo [2, 23] y,

iii. las presas están afectadas por un efecto Allee [12].

Efecto Allee

El efecto Allee es un fenómeno que afecta a algunas poblaciones [40], el cual se caracteriza porque existe una relación positiva entre la tasa de crecimiento y el tamaño de la población en densidades bajas [40], principalmente debido a la cooperación entre los individuos (comportamiento social) [12]. A altas densidades de la misma población, la competencia es más dominante que la cooperación, y en consecuencia el efecto Allee pierde su impacto. Esto implica que la relación se invierte a altas densidades [19].

Este efecto recibe diferentes nombres [30] entre los más usados es depensación [9]. Es importante para comprender por qué algunas poblaciones se enfrentan a la extinción a bajas densidades de población [12, 40].

Es originado por diversas causas [6, 12] y su estudio se ha incrementado en las últimas décadas, debido a la alta probabilidad que las poblaciones afectadas puedan extinguirse [8]. Puede ser descrito mediante distintas expresiones algebraicas, muchas de las cuales son topológicamente equivalentes [19].

En este trabajo se usará la forma matemática más común en la literatura ecológica [19, 20, 23], representada por la ecuación polinomial siguiente:

\frac{dx}{dt} = r (1 − \frac{x}{K} ) (x − m) x.

De acuerdo a la definición de este efecto [12, 40], se debe cumplir que \frac{dx}{dt} > 0 para un cierto rango de valores del parámetro m; claramente esto se cumple si −K < m << K. Cuando m > 0, se tiene un efecto Allee fuerte; existe un efecto Allee débil si m < 0 y un efecto Allee débil especial cuando m = 0 [19, 20, 23]. En Ciencias Pesqueras, los mismos fenómenos se denominan depensación crítica y depensación pura, respectivamente [9, 30].

Además, el efecto Allee puede ser producido por diferentes causas; diversos mecanismos generados por distintas fuentes pueden actuar simultáneamente en una misma población (ver la Tabla 2 en [6] o la Tabla 2.2 en [12]). La combinación de estos fenómenos se denomina efecto Allee múltiple [6].

Una ecuación que describe un efecto Allee doble es dada por:

\frac{dx}{dt} = ( 1 − \frac{x}{K} − \frac{r}{x + n} ) x,

la cual hemos denominado efecto Allee aditivo [1], por la manera de obtenerse [12, 40]. Si se incorpora en el modelo de Leslie-Gower, esta ecuación racional podría generar dinámicas más complejas que las estudiadas en este trabajo, particularmente en la cantidad de ciclos límites [1].

Respuesta funcional

Por otro lado, la respuesta funcional del depredador o función de consumo se refiere al cambio en la densidad de presas atacadas por unidad de tiempo por depredador, cuando cambia la densidad de presas [15], hasta que el tiempo de manipulación comienza a limitar su consumo, pues la curva es acotada [45].

En 1959 el biólogo canadiense Crawford. S. Holling [26], basado en experimentos realizados en laboratorio, describió tres tipos de funciones saturadas, dependiendo solo del tamaño de la población de presas (respuesta funcional presadependiente). Se denominan tipo Holling, I, II o III [45]; más tarde, Robert J. Taylor [42] en 1984 propuso la respuesta funcional Holling tipo IV o no monótonica [42], para describir el fenómeno llamado formación de grupo de defensa.

También han sido formuladas respuestas funcionales dependientes de ambas poblaciones, como son las razón-dependientes (ratio-dependent), de Beddington-DeAngelis, de Hassel-Varley, o de Crowley-Martin [45].

Las respuesta funcionales de Holling tipo III describen el comportamiento depredatorio de especies consumidoras que son generalistas [45]. Estas especies en ausencia o escasez de sus presas favoritas opta por capturar o consumir un alimento alternativo [20, 22, 41]. Gráficamente una respuesta Holling tipo III se representa mediante una curva sigmoidea. [46].

En este trabajo consideramos una respuesta funcional de Holling tipo III o sigmoidea [20, 22] o de forma de S [41], descrita por la función h(x) = \frac{q x^{2}}{x^{2} + a^{2}} , con q y a > 0. En la literatura ecológica, una forma general para describir una respuesta funcional sigmoidea es h (x) = \frac{qx^{n}}{x^{n}+a^{n}}, con q, a > 0 y n ∈ ℕ [41].

Uno de los principales objetivos en el análisis de modelos de depredación es establecer el número de ciclos límites que rodean un punto de equilibrio positivo del tipo centro-foco [11], originados por bifurcaciones de Hopf [10, 16, 17, 27].

En este trabajo, será demostrada la existencia de al menos un ciclo límite en el modelo, el cual es obtenido por bifurcación de Hopf [10, 16, 34]; este ciclo queda en evidencia mediante una simulación numérica.

En la Tabla reftab1, mostraremos un resumen de varios modelos en que la acción de los depredadores es descrita por una respuesta funcional sigmodea, considerando que la población de presas está afectada o no por el efecto Allee. En esta Tabla se puede obsevar que la cantidad de ciclos límites cambia, por causa de este fenómeno o por la expresión algebraica de la respuesta funcional.

Este artículo está organizado del siguiente modo: en la Sección 2, presentamos el modelo y, por medio de un difeomorfismo, se obtiene un modelo topológicamente equivalente al modelo depredador-presa del tipo Leslie-Gower; en la Sección 3, se presentan las principales propiedades del sistema. Simulaciones numéricas se muestran en la Sección 4 y se incluye una Tabla comparando modelos en que se consideran respuestas funcionales sigmoideas. Algunas consecuencias ecológicas y una discusión de los resultados matemáticos se dan en la Sección 5.

2 Proposición del modelo

En este trabajo se analiza una modificación al modelo formulado por Patrick H. Leslie en 1948 [28]. Es representado por un sistema autónomo bidimensional de ecuaciones diferenciales ordinarias (EDO)no lineales del tipo Kolmogorov [15, 18] dado por:

X_{μ} (x, y) : \{ \begin{matrix}\frac{dx}{dt} & = & r (1 − \frac{x}{K} ) (x − m) x − \frac{q x^{2}}{x^{2} + a^{2}} y \\ \frac{dy}{dt} & = & s (1 − \frac{y}{nx+c}) y, \end{matrix}  (1)

donde x = x (t) e y = y (t) son los tamaños poblacionales de las presas y los depredadores respectivamente, para t ≥ 0, medido como biomasa o densidad por unidad de área o volumen, con μ = (r, K, q, a, s, n, m, c) ∈ ℝ_{+}^{6} ×]−K, K[ × ]0, ∞[. Los parámetros tienen los siguientes significados ecológicos:

• r y s indican la tasa de crecimiento intrínseca de la población de presas y depredadores, respectivamente,

• K es la capacidad de carga del medio ambiente para las presas,

• m es el parámetro relacionado al efecto Allee, m > 0, es el umbral de extinción del efecto Allee que afecta a la presa o al mínimo de población viable [12, 21],

• q es la tasa máxima de consumo per capita de los depredadores,

• a es la cantidad de presas necesarias para alcanzar la mitad de q, o tasa de saturación media,

• n representa la calidad energética que proporcionan las presas como alimento para los depredadores,

• c indica el tamaño máximo disponible del alimento alternativo.

Por conveniencia, anotamos a^{2} la tasa de saturación media, pero que también se puede escribir a. El sistema (1) está definido en:

Ω = {(x, y) ∈ ℝ^{2}| x ≥ 0, y ≥ 0} = ℝ_{0}^{+} × ℝ_{0}^{+}.

Equivalencia topológica

Con el objeto de simplificar los cálculos, seguimos la metodología utilizada en [1, 2, 22, 31], efectuando un cambio de variables dado por x = Ku e y = nKv, y un reescalamiento del tiempo dado por t = \frac{(u+C)(u^{2}+A^{2})}{rK}τ [25, 32, 35, 39]. Después de un trabajo algebraico obtenemos un nuevo sistema polinomial de EDO, también del tipo Kolmogorov [15, 18], descrito por:

Y_{η} (u, v) : \{ \begin{matrix}\frac{du}{dτ} & = & u(u + C) ((1 − u)(u − M) (u^{2} + A^{2}) − Quv) \\ \frac{dv}{dτ} & = & Bv(u + C − v)(u^{2} + A^{2}) \end{matrix} (2)

con A^{2} = \frac{a^{2}}{K}, B = \frac{s}{Kr}, C = \frac{c}{nK}, Q = \frac{q}{nr} y M = \frac{m}{K}, donde η = (A, B, C, Q, M) ∈ ]0, 1[ × (ℝ_{0}^{+})^{3} ×]−1, 1[ y está definido en

Ω̄ = {(u, v) ∈ ℝ^{2}| u ≥ 0, v ≥ 0} .

Nota 1 Los sistemas (1) y (2) son toplógicamente equivalentes porque se ha construido el difeomorfismo [10, 14] φ : Ω̄ × ℝ ⟶ Ω × ℝ, tal que

φ(u, v, τ) = (Ku, nKv, \frac{(u+C)(u^{2}+A^{2})}{rK}τ) = (x, y, t) .

El determinante de la matriz Jacobiana de la función φ es:

det Dφ(u, v, τ) = det ( \begin{matrix}K & 0 & 0 \\ 0 & nK & 0 \\ \frac{1}{rK} & 0 & \frac{1}{rK} (u + C) (u^{2} + A^{2}) \end{matrix} )
= \frac{nK}{r} (u + C) (u^{2} + A^{2}) > 0.

Entonces, φ es un diffeomorfismo que preserva la orientación del tiempo. El sistema (2) puede considerarse una extensión continua del sistema (1)

definido en Ω̄ .

En esta subsección asumiremos primero que M > 0, pero más adelante consideraremos el caso M = 0.

Los puntos de equilibrio del sistema (2) o campo vectorial Y_{η} (u, v) son (0, 0), (1, 0), (0, C), (M, 0) cuando M > 0, y los puntos (u_{e}, v_{e}) que están en la intersección de las isoclinas:

v = u + C

v = \frac{(1 − u)(u − M)(u^{2} + A^{2})}{Qu}.

Se obtiene la ecuación polinomial

\begin{matrix}p (u) = u^{4} − (1 + M) u^{3} + (A^{2} + M + Q) u^{2} \\ − (A^{2} (M + 1) − CQ) u + MA^{2} = 0. \end{matrix}  (3)

Por lo tanto, tenemos los siguientes casos: el polinomio p (u) puede cambiar:

a) cuatro veces de signo, si y sólo si, (A^{2} (M + 1) − CQ) > 0.

b) dos veces de signo, si y sólo si, (A^{2} (M + 1) − CQ) ≤ 0.

Luego, para analizar la dinámica del sistema (2), y por consecuencia las del sistema (1) se deben considerar diferentes casos en los cuales es difícil determinar explícitamente las coordenadas de los puntos de equilibrio positivos (o al interior del primer cuadrante).

Notamos que si se considera el efecto Allee débil especial cuando M = 0, la ecuación polinomial (3) se reduce a:

p_{0} (u) = (u^{3} − u^{2} + (A^{2} + Q) u − (A^{2} − CQ)) u = 0  (4)

Una de las soluciones es cero, y las otras son las raíces de la ecuación polinomial

p_{01} (u) = u^{3} − u^{2} + (A^{2} + Q) u − (A^{2} − CQ) = 0  (5)

la cual puede tener:

I) Tres raíces reales positivas, o una, si A^{2} − CQ > 0.

II) Dos raíces reales positivas, o una de multiplicidad 2, o ninguna, si A^{2} − CQ ≤ 0.

Si A^{2} − CQ = 0, se tiene la ecuación polinomial

p_{010} (u) = (u^{2} − u + (A^{2} + Q)) u = 0

una de cuyas raíces es u = 0 y las otras dependen de Δ_{0} = 1 − 4 (A^{2} + Q).

Además, como Q = \frac{A^{2}}{C}, se puede reducir la cantidad de parámetros, obteniendo:

p_{010} (u) = u^{2} − u + A^{2} (1 + \frac{1}{C} ) = 0,  (6)

y Δ_{0} = \frac{1}{C} (C − 4A^{2} (C + 1)).

IIa) Si Δ_{0} > 0 se tienen

u_{10} = \frac{1 − \sqrt{1 − 4 (A^{2} + Q)}}{2} y u_{20} = \frac{1 + \sqrt{1 − 4 (A^{2} + Q)}}{2}

con 0 < u_{10} < u_{20} < 1.

IIb) Si Δ_{0} = 0 se tiene u_{10} = u_{20} = \frac{1}{2} = u_{∗0}.

IIc) Si Δ_{0} < 0 no hay raices reales positivas.

Por conveniencia para obtener la matriz Jacobiana del sistema (2), anotaremos:

H = (1 − u)(u − M) (u^{2} + A^{2}) − Quv

Es fácil ver que la matriz Jacobiana del sistema (2) es:

DY_{η} (u, v) = ( \begin{matrix}DY_{η} (u, v)_{11} & −Qu^{2} (C + u) \\ Bv (2Cu − 2uv + A^{2} + 3u^{2}) & B (A^{2} + u^{2}) (C + u − 2v) \end{matrix} )

con

DY_{η} (u, v)_{11} = \frac{∂}{∂u} (u(u + C)H) = (u + C)H + uH + u(u + C)\frac{∂}{∂u} (H)

siendo

\frac{∂}{∂u} (H) = −4u^{3} + 3 (M + 1) u^{2} − 2 (A^{2} + M) u + (A^{2} (M + 1) − Qv) .

3 Resultados principales

Las siguientes propiedades del sistema (2) son independientes del valor del parámetro M.

Lema 1 El conjunto Γ̄ = {(u, v) ∈ Ω̄ : 0 ≤ u ≤ 1, v ≥ 0} es una región positivamente invariante.

Demostración. Como el sistema (2) es del tipo Kolmogorov [15, 18], entonces los ejes coordenados son conjuntos invariantes.

Asumiendo u = 1, se tiene que

\frac{du}{dτ} = (1 + C) (−Qv) < 0.

Cualquiera que sea el signo de \frac{dv}{dτ} < 0, entonces las trayectorias entran a la region Γ̄.

Notamos que el conjunto Γ = {(x, y) ∈ Ω/0 ≤ x ≤ K, y ≥ 0} es una región positivamente invariante del sistema (1).

Lema 2 Las soluciones son uniformemente acotadas.

Demostración. Para demostrar el acotamiento de las soluciones usaremos el método de Compactificación de Poincaré. Definamos el cambio de variables por zv = 1, wv = u, entonces \frac{dz}{dτ} = \frac{dv}{dτ} \frac{−1}{v^{2}} y \frac{dw}{dτ} = \frac{1}{v} (\frac{du}{dτ} − w \frac{dv}{dτ} ).

Sustituyendo en el sistema (1), se obtiene el nuevo sistema dado por:

U_{η} : \{ \begin{matrix}\frac{dw}{dτ} = −\frac{w}{z^{5}} ( \begin{matrix}−w^{4}z + w^{5} − A^{2}w^{2}z^{3} + A^{2}w^{3}z^{2} + Cw^{4}z − Mw^{4}z \\ −A^{2}Bz^{4} − Bw^{2}z^{2} + Bw^{3}z^{2} − Cw^{3}z^{2} + Mw^{3}z^{2} \\ +Qw^{2}z^{2} + A^{2}Cw^{2}z^{3} − A^{2}Mw^{2}z^{3} + CQwz^{3} \\ +A^{2}BCz^{5} + A^{2}CMz^{5} + A^{2}Bwz^{4} − A^{2}Cwz^{4} \\ +BCw^{2}z^{3} + A^{2}Mwz^{4} + CMw^{2}z^{3} − CMw^{3}z^{2} \\ −A^{2}CMwz^{4} \end{matrix} ) \\ \frac{dz}{dτ} = −\frac{B(A^{2}z^{2} + w^{2})(w + Cz − 1)}{z^{2}}. \end{matrix}

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 29(1): 105–138, Jan–Jun 2022

Notemos que el punto de equilibrio (0, 0) del campo vectorial U_{η} es equivalente al punto (0, ∞) del sistema (2).

Reescalando el tiempo dado por τ = z^{5}T se obtiene el sistema:

Ū_{η} : \{ \begin{matrix}\frac{dw}{dT} = −w ( \begin{matrix}−w^{4}z + w^{5} − A^{2}w^{2}z^{3} + A^{2}w^{3}z^{2} + Cw^{4}z − Mw^{4}z \\ −A^{2}Bz^{4} − Bw^{2}z^{2} + Bw^{3}z^{2} − Cw^{3}z^{2} + Mw^{3}z^{2} \\ +Qw^{2}z^{2} + A^{2}Cw^{2}z^{3} − A^{2}Mw^{2}z^{3} + CQwz^{3} \\ +A^{2}BCz^{5} + A^{2}CMz^{5} + A^{2}Bwz^{4} − A^{2}Cwz^{4} \\ +BCw^{2}z^{3} + A^{2}Mwz^{4} + CMw^{2}z^{3} − CMw^{3}z^{2} \\ −A^{2}CMwz^{4} \end{matrix} ) \\ \frac{dz}{dT} = −z^{3}B (A^{2}z^{2} + w^{2}) (w + Cz − 1) . \end{matrix}

La matriz Jacobiana evaluada en el punto (0, 0) del campo vectorial Ū_{η} es

DŪ_{η} (0, 0) = ( \begin{matrix}0 & 0 \\ 0 & 0 \end{matrix} )

es decir, el punto (0, 0) es no hiperbólico.

Para la desingularización del punto (0, 0) utilizaremos la técnica del blowingup direccional [14], en este caso en la dirección w (horizontal). Definamos:

ϕ_{1}^{w} (a, b) = (a, ab) = (w, z)

y sustituyendo en el último sistema obtenemos:

V_{η} : \{ \begin{matrix}\frac{da}{dT} = −a^{5} ( \begin{matrix}a − Bb^{2} + Qb^{2} − ab + CQb^{3} + Bab^{2} − Cab^{2} \\ +Mab^{2} − A^{2}Bb^{4} + A^{2}ab^{2} − A^{2}ab^{3} + Cab \\ −Mab + BCab^{3} − CMab^{2} + CMab^{3} + A^{2}Bab^{4} \\ +A^{2}Cab^{3} − A^{2}Cab^{4} − A^{2}Mab^{3} + A^{2}Mab^{4} \\ +A^{2}BCab^{5} − A^{2}CMab^{4} + A^{2}CMab^{5} \end{matrix} ) \\ \frac{db}{dT} = a^{4}b (Cb + 1) ( \begin{matrix}a + Qb^{2} − ab + Mab^{2} + A^{2}ab^{2} − A^{2}ab^{3} \\ −Mab − A^{2}Mab^{3} + A^{2}Mab^{4} \end{matrix} ) . \end{matrix}

Reescalando el tiempo considerando α = a^{4}T nos resulta el sistema:

V̄_{η} : \{ \begin{matrix}\frac{da}{dα} & = & −a ( \begin{matrix}a − Bb^{2} + Qb^{2} − ab + CQb^{3} + Bab^{2} − Cab^{2} \\ +Mab^{2} − A^{2}Bb^{4} + A^{2}ab^{2} − A^{2}ab^{3} + Cab \\ −Mab + BCab^{3} − CMab^{2} + CMab^{3} + A^{2}Bab^{4} \\ +A^{2}Cab^{3} − A^{2}Cab − A^{2}Mab^{3} + A^{2}Mab^{4} \\ +A^{2}BCab^{5} − A^{2}CMab^{4} + A^{2}CMab^{5} \end{matrix} ) \\ \frac{db}{dα} & = & b (Cb + 1) ( \begin{matrix}a + Qb^{2} − ab + Mab^{2} + A^{2}ab^{2} \\ −A^{2}ab^{3} − Mab − A^{2}Mab^{3} + A^{2}Mab^{4} \end{matrix} ) . \end{matrix}

La matriz Jacobiana evaluada en el punto (0, 0) del campo vectorial V̄_{η} es:

DV̄_{η} (0, 0) = ( \begin{matrix}0 & 0 \\ 0 & 0 \end{matrix} ) .

De nuevo el punto (0, 0) es no hiperbóico.

Nuevamente utilizaremos el blowing-up direccional [14], en este caso en la dirección a (horizontal). Definamos

ϕ_{2}^{a}(p, s) = (p, ps) = (a, b)

y sustituyendo en el campo vectorial V̄_{η} obtenemos:

W_{η} : \{ \begin{matrix}\frac{dp}{dα} = p^{2} ( \begin{matrix}ps − A^{2}p^{2}s^{2} + A^{2}p^{3}s^{3} + Bps^{2} − Qps^{2} − Bp^{2}s^{2} \\ +Cp^{2}s^{2} − Mp^{2}s^{2} − Cps + Mps + A^{2}Bp^{3}s^{4} \\ −A^{2}Bp^{4}s^{4} − A^{2}Cp^{3}s^{3} + A^{2}Cp^{4}s^{4} + A^{2}Mp^{3}s^{3} \\ −A^{2}Mp^{4}s^{4} − BCp^{3}s^{3} + CMp^{2}s^{2} − CMp^{3}s^{3} \\ −CQp^{2}s^{3} − A^{2}BCp^{5}s^{5} + A^{2}CMp^{4}s^{4} − A^{2}CMp^{5}s \end{matrix} ) \\ \frac{ds}{dα} = ps ( \begin{matrix}−2ps + 2A^{2}p^{2}s^{2} − 2A^{3}p^{3}s^{3} − Bps^{2} + 2Qps^{2} \\ +Bp^{2}s^{2} − 2Cp^{2}s^{2} + 2Mp^{2}s^{2} + 2Cps − 2Mps \\ −A^{2}Bp^{3}s^{4} + A^{2}Bp^{4}s^{4} + 2A^{2}Cp^{3}s^{3} − 2A^{2}Cp^{4}s^{4} \\ −2A^{2}Mp^{3}s^{3} + 2A^{2}Mp^{4}s^{4} + BCp^{3}s^{3} − 2CMp^{2}s^{2} \\ +2CMp^{3}s^{3} + 2CQp^{2}s^{3} + A^{2}BCp^{5}s^{5} − 2A^{2}CMp^{4}s^{4} \\ +2A^{2}CMp^{5}s^{5} + 2 \end{matrix} ) \end{matrix}

Reescalando el tiempo considerando β = pα nos resulta el sistema:

W̄_{η} (p, r) : \{ \begin{matrix}\frac{dp}{dβ} = p ( \begin{matrix}ps − A^{2}p^{2}s^{2} + A^{2}p^{3}s^{3} + Bps^{2} − Qps^{2} − Bp^{2}s^{2} \\ +Cp^{2}s^{2} − Mp^{2}s^{2} − Cps + Mps + A^{2}Bp^{3}s^{4} \\ −A^{2}Bp^{4}s^{4} − A^{2}Cp^{3}s^{3} + A^{2}Cp^{4}s^{4} + A^{2}Mp^{3}s^{3} \\ −A^{2}Mp^{4}s^{4} − BCp^{3}s^{3} + CMp^{2}s^{2} − CMp^{3}s^{3} \\ −CQp^{2}s^{3}−A^{2}BCp^{5}s^{5}+A^{2}CMp^{4}s^{4}−A^{2}CMp^{5}s \end{matrix} ) \\ \frac{ds}{dβ} = s ( \begin{matrix}−2ps + 2A^{2}p^{2}s^{2} − 2A^{3}p^{3}s^{3} − Bps^{2} + 2Qps^{2} \\ +Bp^{2}s^{2} − 2Cp^{2}s^{2} + 2Mp^{2}s^{2} + 2Cps − 2Mps \\ −A^{2}Bp^{3}s^{4} + A^{2}Bp^{4}s^{4} + 2A^{2}Cp^{3}s^{3} − 2A^{2}Cp^{4}s^{4} \\ −2A^{2}Mp^{3}s^{3} + 2A^{2}Mp^{4}s^{4} + BCp^{3}s − 2CMp^{2}s^{2} \\ +2CMp^{3}s^{3} + 2CQp^{2}s^{3} + A^{2}BCp^{5}s^{5} \\ −2A^{2}CMp^{4}s^{4} + 2A^{2}CMp^{5}s^{5} + 2 \end{matrix} ) \end{matrix}

La matriz jacobiana y evaluada en el punto (0, 0) del campo vectorial W̄_{η} es:

DW_{η} (0, 0) = ( \begin{matrix}−1 & 0 \\ 0 & 2 \end{matrix} ) .

Por lo tanto el punto (0, 0) es un punto silla hiperbólico del campo vectorial W̄_{η} (p, r) con variedad estable horizontal y variedad inestable vertical. Por consiguiente, el punto (0, ∞) del campo vectorial Y_{η} (u, v) es un punto silla no hiperbólico con variedad inestable vertical dirigida hacia abajo.

Luego, las soluciones del sistema son acotadas.

Nota 2 Este resultado implica que el modelo está bien planteado (well-posed) en el sentido que si hay un tamaño reducido de la población de presas i.e., x (t) → 0, la población de depredadores no crece desproporcionadamente, esto es, no es posible que y (t) → ∞, cuando t → ∞.

3.1 Naturaleza de los puntos sobre los ejes

La naturaleza de los puntos de equilibrio queda determinada principalmente por el parámetro M, como queda en evidencia en los siguientes lemas:

Lema 3 El punto de equilibrio (0, 0) es

I) silla hiperbólica, si y sólo si, M > 0.

II) repulsor hiperbólico, si y sólo si, M < 0.

III) silla-nodo no-hiperbólico, si y sólo si, M = 0.

Demostración. La matriz Jacobiana es DY_{η} (0, 0) = ( \begin{matrix}−A^{2}M & 0 \\ 0 & A^{2}BC \end{matrix} ).

Claramente,

I) det DY_{η} (0, 0) < 0, si y sólo si, M > 0. Luego, el punto de equilibrio (0, 0) es silla.

II) det DY_{η} (0, 0) > 0, si y sólo si, M < 0.

Luego, la naturaleza del punto de equilibrio (0, 0) depende de la traza,

Se tiene que trDY_{η} (0, 0) = A^{2} (−M + BC) > 0, y por lo tanto, (0, 0) es repulsor hiperbólico.

III) Cuando M = 0, el punto de equilibrio (M, 0) coincide con (0, 0).

Luego, (0, 0) es silla-nodo no-hiperbólico repulsor.

Lema 4 El equilibrio (1, 0) es punto silla hiperbólico para todo η = (A, B, C, Q, M) ∈ ]0, 1[ × (ℝ_{0}^{+})^{3} ×]−1, 1[.

Demostración. La evaluación de la matriz Jacobiana es:

DY_{η} (1, 0) = ( \begin{matrix}− (A^{2} + 1) (C + 1) (1 − M) & −Q (C + 1) \\ 0 & B (A^{2} + 1) (C + 1) \end{matrix} ) .

Por lo tanto,

detDY_{η} (1, 0) = −B (A^{2} + 1)^{2} (C + 1)^{2} (1 − M) .

Luego, el equilibrio (1, 0) es silla hiperbólica para todo valor de M. Lema 5 El punto de equilibrio (M, 0) es

I) repulsor hiperbólico, si y sólo si, M ∈ (]−1, 1[ − {0}).

II) repulsor no-hiperbólico, si y sólo si, M = 0. Demostración. La evaluación de la matriz Jacobiana es

DY_{η}(M, 0 = ( \begin{matrix}M(A^{2} + M)(C + M)(1 − M) & −QM^{2}(C + M) \\ 0 & B(A^{2} + M)(C + M) \end{matrix} ) .

Luego,

det DY_{η} (M, 0) = BM (1 − M) (M + A^{2})^{2} (C + M)^{2} > 0,

para todo −1 < M << 1. Por lo tanto, su naturaleza depende de la traza,

trDY_{η}(M, 0) = M(A^{2} + M)(C + M)(1 − M) + B(A^{2} + M)(C + M). Entonces,

trDY_{η} (M, 0) > 0, si y sólo si,M ∈ (]−1.1[ − {0}) .

Esto es, (M, 0) es repulsor hiperbólico.

Si M = 0, se produce el colapso con el equilibrio (0, 0), y por lo tanto, es silla-nodo repulsor no-hiperbólico. Lema 6 El punto (0, C) es:

I) Atractor hiperbólico, si y sólo si, M > 0.

II) Silla hiperbólica, si y sólo si, M < 0. III) Silla-nodo atractor no-hiperbólico, si y sólo si, M = 0.

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 29(1): 105–138, Jan–Jun 2022

Demostración. La evaluación de la matriz Jacobiana es

DY_{η} (0, C) = ( \begin{matrix}−A^{2}M & 0 \\ A^{2}BC & −A^{2}BC \end{matrix} ) .

Por lo tanto, det DY_{η} (0, C) = A^{4}BCM.

I) Si M > 0, implica que detDY_{η} (0, C) > 0 y su naturaleza depende de

la traza,

trDY_{η} (0, C) = −A^{2} (BC + M) < 0.

Luego, el punto de equilibrio (0, C) es un nodo atractor.

II) Si M < 0, se tiene que detDY_{η} (0, C) < 0 y el punto de equilibrio (0, C) es silla hiperbólica.

III) Como los cambios en la naturaleza de los puntos de equilbrio son suaves, cuando M = 0, se tiene que el equilibrio pasa de ser nodo atractor a silla. Por tanto, en este caso es un silla-nodo atractor.

Dada la dificultad para obtener las coordenadas de los puntos de equilibrio positivos, analizaremos uno de los casos particulares.

En lo que sigue asumiremos que M = 0 y por lo tanto, Q = \frac{A^{2}}{C}, lo cual disminuye la cantidad de parámetros.

3.2 Efecto Allee débil

Asumiendo que la población de presas está afectada por un efecto Allee débil especial se puede determinar explícitamente las coordenadas de los puntos de equilibrio positivos. El sistema (2) se puede reescribir con sólo tres parámetros de la siguiente forma:

Y_{π} (u, v) : \{ \begin{matrix}\frac{du}{dτ} & = & u^{2}(u + C) ((1 − u) (u^{2} + A^{2}) − \frac{A^{2}}{C} v) \\ \frac{dv}{dτ} & = & Bv(u + C − v)(u^{2} + A^{2}) \end{matrix}  (7)

con π = (A, B, C) ∈ ]0, 1[ × (ℝ_{0}^{+})^{2}.

Los puntos de equilibrio de este sistema (7) son: (0, 0), (0.C), (1, 0) y (u_{e}, u_{e} + C) donde u_{e} satisface la ecuación polinomial p_{010} (u) = 0.

I) Cuando Δ_{0} > 0, las raíces se pueden reescribir como:

u_{10} = \frac{1 − \sqrt{\frac{1}{C} (C − 4A^{2} (C + 1))}}{2} y u_{20} = \frac{1 + \sqrt{\frac{1}{C} (C − 4A^{2} (C + 1))}}{2}

Claramente, 0 < u_{10} < u_{20} < 1.

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 29(1): 105–138, Jan–Jun 2022

II) Existen una única raíz real positiva, si sólo si, Δ_{0} = 0, o sea, cuando A^{2} = \frac{C}{4(C+1)}, la cual es dada por u_{∗0} = \frac{1}{2}.

III) Existen raíces complejas conjugadas, si sólo si, Δ_{0} < 0, o A^{2} > \frac{C}{4(C+1)}.

Reiteramos que el parámetro B, no influye en la determinación de los puntos de equilibrio. La matriz Jacobiana del sistema (7) es:

DY_{π}(u, v) = ( \begin{matrix}DY_{π}(u, v)_{11} & −\frac{A^{2}}{C} u^{2}(C + u) \\ Bv(2Cu − 2uv + A^{2} + 3u^{2}) & B(A^{2} + u^{2})(C + u − 2v) \end{matrix} )

siendo DY_{π} (u, v)_{11} = \frac{∂}{∂u} (u^{2}(u + C) ((1 − u) (u^{2} + A^{2}) − Qv)).

3.3 Resultados principales para el sistema (7)

La naturaleza de los puntos de equilibrio sobre los ejes del sistema (7) o campo vectorial Y_{π} (u, v) ha sido establecida en los lemas 2, 3, 4 y 5.

Por lo tanto, se tiene que

• El equilibrio (0, 0) es un repulsor no-hiperbólico.

• El punto (1, 0) es silla hiperbólica para todo valor de parámetros.

• El punto (0, C) es atractor no-hiperbólico.

3.3.1 Naturaleza de los puntos de equilibrio positivos

De acuerdo al análisis anterior se deben estudiar tres casos, según existan dos, uno o ningún punto de equilibrio al interior del primer cuadrante. Sólo nos enfocaremos a la alternativa en que existen dos puntos de equilibrio positivos. Si ellos coinciden se tiene un punto de equilibrio positivo no-hiperbólico.

Caso 1 Existen dos puntos de equilibrio positivos.

La matriz Jacobiana en los puntos de equilibrio positivos (u, u + C) es:

DY_{π} (u, u + C) = ( \begin{matrix}DY_{π} (u, v)_{11} & −\frac{A^{2}}{C} u^{2} (C + u) \\ B (A^{2} + u^{2}) (C + u) & −B (A^{2} + u^{2}) (C + u) \end{matrix} ) .

Sea H_{0} = ((1 − u) (u^{2} + A^{2}) − \frac{A^{2}}{C} v). Luego,

DY_{π} (u, v)_{11} = u^{2}H_{0} + 2u(u + C)H_{0} + u^{2}(u + C)\frac{∂}{∂u}H_{0}

= u^{2}H_{0} + 2u(u + C)H_{0} + u^{2}(u + C) (−A^{2} − 3u^{2} + 2u) .

Por lo tanto, el determinante es

det DY_{π} (u, u + C) = Bu^{2} (A^{2} + u^{2}) (C + u)^{2} \frac{A^{2}C + 3Cu^{2} − 2Cu + A^{2}}{C}

el cual depende del factor T_{1} = A^{2}C + 3Cu^{2} − 2Cu + A^{2}.

La traza es dada por

trDY_{π} (u, u + C) = (u + C) (u^{2} (−A^{2} − 3u^{2} + 2u) − B (A^{2} + u^{2})) ,

la cual a su vez depende del factor

T_{2} = u^{2} (−A^{2} − 3u^{2} + 2u) − B (A^{2} + u^{2}) .

Lema 7 El equilibrio P_{1} = (u_{1}, u_{1} + C) es un punto silla hiperbólico.

Demostración. Sabemos que el signo del detDY_{π} (u, u + C) depende del factor:

T_{1} = 3Cu^{2} − 2Cu + (A^{2}C + A^{2}) .

Considerando la ecuación (6), se tiene que:

A^{2} (1 + \frac{1}{C} ) = u − u^{2}.

Reemplazando se obtiene que

T_{1} = 3u^{2} − 2u + (u − u^{2}) = u (2u − 1).

Evaluando para u = u_{1}, se llega

T_{1} (u_{1}) = −u_{1}\sqrt{C − 4A^{2} (C + 1)} < 0.

Por lo tanto, detDY_{π} (u_{1}, u_{1} + C) < 0. Luego, de acuerdo al Teorema de la traza y el determinante, el punto P_{1} = (u_{1}, u_{1} + C) es un punto silla hiperbólico.

Antes de establecer la naturaleza del punto P_{2} = (u_{2}, u_{2} + C) debemos mostrar algunas propiedades del sistema, relacionadas con las variedades estable e inestable del equilibrio P_{1} y la variedad inestable del punto silla (1, 0).

Se denotan W_{+}^{s} (P_{1}) la variedad estable por arriba, W_{+}^{u} (P_{1}) la variedad inestable por derecha del punto (u_{1}, u_{1} + C), y W_{+}^{u} (1, 0) la variedad inestable del punto silla (1, 0).

Teorema 1 Existe un subconjunto en el espacio de parámetros, para el cual las variedades W_{+}^{s} (P_{1}) y W_{+}^{u} (P_{1}) determinan una curva homoclínica que rodea el punto P_{2} = (u_{2}, u_{2} + C).

Demostración. Usando el Teorema de Existencia y Unicidad de las soluciones [10, 14, 34] y la geometría de las variedades estable e inestable del punto silla P_{1} = (u_{1}, u_{1} + C).

La variedad inestable derecha W_{+}^{u} (P_{1}) del punto silla P_{1} no puede cruzar la recta hacia la derecha de la recta u = 1, que es la frontera de la región de invarianza Γ̄.

Por lo tanto, esta solución deben cambiar su dirección y retornar hacia la izquierda.

A su vez, el α − limite de la variedad estable superior W_{+}^{s} (P_{1}) puede estar:

i) en el punto (∞, 0) del campo vectorial compacto del sistema (7),

ii) el punto P_{2}, cuando este es repulsor, o

iii) un ciclo límite inestable rodeando al punto de equilibrio P_{2} cuando este es atractor.

Luego, existen condiciones en el espacio de parámetros para los cuales W_{+}^{s} (P_{1}) ∩ W_{+}^{u} (P_{1}) ≠ ∅.

Entonces, existe un curva homoclínica, creada por el punto silla P_{1} y rodeando al punto P_{2}.

Teorema 2 Sean W_{+}^{u} (1, 0) la variedad inestable superior del punto silla (1, 0) y W_{+}^{s} (P_{1}) la variedad estable superior del punto (u_{1}, u_{1} + C). Existe un subconjunto de parámetros para los cuales W_{+}^{u} (1, 0) ∩ W_{+}^{s} (P_{1}) ≠ ϕ, formando una curva heteroclínica Σ̄.

Demostración. Usando también el Teorema de Existencia y Unicidad de las soluciones [10, 14, 34] y la geometría de las variedades inestables del punto silla (1, 0) y la variedad estable superior del punto P_{1} = (u_{1}, u_{1} + C).

Como se tiene en el Teorema 1, el α − limite de la variedad estable superior W_{+}^{s} (P_{1}) puede estar en:

i) el punto (∞, 0) del campo vectorial compacto del sistema (7),

ii) el punto P_{2}, cuando este es repulsor, o

iii) un ciclo límite inestable rodeando al punto de equilibrio P_{2} cuando este es atractor.

A su vez, el ω − limite de la variedad inestable superior del punto silla (1, 0), puede ser

i) el punto P_{2}, cuando éste es atractor, o

ii) un ciclo límite estable rodeando al punto de equilibrio P_{2} cuando este es repulsor,

iii) el punto (0, C), o bien

iv) el punto P_{1}, coincidiendo ambas variedades.

Luego, de acuerdo a iv) se tiene que W_{+}^{u} (1, 0) ∩ W_{+}^{s} (P_{1}) ≠ ϕ.

Nota 3 Los puntos con condiciones iniciales que están sobre la curva Σ̄ = W_{+}^{u} (1, 0) ∩ W_{+}^{s} (P_{1}), tienen como ω − limite, el equilibrio (0, C). Esto implica que existen trayectorias con condicioes iniciales muy cercanas, pero a diferente lado de Σ̄, pueden tener conjuntos ω − limite muy alejados. Es decir, las soluciones son altamente sensibles de las condiciones iniciales.

Además, es sabido que el rompimiento de una curva homoclínica o heteroclínica originan ciclos límites no-infinitesimales [16], los que podrían coincidir con los ciclos límites originados por bifuración de Hopf. La estabilidad de cada uno de ellos debe ser determinada [4].

La naturaleza del punto de equilibrio (u_{2}, u_{2} + C) depende, no sólo de los valores propios de la matriz Jacobiana, sino esencialmente de las posiciones relativas de las variedades W_{+}^{u} (P_{1}) y W_{+}^{u} (1, 0), como se establece en el Teorema 3.

Teorema 3 Consideremos los puntos (u^{∗}, v^{u}) ∈ W_{+}^{u} (1, 0) y (u^{∗}, v^{s}) ∈ W_{+}^{u} (P_{1}), con 0 < u_{1} < u^{∗} < 1. La relación entre v^{u} y v^{s}, determina la estabilidad del punto (u_{2}, u_{2} + C) del modo siguiente:

1. Asumiendo que v^{s} > v^{u}, entonces el punto de equilibrio (u_{2}, u_{2} + C) es:

(a) Un atractor hiperbólico, si y solo si, B > \frac{u_{2}^{2}(−A^{2}−3u_{2}^{2}+2u_{2})}{(A^{2}+u_{2}^{2})}.

(b) Un repulsor hiperbólico, si y sólo si, B < \frac{u_{2}^{2}(−A^{2}−3u_{2}^{2}+2u_{2})}{(A^{2}+u_{2}^{2})}.

(c) Un foco débil, si y sólo si, B = \frac{u_{2}^{2}(−A^{2}−3u_{2}^{2}+2u_{2})}{(A^{2}+u_{2}^{2})}.

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 29(1): 105–138, Jan–Jun 2022

2. Suponiendo que v^{s} < v^{u}, entonces el punto de equilibrio (u_{2}, u_{2} + C) es

(a) Un atractor rodeado por un ciclo límite inestable.

(b) Un nodo o foco repulsor hiperbólico y las trayectorias del sistema (7) tienen al punto (0, C) como su ω−limite, el cual es un equilibrio casi globalmente estable [33, 36].

Demostración.

1. Supongamos que v^{s} > v^{u}.

Entonces W_{+}^{s} (P_{1}) la variedad estable superior del punto (u_{1}, u_{1} + C), está por encima de W_{+}^{u} (1, 0) la variedad estable superior del punto (1, 0).

Luego el ω − limite de W_{+}^{u} (1, 0) debe ser:

i) el punto P_{1}, originando una curva heteroclínica, ii) el punto P_{2}, cuando este es atractor (foco o nodo), o iii) un ciclo límite atractor.

Evaluando el factor T_{1} en el punto P_{2} se tiene que

T_{1} = 3Cu_{2}^{2} − 2Cu_{2} + (A^{2}C + A^{2}) .

Siguiendo el mismo procedimiento usado en el caso del punto P_{1}, se tiene que:

T_{1} = u_{2}\sqrt{C − 4A^{2} (C + 1)} > 0.

Luego, detDY_{η} (u_{2}, u_{2} + C) > 0.

Por lo tanto, la naturaleza depende del signo de la traza.

Luego, cuando T_{2} = u_{2}^{2} (−A^{2} − 3u_{2}^{2} + 2u_{2}) − B (A^{2} + u_{2}^{2}) es positivo, negativo o cero, se tiene la tesis.

Si la traza es positiva, de acuerdo con el Teorema de Existencia y Unicidad y por el Teorema de Poincaré-Bendixson, debe existir un ciclo límite rodeando al punto de equilibrio P_{2}.

2. Supongamos que v^{s} < v^{u}.

Entonces W_{+}^{s} (P_{1}) la variedad estable superior del punto P_{1}, está por debajo de W_{+}^{u} (1, 0) la variedad estable superior del punto (1, 0).

(a) Si T_{2} < 0 el punto de equilibrio P_{2} es un atractor hiperbólico. Entonces el ω − limite de W_{+}^{s} (u_{1}, u_{1} + C) es: a1) el punto P_{2}, o bien a2) un ciclo límite atractor rodeando a un ciclo límite inestable y al punto de equilibrio (u_{2}, u_{2} + C), existiendo dos ciclos concéntricos.

(b) Si T_{2} > 0, el punto de equilibrio P_{2} es repulsor (foco o nodo). Entonces, las soluciones con α − limite en las vecindades del punto P_{2} tienen como ω − limite: b1) un ciclo límite atractor (o estable), o bien, b2) el punto (0, C).

Esto completa la demostración.

Corolario 1 Existe una bifurcación de Hopf [10, 14, 27] en el punto de equilibrio P_{2} para el valor de la bifurcación S = (M − 2u_{e2} + 1) u_{e2}.

Demostración. La demostración sigue del Teorema 3, cuando:

det DY_{η} (u_{2}, u_{2} + C)
= Bu_{2}^{2} (A^{2} + u_{2}^{2}) (C + u_{2})^{2} \frac{3Cu_{2}^{2} − 2Cu_{2} + (A^{2}C + A^{2})}{C},

es positivo, y

trDY_{η} (u_{2}, u_{2} + C) = (u_{2} + C) (u_{2}^{2} (−A^{2} − 3u_{2}^{2} + 2u_{2}) − B (A^{2} + u_{2}^{2}))

cambia de signo.

Verificando la condición de transversalidad [10, 16, 34] se tiene que

\frac{∂}{∂B} (trDY_{η} (u_{2}, u_{2} + C)) = − (u_{2} + C) (A^{2} + u_{2}^{2}) < 0.

Por lo tanto, existe al menos un ciclo límite.

Nota 4 Los resultados obtenidos para la naturaleza de los puntos de equilibrio P_{1} = (u_{1}, u_{1} + C), P_{2} = (u_{2}, u_{2} + C) y (0, C) establecen que, dependiendo de las respectivas condiciones en el espacio de parámetros, los tamaños poblacionales de presas y depredadores pueden:

i) oscilar alrededor de tamaños fijos, o

ii) tender a tamaños fijos, o

iii) tender a cero en el caso de las presas (hay extinción), mientras que los depredadores alcanzan un valor máximo poblacional.

Por otra parte, en la Figura 5, se muestra que puede existir un segundo ciclo límite concéntrico al punto P_{2}. En tal caso, existe el fenómeno de triestabilidad, por cuanto, para el mismo subconjunto de parámetros se tiene que P_{2} es local asintóticamente estable, coexistiendo con un ciclo límite atractor y con el punto atractor local no-hiperbólico (0, C). De acuerdo a los tamaños poblacionales iniciales, puede darse cualquiera de las situaciones descritas en i), ii) o iii).

Caso 2 Existe un único punto de equilibrio positivo.

Esto sucede para Δ = 1 − \frac{4A^{2}}{C} (C + 1) = 0, teniendo que A^{2} = \frac{C}{4(C+1)}.

Teorema 4 Bifurcación de Bogdanov-Takens [10, 27]

a) El punto P_{∗} = (u_{∗}, u_{∗} + C), colapso de P_{1} y P_{2} es un silla-nodo atractor.

b) El punto P_{∗} = (u_{∗}, u_{∗} + C), colapso de P_{1} y P_{2} es un punto cúspide, si y solo si, trDY_{η} (u_{∗}, u_{∗} + C) = 0 y detDY_{η} (u_{∗}, u_{∗} + C) = 0.

Demostración.

a) Como u_{∗} = \frac{1}{2} y

det DY_{η} (u_{∗}, u_{∗} + C) = Bu_{2}^{2} (A^{2} + u_{2}^{2}) (C + u_{2})^{2} T_{2}

con T_{2} = u_{∗}\sqrt{Δ}. Es claro que detDY_{η} (u_{∗}, u_{∗} + C) = 0, cuando colapsan los puntos P_{1} y P_{2}, ya que Δ = 0.

Luego, se tiene que P_{∗} es un silla-nodo.

Antes de esta coincidencia, los puntos P_{1} y P_{2} son silla y nodo atractor, respectivamente. Luego, P_{∗} es un silla-nodo atractor [10, 27].

b) Como trDY_{η} (u_{∗}, u_{∗} + C) depende del factor

T_{2} = (u_{∗})^{2} (−A^{2} − 3 (u_{∗})^{2} + 2u_{∗}) − B (A^{2} + (u_{∗})^{2}) ,
T_{2} = (1/2)^{2} (−A^{2} − 3 (1/2)^{2} + 2 (\frac{1}{2})) − B (A^{2} + (1/2)^{2})
= (−A^{2} − \frac{1}{4}) B + (\frac{1}{16} − \frac{1}{4}A^{2}) .

Luego, si B = \frac{(\frac{1}{16} − \frac{1}{4} A^{2})}{(A^{2}+ \frac{1}{4})}, entonces trDY_{η} (\frac{1}{2}, \frac{1}{2} + C) = 0. Por lo tanto, P_{∗} es un punto cúspide.

Nota 5 En el caso en que P_{∗} es un punto cúspide existe una única solución que llega a ese punto. mientras que todas las demás tienden al punto (0, C), el cual es casi global asintóticaente estable [33, 36].

Por ser una condición cerrada (una igualdad) esta situación es difícil de detactarla mediante simulaciones numéricas.

4 Simulaciones numéricas y comparación entre modelos

Cuando existen más de dos puntos de equilibro positivos, el estudio de la dinámica del sistema se complica. En el modelo propuesto no es simple determinar las coordenadas exactas de los cuatro puntos de equilibrio positivos que pueden existir para cierta condición de parámetros. Por eso mostramos simulaciones que esclarecen las dinámicas del modelo si M = 0.

4.1 Simulaciones

Figura 1: Existencia de dos puntos de equilibrios positivos. Para A = 0.3, B = 0.04, C = 0.6, Q = 0.15, M = 0, P_{1} es un punto silla y P_{2} es un foco inestable rodeado de un ciclo límite estable generado vía bifurcación de Hopf. El punto (0, C) es atractor local no hiperbólico.

influencia del efecto allee débil en las presas en un... 127

Figura 2: Existencia de dos puntos de equilibrios positivos. Para A = 0.3, B = 0.01, C = 0.6, Q = 0.15, M = 0, P_{1} es un punto silla y P_{2} es un foco estable. El punto (0, C) es atractor local no hiperbólico.

Figura 3: Existencia de un punto de equilibrio positivo. Para A = 0.289, B = 0, 5, C = 0, 5, Q = 1.666, M = 0, la única singularidad positiva es un punto silla y el punto de equilibrio (0, C) es un atractor no hiperbólico, casi global asintóticamente estable [33, 36].

128 s. valenzuela-figueroa et al.

Figura 4: No hay existencia de puntos de equilibrio positivo. Para A = 0.289, B = 0.125, C = 0, 5, Q = 0.251505, M = 0, el único punto de equilibrio (0, C) es un atractor no hiperbólico, global asintóticamente estable.

Figura 5: Existencia de un punto de equilibrio positivo. Para A = 0.26775, B = 0.109, C = 0.125, Q = 0.251525, M = 0, se tiene que A^{2} − CQ > 0. La única singularidad positiva es un punto local asintóticamente estable rodeado de dos ciclos límites, el interior inestable y el exterior estable. El punto de equilibrio (0, C) es un atractor no hiperbólico.

4.2 Comparación entre modelos

En la Tabla 1 mostramos varios modelos que consideran una respuesta funcional sigmoidea describiendo la acción de los depredadores, destacando el máximo de puntos de equilibrio positivo que pueden existir.

Recordemos que en el modelo propuesto por P. H. Leslie in 1948 [28], se asume la respuesta funcional lineal h (x) = qx; tanto el crecimiento de la presas como el crecimiento de los depredadores son representados por la función logística; entonces, el modelo de Leslie es descrito por el sistema bidimensional de ecuaciones diferenciales ordinarias:

X_{α} (x, y) : \{ \begin{matrix}\frac{dx}{dt} & = & r (1 − \frac{x}{K} ) x − qxy \\ \frac{dy}{dt} & = & s (1 − \frac{y}{nx}) y \end{matrix}

con α = (r, K, q, s, n) ∈ ℝ_{+}^{5}.

La capacidad de soporte del medio ambiente de los depredadores es K_{y} (x) = nx, esto es, se asume proporcional a la población de presas. El sistema tiene un único punto de equilibrio positivo que es global asíntóticamente estable (gas) lo cual se prueba construyendo funciones de Lyapunov adecuadas [24].

En la Tabla 1 que se presenta a continuación, se usa la siguiente notación:

• N: cantidad máxima de puntos de equilibrio positivos,

• C: indica la cantidad de ciclos límites,

• D: son los dominios de los sistemas, los cuales pueden ser los conjuntos:

– Ω_{1} = {(x, y) ∈ ℝ^{2} : x > 0; y ≥ 0},

– Ω_{2} = {(x, y) ∈ ℝ^{2} : x ≥ 0; y ≥ 0},

• μ, η, ν, μ, λ, π y σ indican los respectivos vectores de parámetros positivos,

• f_{1} (x) = r (1 − \frac{x}{K} ) x, es la función logística,

• f_{2} (x) = r (1 − \frac{x}{K} ) (x − m) x, función de crecimiento de las presas incluyendo efecto Allee simple,

• f_{3} (x) = r (1 − \frac{x}{K} ) (1 − \frac{m+b}{x+b}) x, función de crecimiento de las presas incluyendo efecto Allee doble.

Tabla 1: Comparación entre modelos de depredación del tipo Leslie-Gower modificados, con respuesta funcional sigmoidea presa-dependiente.

Respuesta funcional h (x)

Modelo Leslie-Gower modificado N C D Ref.

1) h_{1} (x) = \frac{q x^{2}}{x^{2} + a^{2}}  X_{π} : \{ \begin{matrix}\frac{dx}{dt} & = & f_{1} (x) − \frac{q x^{2}}{x^{2} + a^{2}} y \\ \frac{dy}{dt} & = & s (1 − \frac{y}{nx}) y \end{matrix}  1 2 Ω_{1} [22]
2) h_{1} (x) = \frac{q x^{2}}{x^{2} + a^{2}}  X_{σ} : \{ \begin{matrix}\frac{dx}{dt} & = & f_{1} (x) − \frac{q x^{2}}{x^{2} + a^{2}} y \\ \frac{dy}{dt} & = & s (1 − \frac{y}{nx+c}) y \end{matrix}  3 2 Ω_{2} [44]
3) h_{2} (x) = \frac{q x^{2}}{x^{2} +bx+ a^{2}} X_{μ} : \{ \begin{matrix}\frac{dx}{dt} & = & f_{1} (x) − \frac{q x^{2}}{x^{2} +bx+ a^{2}} y \\ \frac{dy}{dt} & = & s (1 − \frac{y}{nx+c}) y \end{matrix} 3 2 Ω_{2} [38]
4) h_{2} (x) = \frac{q x^{2}}{x^{2} +bx+ a^{2}} X_{λ} : \{ \begin{matrix}\frac{dx}{dt} & = & f_{1} (x) − \frac{q x^{2}}{x^{2} +bx+ a^{2}} y \\ \frac{dy}{dt} & = & s (1 − \frac{y}{nx}) y \end{matrix} 3 4 Ω_{1} [13]
5) h_{1} (x) = \frac{q x^{2}}{x^{2} + a^{2}}  X_{μ} : \{ \begin{matrix}\frac{dx}{dt} & = & f_{2} (x) − \frac{q x^{2}}{x^{2} + a^{2}} y \\ \frac{dy}{dt} & = & s (1 − \frac{y}{nx+c}) y \end{matrix}  4 2 Ω_{2} [46]
6) h_{1} (x) = \frac{q x^{2}}{x^{2} + a^{2}}  X_{η} : \{ \begin{matrix}\frac{dx}{dt} & = & f_{3} (x) − \frac{q x^{2}}{x^{2} + a^{2}} y \\ \frac{dy}{dt} & = & s (1 − \frac{y}{nx}) y \end{matrix}  4 2 Ω_{1} [5]
7) h_{2} (x) = \frac{q x^{2}}{x^{2} +bx+ a^{2}} X_{ν} : \{ \begin{matrix}\frac{dx}{dt} & = & f_{2} (x) − \frac{q x^{2}}{x^{2} +bx+ a^{2}} y \\ \frac{dy}{dt} & = & s (1 − \frac{y}{nx}) y \end{matrix} 4 2 Ω_{1} [37]

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 29(1): 105–138, Jan–Jun 2022

Observamos que en los modelos propuestos la cantidad de puntos de equilibrio positivos aumenta de 1 a 3, cuando se considera un alimento alternativo para los depredadores o cuando se cambia la respuesta funcional h_{1} (x) a h_{2} (x) (diferencia entre el modelo 1 con los modelos 2, 3 y 4). Ese número aumenta a cuatro (4) si se considera efecto Allee (simple o doble) considerando las respuestas funcionales h_{1} (x) o h_{2} (x).

También, la cantidad de ciclos limites puede cambiar aumentando su número de 0 a 1 o 2. En el artículo de Dai et al. [13], se reporta la existencia de hasta cuatro ciclos límites, con expresiones matemáticas complicadas.

La diferencias dinámicas entre los modelos son en general bastante significativas. Por ejemplo, puede implicar la desaparición de la población de presas cuando el punto (0, C) es atractor (hiperbólico o no-hiperbólico). Esta puede ser una situación deseable en alguna interacción de depredación si la presa constituye una plaga.

Hacemos presente que existen otras formas matemáticas para representar una respuesta funcional sigmoidea como son:

i) h (x) = \frac{q x^{m}}{x^{m} + a^{m}}, con m ∈ ℕ, m ≥ 2, a y q > 0.

ii) h (x) = q tanh x, con q > 0.

En general, un modelo de depredador-presa creíble tambián debería poseer propiedades biológicas y ecológicas razonables [7]. Estimamos que los modelos de la Tabla 1, satisfacen los atributos y condiciones establecidas en [7] que lo hacen verosímiles.

Sin embargo, debe tenerse en cuenta que la complejidad dinámica de algunos de los modelos (complejidad matemática), no implica una complejidad ecológica; esto significa que existen interacciones de depredación en la naturaleza donde los depredadores son generalistas, en las cuales no necesariamente se satisfacen todas las propiedades analíticas del modelo estudiado.

5 Conclusiones

En este trabajo hemos establecido las principales propiedades de un modelo de depredación del tipo Leslie-Gower, considerando que los depredadores tienen un alimento alternativo, cuando su presa favorita escasea; hemos asumido también que la población de presas es afectada por el fenómeno denominado efecto Allee [6, 8, 30] y que la acción de los depredadores es descrita por una respuesta funcional sigmoidea o Holling tipo III.

Estudiamos parcialmente un sistema topológicamente equivalente al original, dependiente sólo de cuatro parámetros que se obtuvo mediante un difeomorfismo. En este segundo sistema (y por tanto en el original) existen diferentes situaciones dinámicas según las relaciones obtenidas entre los parámetros.

Dada la dificultad para obtener explícitamente los puntos de equilibrio analizamos un caso particular del efecto Allee, suponiendo M = 0. Establecimos condiciones en el espacio de parámetros para la existencia de dos puntos de equilibrio positivos, uno de los cuales el punto P_{1} = (u_{1}, u_{1} + C), es siempre un punto silla.

Entre los resultados más importantes, se ha demostrado la existencia de una curva separatriz Σ̄, dividiendo el comportamiento de las soluciones o trayectorias del sistema en el plano de fase. Esta curva separatriz es generada por la variedad estable superior W_{+}^{s} (P_{1}) del punto silla (u_{1}, u_{1} + C).

Esto implica que para el mismo conjunto de parámetros, dos soluciones muy cercanas pero en un lado diferente de esa separatriz Σ̄, tendrían ω − limites diferentes y lejanos, lo cual significa que las soluciones son altamente sensibles a las condiciones iniciales.

Todo lo anterior indica que resulta difícil predecir si las poblaciones coexisten o una de ellas se extingue, para un mismo subconjunto de parámetros. Sin embargo, la naturaleza del punto (0, C) implica que la extinción de la población de presas es posible y los depredadores crecen hasta su máxima capacidad de carga C.

Se tiene entonces que para un mismo tamaño poblacional de las presas, pero diferentes tamaños poblacionales de los depredadores, ambas poblaciones podrían coexistir, o bien, la población de presas podría desaparecer, mientras que la población de depredadores se mantiene en el largo plazo. Esta última situacion sería deseable, si las presas constituyen una plaga.

Nuestros resultados también confirman el hecho de que si la población de depredadores es generalista, la extinción de cada especie puede evitarse, para ciertas condiciones iniciales [21, 25, 32, 35]. Consideramos que este modelo, representa bien algunas interacciones del mundo real, dada la variedad de comportamientos dinámicos que describe. Además, estimamos que el sistema estudiado cumple con los atributos y condiciones que se describen en el articulo de Berryman et al. [7] .

El modelo estudiado tiene una rica e interesante dinámica. Esta riqueza queda en evidencia cuando se compara con la dinámica del modelo del tipo Leslie-Gower modificado, en el cual la tasa de crecimiento de las presas es del tipo logístico y los depredadores tienen también un alimento alternativo, pero en el que el efecto Allee en las presas no es asumido. Sin embargo, el modelo estudiado en este trabajo, tiene una dinámica más cercana a otros modelos del tipo Leslie-Gower como son los modelos estudiados en [21] y en [22].

Una cuestión que hemos dejado pendiente se refiere a establecer la unicidad o multiplicidad de los ciclos límites del sistema. La determinación de la cantidad de ciclos límites en modelos de depredación [11] está relacionado con el Problema 16 de Hilbert [16]. Este problema formulado el año 1900, por el matemático alemán David Hilbert, se refiere a la cantidad y posiciones relativas de ciclos límites en sistemas polinomiales de ecuaciones diferenciales ordinarias [16, 17], y aún permanece abierto.

Numéricamente, se muestra la existencia de dos ciclos límites concéntricos al punto P_{2} (en la Figura 5, el interior inestable y el exterior estable. Por tanto, existe el fenómeno de triestabilidad, por cuanto, P_{2} es un equilibrio local asintóticamente estable, coexistiendo con un ciclo límite estable (atractor) y con el punto atractor local no hiperbólico (0, C).

Extensiones de este trabajo para desarrollar en futuros artículos basados en el modelo de Leslie-Gower son: I) Asumir los casos en que el modelo tiene tres o cuatro puntos de equilibrio positivo, esto es, considerar efecto Allee fuerte. II) Considerar el efecto Allee aditivo en lugar del efecto Allee simple asumido en este trabajo. III) Determinar la influencia del exponente n ∈ ℕ. de la respuesta funcional generalizada, considerando ya sea efecto Allee simple o múltiple [6].

Agradecimientos

Los autores agradecen a los árbitros anónimos por expresar comentarios y sugerencias que permitieron mejorar este manuscrito.

Financiamiento

No hay fondos comprometidos para el desarrollo de este artículo, salvo los salarios percibidos por los autores en sus respectivas instituciones.

Referencias

[1] P. Aguirre, E. González-Olivares, E. Sáez, Three limit cycles in a Leslie-Gower predator-prey model with additive Allee effect, SIAM Journal on Applied Mathematics 69(2009) 1244–1262. Doi: 10.1137/070705210

[2] C. Arancibia-Ibarra, E. González-Olivares, A modified Leslie-Gower predator-prey model with hyperbolic functional response and Allee effect on prey, in: R. Mondaini (Ed.) BIOMAT 2010 International Symposium on Mathematical and Computational Biology, World Scientific, Singapore (2011) pp.146–162. Doi: 10.1142/9789814343435_0010

[3] N. Bacaër, A Short History of Mathematical Population Dynamics, Springer, New York NY, 2011. Doi: 10.1007/978-0-85729-115-8

[4] A.D. Bazykin, Nonlinear Dynamics of interacting populations, World Scientific, Singapore, 1998. Doi: 10.1142/2284

[5] R. Becerra-Klix, E. González-Olivares, A Leslie-Gower type predation model considering double Allee effect on prey and a sigmoid functional response, in J. Vigo-Aguiar (Ed.), Proceedings of the 17th International Conference on Computational and Mathematical Methods in Science and Engineering, CMMSE, Cádiz, España,

2017 pp.252–263. In: http:// rd.unir. net/ sisi/ research/

resultados/  1526642421Proceedings_CMMSE_2017_vol_
1_6.desprotegido_%20BEEBOTS.  PDF

[6] L. Berec, E. Angulo, F. Courchamp, Multiple Allee effects and population management, Trends in Ecology & Evolution, 22(2007) 185–191. Doi: 10.1016/j.tree.2006.12.002

[7] A.A. Berryman, A.P. Gutierrez, R Arditi, Credible parsimonious and useful predator-prey models- A reply to Abrams, Gleesson, and Sarnelle, Ecology 76(1995) 1980–1985. Doi: 10.2307/1940728

[8] D.S. Boukal, L. Berec, Single-species models and the Allee effect: Extinction boundaries, sex ratios and mate encounters, Journal of Theoretical Biology 218(2002), 375–394. Doi: 10.1006/jtbi.2002.3084

[9] C.W. Clark, Mathematical Bioeconomic: The Optimal Management of Renewable Resources, 2nd Ed., John Wiley and Sons, New York NY, 1990. Doi: 10.1137/1020117

[10] C. Chicone, Ordinary Differential Equations with Applications no.2, Texts in Applied Mathematics 34, Springer, 2006. Doi: 10.1007/0-387-35794-7

[11] K.-S. Cheng, Uniqueness of a limit cycle for a predator-prey system, SIAM Journal of Mathematical Analysis 12(1981), no. 4, 541–548. Doi: 10.1137/0512047

[12] F. Courchamp, L. Berec, J. Gascoigne, Allee effects in ecology and conservation, Oxford University Press, Oxford, 2008. Doi: 10.1093/acprof:oso/9780198570301.001.0001 [13] Y. Dai, Y. Zhao, B. Sang, Four limit cycles in a predator-prey

system of Leslie type with generalized Holling type III functional response, Nonlinear Analysis: Real World Applications 50(2019) 218–239. Doi: 10.1016/J.NONRWA.2019.04.003

[14] F. Dumortier, J. Llibre, J.C. Artés, Qualitative theory of planar differential systems, Springer, Berlin, 2006. Doi: 10.1007/978-3-540-32902-2

[15] H.I. Freedman, Deterministic Mathematical Model in Population Ecology, Marcel Dekker, New York NY, 1980. Doi: 10.2307/3556198

[16] V.A. Gaiko, Global Bifurcation Theory and Hilbert’s Sixteenth Problem, Mathematics and its Applications 559, Springer, Boston, MA, 2003. Doi: 10.1007/978-1-4419-9168-3

[17] V.A. Gaiko, V. Vuik, Global dynamics in the Leslie-Gower model with the Allee effect, International Journal of Bifurcation and Chaos, 28(2018), no.

12. Doi: 10.1142/S0218127418501511

[18] B-S. Goh, Management and Analysis of Biological Populations, Elsevier Scientific Publ. Co., 1980.

[19] E. González-Olivares, B. González-Yañez, J. Mena-Lorca, R. Ramos-Jiliberto, Modelling the Allee effect: Are the different mathematical forms proposed equivalents?, in: R. Mondaini (Ed.) Proceedings of the International Symposium on Mathematical and Computational Biology BIOMAT 2006, E-papers Serviços Editoriais Ltda., Rio de Janeiro, Brazil (2007), pp. 53–71.

[20] E. González-Olivares, A. Rojas-Palma, Multiple limit cycles in a Gause type predator-prey model with Holling type III functional response and Allee effect on prey, Bulletin of Mathematical Biology 35 (2011) 366–381. Doi: 10.1007/s11538-010-9577-5

[21] E. González-Olivares, J. Mena-Lorca, A. Rojas-Palma, J. D. Flores, Dynamical complexities in the Leslie-Gower predator-prey model as consequences of the Allee effect on prey, Applied Mathematical Modelling 35(2011) 366–381. Doi: 10.1016/j.apm.2010.07.001

[22] E. González-Olivares, P. Tintinago-Ruiz, A. Rojas-Palma, A Leslie-Gower type predator-prey model with sigmoid funcional response, International Journal of Computer Mathematics 93(2015), no. 9, 1895–1909. Doi: 10.1080/00207160.2014.889818 [23] E. González-Olivares, C. Arancibia-Ibarra, A. Rojas-Palma, B. González-

Yañez, Dynamics of a Leslie-Gower predation model considering a generalist predator and the hyperbolic functional response, Mathematical Biosciences and Engineering 16(2019), no. 6, 7995–8024. Doi: 10.3934/mbe.2019403

[24] E. González-Olivares E, A. Rojas-Palma, Global stability in a modified Leslie-Gower type predation model assuming mutual interference among generalist predators, Mathematical Biosciences and Engineering, 17(2020), no. 6, 7708–7731. Doi: 10.3934/mbe.2020392

[25] B. González-Yañez, E. González-Olivares, J. Mena-Lorca, Multistability on a Leslie-Gower type predator-prey model with nonmonotonic functional response, in: R. Mondaini, R. Dilao (Eds.), BIOMAT 2006 - International Symposium on Mathematical and Computational Biology (Manaus, Brazil), World Scientific, Singapour: 2007 pp.359–384. Doi: 10.1142/9789812708779_0023

[26] C.S. Holling, The components of predation as revealed by a study of smallmammal predation of the European pine sawfly, Canadian Entomologist 91(1959) 293–320. Doi: 10.4039/Ent91293-5

[27] Y.A. Kuznetsov, Elements of Applied Bifurcation Theory (3rd ed.), Springer, Cham, 2004. Doi: 10.1007/978-1-4757-3978-7

[28] P.H. Leslie, Some further notes on the use of matrices in population mathematics, Biometrika 35(1948) 213–245. Doi: 10.1093/biomet/35.3-4.213

[29] P.H. Leslie, J.C. Gower, The properties of a stochastic model for the predator-prey type of interaction between two species, Biometrika 47(1960) 219–234. Doi: 10.2307/2333294

[30] M. Liermann, R. Hilborn, Depensation: evidence, models and implications, Fish and Fisheries 2(2001), no. 1, 33–58. Doi: 10.1046/j.1467-2979.2001.00029.x

[31] N. Martínez-Jeraldo, E. Rozas-Torres, E. González-Olivares, Un modelo de depredación del tipo Leslie-Gower considerando depredadores generalistas y efecto Allee en las presas, Selecciones Matemáticas 8(2021), no. 1, 147–160. Doi: 10.17268/sel.mat.2021.01.14 [32] J. Mena-Lorca, E. González-Olivares, B. González-Yañez, The Leslie-

Gower predator-prey model with Allee effect on prey: A simple model with a rich and interesting dynamics, in: Mondaini R. (Ed.), Proceedings of the 2006 International Symposium on Mathematical and Computational Biology BIOMAT 2006, E-papers Serviços Editoriais Ltda., Rio de Janeiro, Brazil (2007) pp. 105–132. Doi: 10.1142/6483

[33] P. Monzón, Almost global attraction in planar systems, Systems & Control Letters 54(2005) 753–758. Doi: 10.1016/j.sysconle.2004.11.014

[34] L. Perko, Differential equations and dynamical systems(3rd ed.), Springer, New York NY, 2001. Doi: 10.1007/978-1-4613-0003-8

[35] L. Puchuri-Medina, E. González-Olivares, A. Rojas-Palma, Multistability in a Leslie-Gower type predation model with rational nonmonotonic functional response and generalist predators, Computational and Mathematical Methods 2(2020). Doi: 10.1002/cmm4.1070

[36] A. Rantzer, A dual to Lyapunov’s stability theorem, Systems and Control Letters 42(2001), no. 3, 161–168. Doi: 10.1016/S0167-6911(00)00087-6

[37] F.J. Reyes-Bahamón, Sobre la dinámica de algunos modelos depredadorpresa tipo Leslie con respuesta funcional no monótona y efecto Allee en las presas, Tesis de Maestría en Ciencias Matemática Aplicada, Universidad Nacional de Colombia, 2017. In: https:// repositorio. unal.

edu. co/handle/ unal/ 59746

[38] A. Rojas Palma, P.C. Tintinago-Ruiz, E. González-Olivares, Bifurcaciones en un modelo de depredación tipo Leslie-Gower modificado con una respuesta funcional Holling tipo III racional, en preparación, 2021.

[39] E. Sáez, E. González-Olivares, Dynamics on a Predator-prey Model, SIAM J. Applied Mathematics 59(1999) 1867–1878. Doi: https:// doi.

org/ 10.1137/ S0036139997318457 10.1137/S0036139997318457

[40] P.A. Stephens, W.J. Sutherland, R.P. Freckleton, What is the Allee effect?, Oikos 87(1999) 185–190. Doi: 10.2307/3547011

[41] J. Sugie, K. Miyamoto, K. Morino, Absence of limits cycle of a predatorprey system with a sigmoid functional response, Applied Mathematical Letter 9(1996), no. 4, 85–90. Doi: 10.1016/0893-9659(96)00056-0

[42] R.J. Taylor, Predation, Chapman and Hall, London, 1984. In: https://

www. worldcat. org/ title/ predation/ oclc/ 12453276

[43] P. Tintinago-Ruiz, L. Restrepo-Alape, E. González-Olivares, Consequences of weak Allee effect in a Leslie-Gower type predator-prey model with a generalized Holling type III functional response, in: G. Olivar Tost and O. Vasilieva (Eds.) Analysis, Modelling, Optimization, and Numerical Technique, Springer Proceedings in Mathematics & Statistics, 121(2015), no. 6, 89–103. Doi: 10.1007/978-3-319-12583-1_6

[44] P.C. Tintinago-Ruiz, E. González-Olivares, A. Rojas-Palma, Dinámicas de un modelo de depredación del tipo Leslie-Gower considerando respuesta funcional sigmoidea y alimento alternativo para los depredadores, en preparación, 2021.

[45] P. Turchin, Complex Population Dynamics: A Theoretical/Empirical Synthesis, Princeton University Press, Princeton NJ, 2003. Doi: 10.1515/9781400847280

[46] S. Valenzuela-Figueroa, Modelos de depredación del tipo Leslie-Gower con respuesta funcional sigmoidea y efecto Allee en las presas, Tesis Magister en Matemática, Pontificia Universidad Católica de Val-

paraíso, 2013. En: http:// repositorio.  conicyt. cl/handle/

10533/ 194403