Journal homepage    http://revistas.unitru.edu.pe/index.php/SSMM

Estabilidad en un modelo de depredación del tipo Leslie-Gower modificado considerando competencia entre los depredadores


Stability in a modified Leslie-Gower type predation model considering competence among predators

Eduardo González-Olivares

Gallegos-Zuñiga

Received, may. 08, 2020 - Accepted, Jun. 15, 2020

http://dx.doi.org/10.17268/sel.mat.2020.01.02

Resumen:

En la literatura ecológica, la interferencia (autointerferencia) o competencia entre depredadores para realizar la captura de sus presas, ha sido modelada por diferentes formulaciones matemáticas.

En este trabajo, estableceremos las propiedades dinámicas de modelos de depredación del tipo Leslie-Gower, en el que se incorpora una de estas formas, descrita por la función $b\left ( y\right ) =y^{\alpha }$, con $0<\alpha <1$.

La principal dificultad del análisis se debe a que esta función no es diferenciable para $y=0$, y la matriz Jacobiana no está definida en los puntos de equilibrio sobre el eje horizontal (eje $x$).

Previamente mostramos un resumen con las propiedades fundamentales del modelo del tipo Leslie-Gower, con el objeto de efectuar un adecuado análisis comparativo con modelos donde se incorpora la auto-interferencia entre los depredadores.

Para reforzar nuestros resultados se muestran algunas simulaciones numéricas.

Abstract:

In the ecological literature, the interference (self-interference) or competition among predators to effect the harvesting of their prey has been modeled by different mathematical formulations.

In this work, we will establish the dynamical properties of the Leslie-Gower type predation model, in which is incorporated one of these form, described by the function $b\left ( y\right ) =y^{\alpha }$, with $0<\alpha <1$.

The main difficulty of the analysis is due to this function is not differentiable for $y=0$, and the Jacobian matrix is not defined in the equilibrium points over the horizontal axis ($x-axis$).

Previously, we showed a summary with the fundamental properties of the Leslie-Gower type model, in order to carry out an adequate comparative analysis with models where self-interference between predators is incorporated.

To reinforce our results, some numerical simulations are shown.

Introducción

La relación dinámica entre los depredadores y sus presas ha sido y seguirá siendo uno de los temas dominantes en Ecologıa Teórica (Ecologıa Matemática) y particularmente en Dinámica Poblacional. Esto es debido a su existencia universal [8] y porque permite una mejor comprensión del comportamiento de las cadenas alimenticias o redes tróficas [30,38].

El primer modelo depredador-presa fue propuesto por el matemático italiano Vito Volterra [6,7] en una conocida monografıa publicada en 1926 [39], siendo descrito por un sistema no-lineal de Ecuaciones Diferenciales Ordinarias (EDO); este modelo coincidió con un modelo bidimensional para interacciones bioquımicas propuesto anteriormente por el matemático-fısico estadounidense Alfred J. Lotka: por esto el sistema de EDOS,es conocido con el nombre de modelo de Lotka-Volterra [7,30,31,38].

La caracterıstica principal de este primer modelo es que el único punto de equilibrio positivo para toda condición de parámetros es un centro [10,14], es decir, todas las trayectorias son órbitas cerradas concéntricas alrededor de ese punto. Esto implica que los tamaños poblacionales de los depredadores y sus presas oscilarıan permanentemente alrededor de ese punto para cualquier condición inicial [31]. Este comportamiento de las soluciones del sistema, fue bastante cuestionado porque en la naturaleza no se encuentran interacciones depredador-presa en que existan siempre comportamientos periódicos de los tamaños poblacionales de ambas especies, cualquiera sean las cantidades iniciales consideradas [30,38].

A partir del trabajo de Volterra, se sucedieron distintas propuestas para enfrentar y resolver las diversas objeciones formuladas [34]. Una de las primeras propuestas para solucionar algunas de las objeciones al modelo de Lotka-Volterra fue formulada por el biólogo ruso Georgii F. Gause (1910-1986), quien en 1934 [20] propuso un modelo que toma en cuenta la competencia intraespecıfica de las presas [20], reemplazando el crecimiento Malthusiano incorporado en el modelo de Lotka-Volterra, por la función de crecimiento logıstico (también llamada de ecuación de Verhulst-Pearl [6]).

Un resultado importante es el Teorema propuesto por el matemático ruso Andrei N. Kolmogorov en 1936 [7,30], que establece condiciones en el sistema bidimensional de ecuaciones diferenciales no-lineales, para asegurar la existencia de una única solución periódica estable [13] (matemáticamente, un ciclo lımite atractor). Para el fenómeno de oscilaciones de los tamaños poblacionales de presa y depredadores existe suficiente evidencia en la naturaleza [35] (llamados ciclos ecológicamente estables [13]).

Una alternativa diferente es el modelo propuesto por Patrick Holt Leslie en 1948 [28], que no se ajusta al esquema del modelo de Lotka-Volterra [38]. A diferencia de los modelos compartimentados del tipo Gause [16,38], basados en un principio de trasferencia de masa o energıa [8], el modelo de Leslie se caracteriza porque la ecuación de crecimiento de los depredadores es del tipo logıstico [6], considerando implıcitamente que existe competencia o auto-interferencia entre los depredadores

Leslie asumió que la capacidad de soporte del medio ambiente convencional (environmental carrying capacity) $K_{y}$ es proporcional a la abundancia de presas $x$ [30], esto es, se asume que $K_{y}=K(x)=nx$, donde $x$ representa el tamaño poblacional de las presas. Este supuesto es también asumido en el modelo de Leslie derivado, llamado modelo de May-Holling-Tanner [4,18,30,33,38].

Sin embargo, cuando el depredador es generalista y no existen presas disponibles, los depredadores pueden cambiar a una fuente de alimento diferente [3,25,37]. En este caso, la capacidad de soporte del medio ambiente de los depredadores se puede expresar como $K(x)=nx+c$, donde el parámetro $c>0$, indica la cantidad de alimento alternativo disponible para los depredadores [37]. Obtenemos ası, un modelo de Leslie-Gower modificado [5].

Por otra parte, la acción de los depredadores en la interacción se denomina respuesta funcional de los depredadores o función de consumo [16,30]; se refiere al cambio en la densidad de presa atacada por unidad de tiempo por depredador, cuando la densidad de presa cambia [16]. Se clasifican en varios tipos, dependiendo del tamaño de la población de presas o del tamaño de ambas poblaciones.

El Ecólogo canadiense Crawford S. Holling (1930-2019) en 1959 [26] describió tres tipos de funciones saturadas, basadas en experimentos realizados en laboratorio, las cuales consideró sólo dependientes del tamaño de la población de presas (conocidas como respuesta funcional presa-dependiente).

Posteriormente, Robert J. Taylor en 1984 [36] propuso la respuesta funcional no monotónica, usualmente usada para describir el fenómeno ecológico denominado formación de grupo de defensa [23,37]. Las respuestas funcionales presa-dependiente se nombran como respuestas funcionales del tipo Holling, I, II, IV [36] y existen diversas formas matemáticas para describirlas [30,38].

También han sido propuestas respuestas funcionales que son dependientes de ambas poblaciones tales como la de Beddington-DeAngelis [21,40], las razón-dependiente [38], etc.

En el modelo de Leslie [29], la respuesta funcional del depredador es expresada por la función lineal $h(x)= qx$, también usada en el clásico modelo de Lotka-Volterra [30,38]. Sin embargo, esta respuesta funcional no corresponde a ninguno de los tipos propuestos por Holling [26], por no ser acotada.

En este trabajo estudiaremos el modelo de Leslie-Gower considerando el fenómeno denominado interferencia o competición entre los depredadores (CED), lo cual significa que dos o más depredadores (de la misma especie) se molestan al capturar las presas y por eso las presas pueden evitar la depredación [16].

Una situación similar se puede asumir en modelos Bioeconómicos de acceso abierto, cuando los pescadores se dirigen a un mismo lugar con sus botes o embarcaciones pesqueras para capturar un cardumen de peces; los botes pueden chocar entre ellos [11] e incluso originarse agresiones entre las tripulaciones (efectuando disparos, rotura de implementos de pesca, etc.), como ha sucedido en diversos lugares del planeta [11,12].

Existen tres formas para representar matemáticamente la CED, siendo una de ellas la respuesta funcional de Beddington-DeAngelis [21,38]. Asumiremos la expresión formulada por el matematico canadiense Herbert I. Freedman en 1979 [17], modificando el supuesto de los modelos habituales. Asumió un cambio en la función de depredación proponiendo la función [17]:

$\qquad \qquad B\left( x,y\right) =$ $h\left( x\right) y^{\alpha }$, con $% 0<\alpha <1$.

donde $\alpha$ es la constante de interferencia mutua, y $h\left( x\right)$ es la respuesta funcional del depredador dependiente sólo de la población de presas.

En el siguiente gráfico vemos el comportamiento de la función $% b\left( y\right) =y^{\alpha }$, con $0<\alpha <1$, para diferentes valores de $\alpha$.

Figura: Graficos de la función $b\left( y\right) =y^{\protect\alpha }$, para $\protect\alpha =0.25$, $0.5$, $0.75$ y $\protect\alpha =1$.

Se observa que si $\alpha$ $\rightarrow 0$, la gráfica de la función crece en forma mas abrupta. Esto significa que la competición entre los depredadores es más intensa, mientras que para valores grandes de $% \alpha$ $\rightarrow \infty$ la competición va disminuyendo. Cuando $% \alpha =1$, se tiene el modelo original de Leslie-Gower [29].

Este trabajo está enfocado en el análisis de un modelo derivado del modelo de Leslie-Gower, asumiendo que los depredadores compiten entre ellos y no disponen de un alimento alternativo.

El resto del artıculo está organizado como sigue: En la Sección 2 haremos la formulación de los modelos en estudio y las propiedades básicas del modelo de Leslie-Gower. En la sección 3 mostraremos las propiedades principales del modelo asumiendo CED. En la Sección 4 presentaremos algunas simulaciones computacionales y en la Sección 5 discutiremos las consecuencias de la CED en el modelo del tipo Leslie-Gower modificado

Preliminares

En esta sección, presentaremos las propiedades fundamentales del modelo de Leslie (o de Leslie-Gower) propuesto en 1948 [28] y el modelo de Leslie-Gower modificado considerando alimento alternativo para los depredadores.

En adelante denotaremos por $x=x(t)$ e $y=y(t)$ los tamaños poblacionales de presas y depredadores, respectivamente, medidos en número de individuos, biomasa o densidad por unidad de área o volumen, en cualquier instante $t\geq 0$. Además, $x\left( 0\right) >0$ e $y\left( 0\right) >0$ (No hay generación espontánea en el medio ambiente).

Modelo de Leslie

Es descrito por el sistema de ecuaciones diferenciales no lineales del tipo Kolmogorov [16,22]:

$$L_{\nu }\left( x,y\right) :\left\{ \begin{array}{rcl} \frac{d... ... \frac{dy}{dt} & = & s(1-\frac{y}{nx})y\end{array}% \right. ,$$ (1)

Los parámetros son todos positivos, esto es, $\mu =(r,k,q,s,n)\in \mathbb{R}_{+}^{5}$, teniendo los siguientes significados ecológicos:

$\allowbreak$

$r$ es la tasa intrınseca del crecimiento de las presas;

$K$ es la capacidad de soporte del medio ambiente para las presas;

$q$ es la tasa máxima de consumo per capita, es decir, la máxima cantidad de presas necesarias que pueden ser comidas por cada depredador en cada unidad de tiempo;

$s$ es la tasa intrínsica del crecimiento de los depredadores;

$n$ es la calidad de las presas como alimento para los depredadores.

Propiedades del modelo de Leslie-Gower

Las principales propiedades dinámicas del modelo de Leslie-Gower descrito por el sistema (2.1) o campo vectorial $L_{\nu }\left( x,y\right)$ son:

1. Su dominio es el conjunto

                    $\Omega =\{(x,y)\in \mathbb{R}^{2}: x>0, y\geq 0\}=\mathbb{R}^{+}\times \ \mathbb{R}_{0}^{+}$

   Claramente, el sistema $(2.1)$ no está definido en $x=0$, pero la isoclınica (isoclina) de los depredadores pasa por el punto $\left( 0,0\right)$. Sin embargo estudiaremos su incidencia en la dinámica del sistema (2.1).

2. El conjunto

            $\qquad \Gamma =\left\{ \left( x,y\right) \in %TCIMACRO{\U{211d} }% \mathbb{R} %EndExpansion ^{2}:0\leq x\leq K\text{, }y\geq 0\right\}$,

   es una región positivamente invariante.

3. Las soluciones son acotadas.

           Se puede demostrar usando la compactificación de Poincaré [32] o usando desigualdades diferenciales [9].

   Esta propiedad muestra que el modelo está bien formulado [8].

4. El conjunto

                    $B=\left\{ \left( x,y\right) \in %TCIMACRO{\U{211d} }% \mathbb{R} %EndExpansion ^{2}:0\leq x\leq K\text{, }0\leq y\leq nx\right\}$,

   es también una región positivamente invariante.

           Notamos que si $y>nx$, entonces, $\frac{dy}{dt}<0$. Esto implica que la población de depredadores va a la extinción.

           Si $r\left( 1-\frac{x}{K}\right) -qy>0$, la población de presas crece.

           Si $r\left( 1-\frac{x}{K}\right) -qy\leq 0$, la población de presas se extingue.

5. Los puntos de equilibrio del sistema (2.1) o singularidades del campo vectorial $L_{\nu }\left( x,y\right)$ son $\left( K,0\right)$ y $% \left( x_{e},nx_{e}\right)$ con $x_{e}=\frac{rK}{r+nqK}>0$.

   Por lo tanto $\left( x_{e},nx_{e}\right)$, siempre existe en el interior del primer cuadrante $Int\left( %TCIMACRO{\U{211d} }% \mathbb{R} %EndExpansion _{0}^{+}\right) ^{2}$.

6. La matriz Jacobiana o matriz de la comunidad [38] es

                    $DL_{\nu }\left( x,y\right) =\left( \begin{array}{cc} -\frac{1}{K}\left( 2r... ...y^{2} & -\frac{1}{n}\frac{s}{x}\left( 2y-nx\right) \end{array}% \right)$.

   Notamos que la matriz Jacobiana no está definida para $x=0$.

7. La singularidad $\left( K,0\right)$ es punto silla hiperbólico para cualquier conjunto de parámetros.

           La matriz Jacobiana evaluada en el punto $\left( K,0\right)$ es

                    $DL_{\nu }\left( K,0\right) =\left( \begin{array}{cc} -r & -qK \\ 0 & s\end{array}% \right)$.

   y det $DL_{\nu }\left( K,0\right) =-rs<0$.

   Luego, aplicando el Teorema de la traza y el determinante [15,31], se obtiene la tesis.

8. Naturaleza del punto $\left( 0,0\right)$

   Aunque el modelo de Leslie-Gower no está definido en el punto $\left( 0,0\right)$, estudiaremos su influencia en la dinámica del modelo. Mediante una reparametrización y un reescalamiento del tiempo (variable independiente), estudiaremos un sistema topológicamente equivalente.

   Sean $x=Ku$, $y=nx$ y $\tau =rt$. Después de un proceso algebraico se obtiene el nuevo sistema

   $$V_{\eta }\left( u,v\right) :\left\{ \begin{array}{rcl} \frac{... ...{dv}{d\tau } & = & S\left( u-v\right) v\end{array}% \right. ,$$ (2)

   
   

   donde $\eta =\left( Q,S\right) \in \mathbb{R}_{+}^{2}$, con $Q=\frac{qnK}{r}$ y $S=\frac{s}{r}$. El sistema (2.2) o campo vectorial $V_{\eta }\left( u,v\right)$ es una extensión continua del sistema (2.1).

   Los puntos de equilibrio son $\left( 0,0\right)$, $\left( 1,0\right)$ y $% \left( u_{e},u_{e}\right)$ con $u_{e}$ satisface la ecuación de la isoclinas.

   Obtenemos la ecuación polinomial

   $\qquad 1-u-Qu=0$, es decir, $u_{e}=\frac{1}{Q+1}$.

   La matriz Jacobiana del sistema (2.2) es

            $DV_{\eta }\left( u,v\right) =\left( \begin{array}{cc} -u\left( 3u+2Qv-2\ri... ...\allowbreak Sv & \allowbreak S\left( u-2v\right)% \end{array}% \right)$

   Luego,

            $DV_{\eta }\left( 0,0\right) =\left( \begin{array}{cc} 0 & 0 \\ 0 & 0\end{array}% \right)$

   Para desingularizar el punto $\left( 0,0\right)$ usaremos el método del blowing-up direccional

   Sean $u=m$ y $v=mp$. Luego,

            $\frac{dm}{d\tau }=\frac{du}{d\tau }$, $\frac{dv}{d\tau }=m\frac{dp}{% d\tau }+p\frac{dm}{d\tau }$

   Reemplazando y factorizando se obtiene

            $V_{\eta }\left( m,p\right) :\left\{ \begin{array}{rcl} \frac{dm}{d\tau } &... ...frac{dp}{d\tau } & = & mp\left( S+m-Sp+Qmp-1\right)% \end{array}% \right. ,$

   Reescalando la variable independiente (el tiempo) por $\lambda =m\tau$, se obtiene

            $\frac{dm}{d\tau }=\frac{dm}{d\lambda }\frac{d\lambda }{d\tau }$, $% \frac{dp}{d\tau }=\frac{dp}{d\lambda }\frac{d\lambda }{d\tau }$.

   Por lo tanto, se obtiene el sistema del tipo de Kolmogorov

   $\qquad \bar{V}_{\eta }\left( m,p\right) :\left\{ \begin{array}{rcl} \frac{... ...ac{dp}{d\lambda } & = & p\left( S+m-Sp+Qmp-1\right)% \end{array}% \right. ,$

   La matriz Jacobiana del campo vectorial $\bar{V}_{\eta }\left( m,p\right)$ es

            $D\bar{V}_{\eta }\left( m,p\right) =\left( \begin{array}{cc} \allowbreak 1-... ...} \\ p\left( Qp+1\right) & \allowbreak S+m-2Sp+2Qmp-1\end{array}% \right)$

   La matriz Jacobiana evaluada en $\left( 0,0\right)$ es

            $D\bar{V}_{\eta }\left( 0,0\right) =\left( \begin{array}{cc} 1 & 0 \\ 0 & S-1\end{array}% \right)$

   Luego, el punto $\left( 0,0\right)$ del campo vectorial $\bar{V}_{\eta }\left( m,p\right)$ es

   I)	una silla hiperbólica, si y sólo si, $S<1$.
   II)	Un nodo repulsor, si y sólo si, $S<1$.
   III)	Una silla-nodo repulsora, si y sólo si, $S=1$.

   Entonces, por el blowing-down, el punto de equilibrio no-hiperbólico $% \left( 0,0\right)$ del campo vectorial $V_{\eta }\left( m,p\right)$ es una silla-nodo repulsor para todo valor de $S$.

   Por tanto, el punto de equilibrio no-hiperbólico $\left( 0,0\right)$ del campo vectorial $V_{\eta }\left( u,v\right)$ es una silla-nodo repulsor para todo valor de $S$.

   En consecuencia, el punto $\left( 0,0\right)$ de la extensión continua del modelo de Leslie-Gower tiene un comportamiento de silla-nodo repulsor

   
   

   El sistema (2.2) se puede utilizar para demostrar todas las propiedades del sistema (2.1).

   Por ejemplo, la matriz Jacobiana del campo vectorial $V_{\eta }\left( u,v\right)$ o sistema (2.2) evaluada en el punto de equilibrio $\left( 1,0\right)$ , equivalente al punto $\left( K,0\right)$ del sistema (2.1), es

            $DV_{\eta }\left( 1,0\right) =\left( \begin{array}{cc} -1 & -Q \\ 0 & S\end{array}% \right)$

   Luego, det $DV_{\eta }\left( 1,0\right) =-S<0$, y el punto $\left( 1,0\right)$ es una silla hiperbólica.

9. La singularidad $\left( x_{e},nx_{e}\right)$ que está siempre al interior del primer cuadrante $Int\left( %TCIMACRO{\U{211d} }% \mathbb{R} %EndExpansion _{0}^{+}\right) ^{2}$ es global asintóticamente estable (g.a.e).

   La matriz Jacobiana del campo vectorial $V_{\eta }\left( u,v\right)$ topológicamente equivalente $\left( x_{e},nx_{e}\right)$ en el punto $\left( u,u\right)$ es

            $DV_{\eta }\left( u,u\right) =\left( \begin{array}{cc} -u\left( 3u+2Qu-2\ri... ... & -Qu^{2} \\ \allowbreak Su & \allowbreak -Su\end{array}% \right)$.

   Luego,

           det $DV_{\eta }\left( u,u\right) =\allowbreak Su^{2}\left( 3\left( Q+1\right) u-2\right)$, y

           tr $DV_{\eta }\left( u,u\right) =-u\left( S+3u+2Qu-2\right)$.

   En el punto $\left( \frac{1}{Q+1},\frac{1}{Q+1}\right)$, se tiene

           det $DV_{\eta }\left( \frac{1}{Q+1},\frac{1}{Q+1}\right) =\allowbreak S\left( \frac{1}{Q+1}\right) ^{2}>0$, y

           tr $DV_{\eta }\left( \frac{1}{Q+1},\frac{1}{Q+1}\right) =\allowbreak -% \frac{S+QS+1}{\left( Q+1\right) ^{2}}<0$.

   Por tanto, el punto de equilibrio $\left( \frac{1}{Q+1},\frac{1}{Q+1}\right)$ es un atractor local.

   Considerando el Teorema de Poincaré-Bendixson válido sólo para sistemas bidimensionales, se tiene que el punto de equilibrio $\left( \frac{1}{Q+1},\frac{1}{Q+1}\right)$ es g.a.e.

   En consecuencia $\left( x_{e},nx_{e}\right)$ es g.a.e., en el sistema (2.1).

   Observamos que la matriz Jacobiana evaluada en el punto $\left( x_{e},nx_{e}\right)$ es

                    $DL_{\nu }\left( x,nx\right) =\allowbreak \left( \begin{array}{cc} \frac{Kr-2rx-Knqx}{K} & -qx \\ ns & -s\end{array}% \right)$,

   
   

   También la estabilidad global se puede probar usando la función de Lyapunov [22] descrita por

                    $V\left( x,y\right) =\ln \left( \frac{x}{x_{e}}\right) +\frac{% x_{e}}{x}+m\left( \ln \left( \frac{y}{y_{e}}\right) +\frac{y_{e}}{y}\right)$,

   con $m>0$, una constante adecuada [27].

   Esta función es muy parecida a una de las propuestas por B-S. Goh en 1980 [22].

10. El modelo de Leslie-Gower no tiene ciclos.

    Se aplica el criterio de Bendixson-Dulac [11,22,32] usando la función $g\left( u,v\right) =\frac{1}{u^{2}v}$ en el sistema (2.2).

Modelo de Leslie-Gower modificado

Se denomina modelo de Leslie-Gower modificado o que se ajusta al esquema de Leslie a un modelo descrito por un sistema EDO generalizado

$$G_{\sigma }\left( x,y\right) :\left\{ \begin{array}{rcl} \fra... ...& s(1-\frac{y}{K_{y}\left( x\right) })y\end{array}% \right. .$$ (3)

donde

$f\left( x\right)$ es la tasa de crecimiento per capita de las presas,

$h\left( x,y\right)$ es la respuesta funcional o tasa de consumo de los depredadores, dependiente sólo de la población de presas (respuesta funcional presa dependiente) o de ambas poblaciones,

$K_{y}\left( x\right)$ es la capacidad de soporte de los depredadores, dependiente del tamaño poblacional de las presas, y

$s$ es la tasa de crecimiento intrınseco de los depredadores.

A partir del sistema (2.3) se pueden generar diversos modelos considerando diversos fenómenos ecológicos [1,18,24,37]

Planteamiento del modelo

En este sección presentaremos el modelo derivado del modelo de Leslie que será analizados en la siguiente sección.

Modelo con CED

El modelo a estudiar es descrito por el sistema autónomo bidimensional de ecuaciones diferenciales no lineales del tipo Kolmogorov [16,22]:

$$X_{\mu }\left( x,y\right) :\left\{ \begin{array}{rcl} \frac{d... ...} & = & s\left( 1-\frac{y}{nx}\right) y\end{array}% \right. ,$$ (4)

donde $x=x\left( t\right)$ e $y=y\left( t\right)$ y los parámetros son todos positivos teniendo los mismos significados ecológicos anteriores. El parámetro $\alpha \in ]0,1[$, indica el nivel de interferencia entre los depredadores. Por lo tanto, $\mu =\left( r,K,q,s,n,\alpha \right) \in \mathbb{R}_{+}^{5}\times ]0,1[$.

Hacemos presente que el sistema (3.1) no está definido en $x=0$, es decir, su dominio es el conjunto

$\qquad \qquad \Omega =\left\{ \left( x,y\right) \in \mathbb{R}^{2}: x>0,\ y\geq 0\right\} =\mathbb{R}^{+}\times \mathbb{R}_{0}^{+}$

Con el objeto de simplificar los cálculos seguimos la metodologıa desarrollada en [1,2,18,19,24]; efectuamos un cambio de variables y un reescalamiento del tiempo como se muestra en la siguiente proposición.

Demonstración. Sea $x=Ku$ e $y=nKv$ entonces $\frac{dx}{dt}=K\frac{du}{dt}$ y $\frac{dy}{dt}% =nK\frac{dv}{dt}$.

Reemplazando y factorizando en el sistema (3.1) tenemos:

$\qquad V_{\mu }\left( u,v\right) :\left\{ \begin{array}{rcl} \frac{du}{dt}... ...\\ \frac{dv}{dt} & = & s\left( 1-\frac{v}{u}\right) v\end{array}% \right.$

Este último sistema no está definido para $u=0$.

Para subsanar esta dificultad, realizamos un reescalamiento del tiempo, dado por $\tau =\frac{r}{u}t$.

Por regla de la cadena obtenemos

         $\frac{du}{dt}=\frac{du}{d\tau }\frac{d\tau }{dt}=\frac{r}{u}\frac{du}{d\tau }$ y $\frac{dv}{dt}=\frac{dv}{d\tau }\frac{d\tau }{dt}=\frac{r}{u}% \frac{dv}{d\tau }$

Reemplazando en ell campo vectorial $V_{\mu }\left( x,y\right)$ tenemos

$\qquad \bar{V}_{\mu }\left( x,y\right) :\left\{ \begin{array}{rcl} \frac{d... ...\\ \frac{dv}{dt} & = & \frac{s}{r}\left( u-v\right) v\end{array}% \right.$

Definiendo $Q=\frac{q\left( nK\right) ^{\alpha }}{r}$ y $S=\frac{s}{r}$ obtenemos el sistema (3.2). $\qedsymbol$

Nota 1   Hemos construido la función $\varphi :\overline{\Omega }\times \mathbb{R}% \longrightarrow {\Omega \times \mathbb{R}}$ dada por

$\qquad \varphi \left( u,v,\tau \right) =\left( Ku.nKv,\frac{u}{r}\tau \right) =\left( x,y,t\right)$

La matriz Jacobiana de la transformación $\varphi$ es

$\qquad D\varphi (u,v,\tau )=\left( \begin{array}{ccc} K & 0 & 0 \\ 0 & nK & 0 \\ \frac{1}{r} & 0 & \frac{u}{r}% \end{array}% \right)$

y det $D_{\varphi }=\frac{nK^{2}u}{r}>0$.

Entonces, $\varphi$ es un difeomorfismo que preserva la orientación del tiempo, por lo cual el campo vectorial $X_{\mu }$ es topológicamente equivalente al campo vectorial $U_{\eta }=\varphi \circ X_{\mu }$.

De esta forma, el sistema (3.2) es una extensión continua del sistema (3.1), para $x=0$, en particular en el punto $(0,0)$.

Puntos de equilibrio

Los puntos crıticos (estacionarios o de equilibrio) del sistema se obtienen considerando $\frac{du}{d\tau }=0=\frac{dv}{d\tau }$, los cuales corresponden a las soluciones de las ecuaciones:

$\qquad \left\{ \begin{array}{rcl} u^{2}\left( 1-u-Qv^{\alpha }\right) & = & 0 \\ S\left( u-v\right) v & = & 0\end{array}% \right.$

De esta forma, los puntos de equilibrio sobre los ejes son $(0,0)$, $(1,0)$ y $% \left( u_{e},v_{e}\right)$ que satisface la ecuaciones de las isoclinas (nulclinas) $v=u$ y $v=\left( \frac{1-u}{Q}\right) ^{\frac{1}{\alpha }}$.

Demonstración. Como el punto de equilibrio positivo $\left( u_{e},u_{e}\right)$, es el punto de intersección de las isoclinas, la ordenada $u_{e}$ satisface la relación

$\qquad u_{e}+Q{u_{e}}^{\alpha }-1=0$.

Luego, $u_{e}$ es un cero (o raız) de la función $f\left( u\right) =u+Qu^{\alpha }-1$ [19].

Observemos que $f\left( u\right)$ es una función continua para todo $% u\geq 0$.

Como $f\left( 1\right) =Q>0$ y $f(0)=-1<0$ entonces, por el Teorema del Valor intermedio, $f\left( u\right)$ tiene un cero en el intervalo abierto $% \left] 0,1\right[$.

Además, la derivada de $f$ con respecto a $u$ es $f\left( u\right) =1+Q{% \alpha }\frac{1}{u^{1-\alpha }}>0$ para $u>0$. Por lo tanto, $f\left( u\right)$ es estrictamente creciente para todo $u>0$.

De esta forma, existe un único cero para la función $f$para todo $% u>0$, y por consecuencia, existe un único punto de equilibrio positivo $% \left( u_{e},u_{e}\right)$ para el sistema (3.2) al interior del primer cuadrante. $\qedsymbol$

Matriz Jacobiana

Para analizar la naturaleza de cada uno de los puntos de equilibrio, de acuerdo al Teorema de Hartman y Grossman (o Teorema de Linealización) [10,15], se debe analizar la matriz Jacobiana [4] (o matriz de la comunidad [38]) la cual es

$$DU_{\eta }(u,v)=\left( \begin{array}{cc} u\left( 2-3u-2Qv^{\... ...ha -1} \\ -Sv^{2} & S\left( u-2v\right)% \end{array}% \right)$$

Notamos que $DU_{\eta }(u,v)$ no está definida para $v=0$.

Resultados Principales

Para el sistema (3.2) tenemos las siguientes propiedades:

Demonstración. Como el modelo es del tipo Kolmogorov [16], entonces los ejes coordenados son conjuntos invariantes [10].

Considerando $u=1$ se tiene que $\frac{du}{d\tau }=-Qv^{\alpha }<0$ y cualquiera sea el signo de $\frac{dv}{d\tau }$, las trayectorias con condiciones iniciales en el plano de fase, fuera de la región $\overline{% \Gamma }$ entran en ella, mientras las que tienen condiciones iniciales dentro de $\overline {\Gamma }$, permanecen allı $\overline {\Gamma }$.

Por lo tanto, se tiene que $\overline {\Gamma }$ es una región de invariancia para el sistema (3.2). $\qedsymbol$

En la siguiente figura se muestra la región de invarianza $\overline{% \Gamma }$.

Figura 4.1: Región de invariancia $\overline {\Gamma }$.

Demonstración. Usando la compactificación de Poincaré [32].

Sean $X=\frac{u}{v}$ e $Y=\frac{1}{v}$, entonces,

$\qquad \frac{dX}{d\tau }=\frac{1}{v}\frac{du}{d\tau }-\frac{u}{v^{2}}\frac{% ... ...d\tau }=\frac{1}{v^{2}}\left( v\frac{du}{d\tau }-u\frac{dv}{d\tau }% \right)$

$\qquad \frac{dY}{d\tau }=-\frac{1}{v^{2}}\frac{dv}{d\tau }$

Luego, $u=\frac{X}{Y}$ y $v=\frac{1}{Y}$ obteniendo

$\qquad \left\{ \begin{array}{rcl} \frac{dX}{d\tau } & = & Y\left( \frac{du... ...\ \frac{dY}{d\tau } & = & -Y^{2}\frac{dv}{d\tau }% \end{array}% \right.$

Luego, se tiene

$\qquad U_{\eta }\left( X,Y\right) :\left\{ \begin{array}{rcl} \frac{dX}{d\... ...\frac{X}{Y}-\frac{1}{Y}\right) \frac{1}{Y}\right]% \end{array}% \right.$

Simplificando obtenemos

$$\bar{U}_{\eta }\left( X,Y\right) :\left\{ \begin{array}{rcl} ... ...dY}{d\tau } & = & -S\left( X-1\right)% \end{array}% \right.$$ (5)

Este último sistema no está definido para $Y=0.$ Para solucionar esta dificultad, se considera un reescalamiento del tiempo dado por $T=\frac{% 1}{Y^{2}}\tau$. Entonces, por regla de la cadena

$\qquad \frac{dX}{d\tau }=\frac{dX}{dT}\frac{dT}{d\tau }=\frac{1}{Y^{2}}% \frac{dX}{dT}$ y $\frac{dY}{d\tau }=\frac{dY}{dT}\frac{dT}{d\tau }=% \frac{1}{Y^{2}}\frac{dY}{dT}$.

Aplicando estas relaciones al sistema $\left(4.1\right)$ se obtiene

$\qquad \widetilde{U}_{\eta }\left( X,Y\right) :\left\{ \begin{array}{rcl} ... ...] \\ \frac{dY}{dT} & = & -SY^{2}\left( X-1\right)% \end{array}% \right.$

$\frac{d}{dY}\left( X\left( \left( Y-X-QY^{1-\alpha }\right) -SY\left( X-1\rig... ...k X\left( S-\frac{Q}{Y^{\alpha }}-SX+% \frac{Q}{Y^{\alpha }}\alpha +1\right)$

La matriz jacobiana del campo vectorial $\widetilde{U}_{\eta }$ es

$$D\widetilde{U}_{\eta }\left( X,Y\right) =\left( \begin{array}{... ...ght) \\ -SY^{2} & - 2SY\left( X-1\right)% \end{array}% \right)$$

La cual no está definida para $Y=0$.

Aplicando el método de blowing up direccional, mediante el cambio de variables $X=rs^{2}, Y=s$ entonces

$\qquad \frac{dX}{dT}=2sr\frac{ds}{dT}$ y $\frac{dY}{dT}=\frac{ds}{dT}$

De esta forma, obtenemos el sistema

$$\overline{U}_{\eta }\left( r,s\right) :\left\{ \begin{array}{r... ... \right) \\ \frac{ds}{dT} & = & \frac{dY}{dT}% \end{array}% \right.$$

después de algunas manipulaciones algebraicas llegamos a

$$\overline{U}_{\eta }\left( r,s\right) :\left\{ \begin{array}{... ...frac{ds}{dT} & = & Ss^{2}(1-rs^{2})% \end{array}% \right. ,$$ (6)

Hacemos un nuevo reescalamiento de la variable independiente (el tiempo) dado por $\lambda =sT$, entonces

         $\frac{dr}{dT}=\frac{dr}{d\lambda }\frac{d\lambda }{dT}=s\frac{dr}{% d\lambda }$ y $\frac{ds}{dT}=\frac{ds}{d\lambda }\frac{d\lambda }{dT}% =s\frac{ds}{d\lambda }$.

Aplicando estas relaciones al sistema (4.2) obtenemos

$$\overline{U}_{\eta }\left( r,s\right) :\left\{ \begin{array}{... ...da } & = & Ss\left( 1-rs^{2}\right)% \end{array}% \right. ,$$ (7)

La matriz Jacobiana asociada al sistema del tipo Kolmogorov (4.3) que es

$$D\overline{U}_{\eta }(r,s)=\left( \begin{array}{cc} rs\left(... ...-\alpha )rs^{1-\alpha } \\ -Ss^{3} & -3rs^{2}+S\end{array}% \right)$$   ,$$$$

Luego,

$$D\overline{U}_{\eta}(0,0) =\left( \begin{array}{cc} -S & 0 \\ 0 & S\end{array}% \right)$$   ,$$$$

Por consiguiente, el punto $(0,0)$ del campo vectorial $\overline{U}_{\eta }\left( r,s\right)$ es una silla hiperbólica. Luego es una silla no-hiperbólica del campo vectorial $\widetilde{U}_{\eta }\left( X,Y\right)$. Entonces, el punto $\left( 0,\infty \right)$ es una silla no-hiperbólica del campo vectorial compactificado de $U_{\eta }\left( u,v\right)$.

Por tanto, las trayectorias son acotadas. $\qedsymbol$

Naturaleza de los puntos de equilibrio sobre los ejes

La naturaleza de los puntos sobre los ejes no puede ser determinadas por los métodos usuales en el sistema (3.2), porque la matriz Jacobiana no está definida para $v=0$.

Para el estudio de estos puntos haremos un nuevo cambio de variables en el sistema (3.2).

        Sea $v=z^{\frac{1}{\alpha }}$, luego $\frac{dv}{d\tau }=\frac{1}{% \alpha }z^{\frac{1}{\alpha }-1}\frac{dz}{d\tau }$

Claramente si $v=0$, entonces $z=0$. Además, $\frac{d}{dz}\left( z^{% \frac{1}{\alpha }}\right) =\allowbreak \frac{1}{\alpha }z^{\frac{1}{\alpha }% -1}$ está definida para todo $0<\alpha <1$.

Reemplazando tenemos que

$\qquad Z_{\eta }\left( u,z\right) :\left\{ \begin{array}{rcl} \frac{du}{d\... ...{% \frac{1}{\alpha }}\right) z^{\frac{1}{\alpha }}% \end{array}% \right.$

o bien,

$$Z_{\eta }\left( u,z\right) :\left\{ \begin{array}{rcl} \frac{... ...\left( u-z^{\frac{1}{\alpha }}\right) z\end{array}% \right. ,$$ (8)

donde $\eta =\left( Q,S,\alpha \right) \in \mathbb{R}_{+}^{2}\times ]0,1[$. Su dominio es:

         $\hat{\Omega}=\left\{ \left( u,z\right) \in \mathbb{R}^{2}: u\geq 0, z\geq 0\right\} = \left( \mathbb{R}_{0}^{+}\right) ^{2}$.

Los puntos de equilibrio sobre los ejes son $(0,0),(1,0)$. Además, el punto $\left( u_{e},z_{e}\right)$ debe satisfacer las ecuaciones de las isoclinas (nulclinas),

$\qquad z=\frac{1}{Q}\left( 1-u\right)$ recta de pendiente negativa que pasa por los puntos $\left( 0,\frac{1}{Q}\right)$ y $\left( 1,0\right)$ y,

         $z=u^{\alpha }$, función creciente no acotada, por lo que se puede asegurar la existencia de un punto de intersección $\left( u_{e},z_{e}\right)$, con $0<u_{e}<1$.

La matriz Jacobiana es

$\qquad DZ_{\eta }\left( u,z\right) =\left( \begin{array}{cc} -u\left( 3u+2... ...{\alpha }}-\alpha u+z^{\frac{1}{% \alpha }}\right) \end{array}% \right)$,

definida para todo $\left( u,z\right) \in \hat{\Omega}$.

Para el sistema (4.4) tenemos las siguientes propiedades

Demonstración. Claramente,

         $DZ_{\eta }\left( 0,0\right) =\left( \begin{array}{cc} 0 & 0 \\ 0 & 0\end{array}% \right)$

Consideremos el blowing-up direccional [10,14] dado por $u=r$ y $z=rs$, donde $r$ y $s$ son funciones de $\tau$. Se tiene $% \frac{dz}{d\tau }=s\frac{dr}{d\tau }+r\frac{ds}{d\tau }$. Luego, $\frac{ds}{d\tau }% =\frac{1}{r}\left( \frac{ds}{d\tau }-s\frac{dr}{d\tau }\right)$

Reemplazando en el sistema $\left( 4.4\right)$ se tiene el nuevo sistema

$\qquad \bar{Z}_{\eta }\left( r,s\right) :\left\{ \begin{array}{rcl} \frac{... ...pha }}\right) rs-sr^{2}\left( 1-r-Qrs\right) \right) \end{array}% \right.$

factorizando obtenemos:

$\qquad \bar{Z}_{\eta }\left( r,s\right) :\left\{ \begin{array}{rcl} \frac{... ...% \frac{1}{\alpha }-1}s^{\frac{1}{\alpha }}\right) \end{array}% \right.$

Reescalando el tiempo (variable independiente) $T=r\tau$, se tiene:

$\qquad \hat{Z}_{\eta }\left( r,s\right) :\left\{ \begin{array}{rcl} \frac{... ... r^{\frac{1}{\alpha }-1}s^{\frac{1}{\alpha }}\right) \end{array}% \right.$

el cual es un sistema del tipo Kolmogorov y cuya matriz Jacobiana es:

$\qquad D\hat{Z}_{\eta }\left( r,s\right) =\left( \begin{array}{cc} \allowb... ...a }-2}s^{\frac{1}{\alpha }}\alpha +1\right) & S\alpha -1\end{array}% \right)$.

Evaluada en el punto $\left( 0,0\right)$ se obtiene

$\qquad D\hat{Z}_{\eta }\left( 0,0\right) =\left( \begin{array}{cc} 1 & 0 \\ 0 & S\alpha -1\end{array}% \right)$.

Luego, $\left( 0,0\right)$ del campo vectorial $\hat{Z}_{\eta }\left( r,s\right)$ es un punto:

i) repulsor hiperbólico, si y sólo si, $S\alpha -1>0$,

ii) silla-nodo hiperbólico, si y sólo si, $S\alpha -1<0$, repulsor por el $eje-x$, y

iii) silla-nodo no-hiperbólico, si y sólo si, $S\alpha -1=0$, repulsor por el $eje-x$.

Por lo tanto, $\left( 0,0\right)$ del campo vectorial $\bar{Z}_{\eta }\left( r,s\right)$ es silla-nodo no-hiperbólico repulsor, para todo $% S>0$ y $0<\alpha <1$.

Por el blowing-down, $\left( 0,0\right)$ del campo vectorial $V_{\eta }\left( u,v\right)$ es silla-nodo no-hiperbólico repulsor, para todo $% \eta =\left( Q,S,\alpha \right) \in \mathbb{R}_{+}^{2}\times ]0,1[$. $\qedsymbol$

Demonstración. La matriz Jacobiana evaluada en $\left( 1,0\right)$ es

$\qquad DZ_{\eta }\left( 1,0\right) =\left( \begin{array}{cc} -1 & -Q \\ 0 & S\alpha \end{array}% \right)$

Luego, el punto de equilibrio $\left( 1,0\right)$ del sistema (4.4) o singularidad del campo vectorial $Z_{\eta }\left( u,z\right)$ es una silla hiperbólica para todo $\eta =\left( Q,S,\alpha \right) \in \mathbb{R}_{+}^{2}\times ]0,1[$.

En consecuencia, el punto de equilibrio $\left( 1,0\right)$, del sistema (3.2) es una silla no-hiperbólica, para todo $\eta =\left( Q,S,\alpha \right) \in \mathbb{R}_{+}^{2}\times ]0,1[$. $\qedsymbol$

Naturaleza del punto de equilibrio positivo

Recordamos que el punto de equilbrio positivo, $\left( u_{e},v_{e}\right)$ del campo vectorial $U_{\eta }(u,v)$ o sistema (3.1) satisface las siguientes relaciones:

$$\left\{ \begin{array}{rcl} u_e & = & v_e \\ Q{v_e}^{\alpha} & = & 1-u_e\end{array}% \right.$$

Demonstración. Evaluando la matriz jacobiana en el punto $(u_{e},v_{e})$ y usando la relación $u_{e}=v_{e}$ obtenemos:

$\qquad DU_{\eta }\left( u_{e},v_{e}\right) =\left( \begin{array}{cc} u_{e}... ...{e}}% ^{1+\alpha } \\ S{u_{e}} & -Su_{e}% \end{array}% \right)$

Finalmente, aplicando la relación $Q{u_{e}}^{\alpha }=1-u_{e}$, se obtiene la matriz

$\qquad DU_{\eta }\left( u_{e},v_{e}\right) =\left( \begin{array}{cc} -{u_{... ...pha {u_{e}}(1-u_{e}) \\ S{u_{e}} & -Su_{e}% \end{array}% \right)$,

obteniendo que

        det $DU_{\eta }\left( u_{e},v_{e}\right) =S{u_{e}}^{3}+\alpha S{u_{e}}% ^{2}\left(... ...{e}\right) =\allowbreak Su_{e}^{2}\left( \alpha +u_{e}-\alpha u_{e}\right) >0$, pues $1-u_{e}>0$.

Por otra parte,

        tr $DU_{\eta }\left( u_{e},v_{e}\right) =-{u_{e}}^{2}-Su_{e}=-u_{e}% \left( u_{e}+S\right) <0$.

De esta forma, $\left( u_{e},v_{e}\right)$ es local asintóticamente estable. $\qedsymbol$

Nota 2   Recordando que el punto $\left( 0,0\right)$ y el punto $\left( 1,0\right)$ son silla-nodo repulsor, por Teorema de Poincaré-Bendixson el único punto de equilibrio positivo $\left( u_{e},v_{e}\right)$ es globalmente asintóticamente estable (gae). Sin embargo, reforzaremos esta conclusión con los siguientes teoremas.

Demonstración. Consideremos la función de Lyapunov [22,27] dada por

$$V\left( u,v\right) =S\int_{u_{e}}^{u}\frac{\xi -u_{e}}{{\xi }^{2}... ...{v_{e}}^{v}\frac{{\omega }^{\alpha }-{v_{e}}^{\alpha }}{\omega }% d\omega$$

Definamos

$\qquad h\left( u\right) =S\int_{u_{e}}^{u}\frac{\xi -u_{e}}{{\xi }^{2}}d\xi$ y $j\left( v\right) =Q\int_{v_{e}}^{v}\frac{{\omega }^{\alpha }-{v_{e}}^{\alpha }}{\omega }d\omega$

Entonces $V\left( u,v\right) =h\left( u\right) +j\left( v\right)$

Si $u_{e}<u$ entonces $h\left( u\right) >0$ puesto que $\frac{\xi -u_{e}}{{% \xi }^{2}}>0$ para todo $\xi \in (u_{e},u)$.

Si $u<u_{e}$ entonces $h(u)=-S\int_{u}^{u_{e}}\frac{\xi -u_{e}}{{\xi }^{2}}% d\xi$ y como $\frac{\xi -u_{e}}{{\xi }^{2}}<0$ para todo $\xi \in (u,u_{e})$.

Entonces, $-\frac{u_{e}-\xi }{{\xi }^{2}}>0$ para todo $\xi \in (u,u_{e})$.

Luego, por la linealidad de la integral $h\left( u\right) =S\int_{u}^{u_{e}}% \frac{u_{e}-\xi }{{\xi }^{2}}d\xi >0$.

Por lo tanto, $h\left( u\right) >0$ para todo $u\neq u_{e}$.

Notemos que $h\left( u\right)$ es la función nula si y sólo si $% u=u_{e}$

De forma análoga ocurre con $j\left( v\right)$.

Si $v_{e}<v$ entonces $j\left( v\right) >0$ puesto que $\frac{{\omega }% ^{\alpha }-{v_{e}}^{\alpha }}{\omega }>0$ para todo $\omega \in \left( v_{e},v\right)$.

Si $v<v_{e}$ entonces $j\left( v\right) =-Q\int_{v}^{v_{e}}\frac{{\omega }% ^{\alpha }-{v_{e}}^{\alpha }}{\omega }d\omega$ y como $\frac{{\omega }% ^{\alpha }-{v_{e}}^{\alpha }}{\omega }<0$ entonces $j\left( v\right) >0.$

Por lo tanto, $j\left( v\right) >0$ para todo $v\neq v_{e}$, y $j\left( v\right)$ es la función nula si y sólo si $v=v_{e}$.

Ası, $V\left( u,v\right) >0$ para todo $\left( u,v\right) \neq \left( u_{e},v_{e}\right)$ y $V(u,v)=0$ si y sólo si $(u,v)=\left( u_{e},v_{e}\right)$.

Luego, $V(u,v)$ es definida positiva en una vecindad perforada de $\left( u_{e},v_{e}\right)$ al interior del primer cuadrante.

Por otra parte, la derivada de $V(u,v)$ respecto al tiempo es

$$\frac{dV\left( u,v\right) }{d\tau }=S\left( \frac{u-u_{e}}{u^{2}}... ...left( \frac{% v^{\alpha }-{v_{e}}^{\alpha }}{v}\right) v\left( u-v\right)$$

$=S\left( u-u_{e}\right) \le... ...ha }\right) +SQ\left( v^{\alpha }-{v_{e}}^{\alpha }\right) \left( u-v\right)$

Usando las ecuaciones que satisface el punto de equilibrio $\left( u_{e},v_{e}\right)$

$\qquad 1-u_{e}-Qv^{\alpha }=0$ y $u_{e}-v_{e}=0$

tenemos que

$\frac{dV\left( u,v\right) }{d\tau }=S\left( u-u_{e}\right) \left( 1-u-Qv^{\al... ...lpha }-{v_{e}}^{\alpha }\right) \left( u-v-\left( u_{e}-v_{e}\right) \right)$

$\qquad =S(u-u_{e})(-(u-u_{e})-Q(v^{\alpha }-{v_{e}}^{\alpha }))+SQ(v^{\alpha }-{v_{e}}^{\alpha })((u-u_{e})-(v-v_{e}))$

$\qquad =-S(u-u_{e})^{2}-SQ(u-u_{e})(v^{\alpha }-{v_{e}}^{\alpha })+SQ(v^{\alpha }-{v_{e}^{\alpha }})(u-u_{e})-SQ(v^{\alpha }-{v_{e}}^{\alpha })(v-v_{e})$

$\qquad =-S(u-u_{e})^{2}-SQ(v^{\alpha }-{v_{e}}^{\alpha })(v-v_{e})$

Como $0<\alpha <1$, entonces existe un función $p\left( v,v_{e}\right)$ tal que

$$\left( v-v_{e}\right) =\left( v^{\alpha }-{v_{e}}^{\alpha }\right) p\left( v,v_{e}\right)$$

donde $p\left( v,v_{e}\right) )$ es una función de producto de potencias de $v$ y $v_{e}$ con coeficientes positivos.

Luego,

$$\frac{dV\left( u,v\right) }{d\tau }=-S\left( u-u_{e}\right) ^{2}-SQ\left( v^{\alpha }-{v_{e}}^{\alpha }\right) ^{2}p\left( v,v_{e}\right)$$

Ası, $\frac{dV\left( u,v\right) }{d\tau }\leq 0$ para todo $\left( u,v\right)$ en el primer cuadrante.

Ahora bien, como

$\qquad \qquad \left( u-u_{e}\right) ^{2}>0$ y $Q\left( v^{\alpha }-{v_{e}}% ^{\alpha }\right) ^{2}p\left( v,v_{e}\right) >0$.

$\qquad \frac{dV\left( u,v\right) }{d\tau }=0 \Leftrightarrow -\left( u-u_{... ...-Q\left( v^{\alpha }-{v_{e}}^{\alpha }\right) ^{2}p\left( v,v_{e}\right) =0$

$\qquad \Leftrightarrow \left( u-u_{e}\right) ^{2}=0 \wedge \left( v^{\alpha }-{v_{e}}^{\alpha }\right) ^{2}p\left( v,v_{e}\right) =0$

Notemos que $p\left( v,v_{e}\right)$ es una función positiva puesto que $v>0$ y $v_{e}>0$.

Luego, como cada sumando de $p\left( v,v_{e}\right)$ es positivo, no existe una raiz positiva de $p(v,v_{e})$.

Por lo tanto, para $v>0$ tenemos

$\qquad \left( v^{\alpha }-{v_{e}}^{\alpha }\right) ^{2}p\left( v,v_{e}\right)... ...0 \Leftrightarrow v^{\alpha }={v_{e}}^{\alpha }\ \Leftrightarrow v=v_{e}$

Entonces, al interior del primer cuandrante, $\frac{dV\left( u,v\right) }{% d\tau }$ se anula si y sólo si $\left( u,v\right) =\left( u_{e},v_{e}\right)$.

Luego, $\frac{dV\left( u,v\right) }{d\tau }$ es definida negativa en una vecindad perforada de $\left( u_{e},v_{e}\right)$ al interior del primer cuadrante.

Finalmente, por el teorema de Lyapunov [32], $\left( u_{e},v_{e}\right)$ es global asintóticamente estable (g.a.e). $\qedsymbol$

Demonstración. Consideremos la función $g\left( u,v\right) =\frac{1}{u^{2}v}$.

Observemos que $g\left( u,v\right)$ es continuamente diferenciable al interior del primer cuadrante.

Ahora bien, considerando

$\qquad \frac{d}{du}\left[ P\left( u,v\right) g\left( u,v\right) \right] +% \frac{d}{dv}\left[ M\left( u,v\right) g\left( u,v\right) \right]$

con $P\left( u,v\right) =u^{2}\left( 1-u-Qv^{\alpha }\right)$ y $M\left( u,v\right) =Sv\left( u-v\right)$, se tiene que:

                 $\frac{d}{du}\left[ \frac{u^{2}\left( 1-u-Qv^{\alpha }\right) }{u^{2}v}\right] +\frac{d}{dv}\left[ \frac{Sv\left( u-v\right) }{u^{2}v}% \right]$

$\qquad \qquad =\frac{d}{du}\left[ \frac{1-u-Qv^{\alpha }}{v}\right] +\frac{d}{dv}\left[ \frac{S\left( u-v\right) }{u^{2}}\right]$

$\qquad \qquad =-\frac{1}{v}-\frac{S}{u^{2}}<0$,

para todo $\left( u,v\right)$ al interior del primer cuadrante.

Por lo tanto, por el criterio de Bendixson-Dulac [10,32], no existen ciclos al interior del primer cuadrante. $\qedsymbol$

Simulaciones numéricas

Para reforzar los resultados analıticos obtenidos para el sistema (3.2) se presentarán algunas simulaciones, indicando los valores de los parámetros utilizados.

Figura 5.1: Para los valores de parámetros $Q= \frac {1}{\sqrt {2}}$, $S=% \frac{1}{\sqrt{2}}$, $\alpha =\frac {1}{\sqrt {2}}$, el punto $(u_e,v_e) \approx (0.5416,0.5416)$ foco estable.

Figura 5.2: Para los valores de parámetros $Q=\frac {30}{\sqrt {2}}$. $S=% \frac{2}{\sqrt{2}}$, $\alpha =\frac {1}{\sqrt {2}}$ ell punto $(u_e,v_e) \approx (0.0131,0.0131)$ foco estable.

Figura 5.3: Para los valores de parámetros $Q=0.0001$, $S=0.5$, $\alpha =0.9999$ el punto $(u_{e},v_{e})\approx (0.9999,0.9999)$ nodo estable.

Conclusiones

En este trabajo se realizó el estudio de la dinámica de un modelo depredador presa del tipo Leslie-Gower considerando competencia (o interferencia) entre los depredadores CED, descrito por el sistema (3.1).

Hemos detallado también las principales propiedades del modelo de Leslie-Gower propuesto en 1948 [28,29] (sistema (2.1), con el objeto de efectuar un estudio comparativo sobre la influencia que tiene en el modelo original, la forma matemática de describir la CED.

Con el fin de simplificar los cálculos, se efectuó una reparametrización y un reescalamiento del tiempo se ha obtenido el sistema (3.2) topológicamente equivalente a (3.1).

Las principales caracterısticas del modelo de tipo Leslie-Gower modificado con CED son:

• Existe una región positivamente invariante.

• Las soluciones son acotadas, es decir el sistema está bien propuesto (well-posed [8]).

• El modelo con CED no está definido en $\left( 0,0\right)$ y tampoco la matriz Jacobiana del sistema topológicamente equivalente (3.1) no está definida para $% y=0$. La naturaleza de los puntos de equilibrio sobre el $eje-x$, se determinaron usando el método del Blowing-up, en el sistema (3.2), cualitativamente equivalente al sistema (3.1).

• El punto de equilibrio $\left( K,0\right)$ es una silla no-hiperbólica.

• Existe una única singularidad positiva (al interior del primer cuadrante), la que es global asintóticamente estable (gae), para todo valor de los parámetros $\left( Q,S,\alpha \right) \in \mathbb{R}% _{+}^{2}\times ]0,1[$.

• No posee ciclos (ni ciclos lımites) al interior del primer cuadrante ( $Int\left( %TCIMACRO{\U{211d} }% \mathbb{R} %EndExpansion _{0}^{+}\right) ^{2}$).

Por lo tanto, ambos modelos, con o sin CED, tienen la misma dinámica, pues en ambos los puntos de equilibrio tienen las mismas caracterısticas (naturaleza), salvo la hiperbolicidad de los puntos de equilibrio $% \left( K,0\right)$ en cada uno de los modelos.

Tanto en el sistema (3.1), como en el modelo de Leslie-Gower sin interferencia entre los depredadores descrito por el sistema (2.1), existe un único punto de coexistencia de las especies, que es atractor global asintóticamente estable (gae). Podemos afirmar entonces que ambos modelos son estructuralmente estables y además no existen bifurcaciones.

Esto significa que cualquiera sea el tamaño inicial de cada de las poblaciones, a medida que transcurra el tiempo, los tamaños poblacionales se estabilizan en dicho punto de equilibrio. Esta propiedad de ambos modelos es muy deseable por algunos ecólogos, o agencias que administran la explotación de recursos renovables.

Por ejemplo, si se considera que un cardumen de peces es capturado por muchos interesados al mismo tiempo (pescadores artesanales o industriales, esto es, los seres humanos actuando como depredadores [12]), la pesquerıa serıa sustentable. Así se deduce en nuestro análisis, asumiendo la CED descrita por la función $B\left( x,y\right) =$ $h\left( x\right) y^{\alpha }$, con $% 0<\alpha <1$.

Sin embargo, no podemos concluir que la CED induce la estabilidad de la interacción depredador-presa, pues sabemos que en muchas de las interaciones del mundo real, existen extinciones de una o ambas poblaciones [12], o bien oscilaciones periódicas de los tamaños poblacionales, las que pueden ser explicados por la existencia de ciclos límites [1,2,13] en el sistema dinámico que describe al modelo.

Estas diferentes situaciones son obtenidas como alternativas dinámicas posibles, en modelos que consideran la CED descrita por otras formas matemáticas. Tales casos se obtienen al asumir una respuesta funcional del tipo Beddington-DeAngelis [21] en el modelo de Leslie-Gower [40].

En este modelo, el punto de equilibrio positivo, no necesariamente es gae, cuando existe, y además existen ciclos lımites alrededor de un punto de equilibrio positivo, lo cual implica oscilaciones de los tamaños poblacionales [40]. Para un mismo conjunto de parámetros, coexisten comportamientos dinámicos distintos, debido a la aparición de curvas separatrices en el plano de fase [40]. Así, las trayectorias o soluciones del sistema son altamente sensibles a las condiciones iniciales.

De este modo surge la disyuntiva, sobre cual es el modelo más eficiente para representar una determinada interacción de depredación en la naturaleza. Nosotros creemos que una misma interacción puede ser representada por un modelo determinado (cualquiera sea su naturaleza matemática) en un momento o tiempo definido, pero su validez dependerá fuertemente de las hipótesis subyacentes, que usualmente no son declaradas.

Teniendo presente que la compledjidad dinámica de un sistema de EDO, no implica la complejidad ecológica del modelo descrito por ese sistema, no podríamos afirmar fehacientemente, que un modelo es mejor que otro. Sin embargo, situaciones ecológicas más complejas (por ejemplo considerar tiempo de gestación, distribución no-homogénea de las poblaciones o clasificados por edades y sexo, etc.) exigen el empleo de herramientas matemáticas más sofisticadas que las EDO.

Referencias Bibliográficas

1

Aguirre P., González-Olivares E., Sáez E. Two limit cycles in a Leslie.Gower predator-prey model with additive Allee effect. Nonlinear Analysis: Real World Applications. 2009; 10:1401-1416.

2

Aguirre P, González-Olivares E, Sáez E. Three limit cycles in a Leslie-Gower predator-prey model with additive Allee effect. SIAM Journal on Applied Mathematics. 2009; 69:1244-1262.

3

Arancibia-Ibarra C., and González-Olivares E. 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 Co.Pte. Ltd., Singapore. 2011; 146-162.

4

Arrowsmith D.K, Place C.M. Dynamical System. Differential equations, maps and chaotic behaviour. Chapman and Hall; 1992.

5

Aziz-Alaoui M.A., Daher Okye M. Boundedness and global stability for a predator-prey model with modified Leslie-Gower and Holling-type II schemes. Applied Mathematics Letters. 2003; 16:1069-1075.

6

Bacaër N. A short history of Mathematical Population Dynamics. Springer-Verlag. 2011.

7

Bazykin A. Nonlinear Dynamics of interacting populations. World Scientific Publishing Co. Pte. Ltd., 1998.

8

Berryman A. A., Gutierrez A.P., Arditi R. Credible, parsimonious and useful predator-prey models. A reply to Abrams, Gleeson and Sarnelle, Ecology. 1995; 76:1980-1985.

9

Birkhoff G.D, Rota G.C. Ordinary Differential Equations (Fourth Edition), John Wiley and Sons, New York; 1989.

10

Chicone C. Ordinary differential equations with applications. 2nd edition, Texts in Applied Mathematics 34, Springer, 2006.

11

Clark C.W. Mathematical Bioeconomic: The optimal management of renewable resources. (2nd edition). John Wiley and Sons, 1990.

12

Clark C.W. The worldwide crisis in fisheries: Economic models and human behavior. Cambridge Univerity Press, 2006.

13

Coleman C.S. Hilbert's 16th. Problem: How many cycles? In: Braun M., Coleman C.S., and Drew D. (Eds.). Differential Equations Model, Springer-Verlag. 1983; 279-297.

14

Dumortier F., Llibre J., Artıs J.C. Qualitative theory of planar differential systems, Springer; 2006.

15

Epstein J.M. Nonlinear Dynamics, Mathematical Biology, and Social Science, Addison-Wesley, 1997.

16

Freedman H.I. Deterministic Mathematical Model in Population Ecology, Marcel Dekker, 1980.

17

Freedman H.I. Stability analysis of a predator-prey system with mutual interference and density-dependent death rates. Bulletin of Mathematical Biology. 1979; 41:67-78.

18

Gallego-Berrío L.M. Consecuencias del efecto Allee en el modelo de depredación de May-Holling-Tanner. Tesis de Maestrıa, Maestrıa en Biomatemáticas Universidad del Quindıo, Armenia, Colombia, 2004.

19

Gallegos-Zuñiga J. Modelo depredador-presa del tipo Leslie-Gower considerando interferencia entre los depredadores. Trabajo para optar al grado de Licenciado en Matemática, Instituto de Matemáticas, Pontificia Universidad Católica de Valparaıso, Chile 2014.

20

Gause G.F. The Struggle for existence. Dover, 1934.

21

Geritz S., Gyllenberg M. A mechanistic derivation of the DeAngelis-Beddington functional response. Journal of Theoretical Biology. 2012; 314:106-108.

22

Goh B-S. Management and Analysis of Biological Populations. Elsevier Scientific Publishing Company, 1980.

23

González-Yañez B., González-Olivares E., Mena-Lorca J. Multistability on a Leslie-Gower Type predator-prey model with nonmonotonic functional response. In R. Mondaini and R. Dilao (eds.), BIOMAT 2006 - International Symposium on Mathematical and Computational Biology, World Scientific Co. Pte. Ltd., 2007; 359-384.

24

González-Olivares E., Mena-lorca J., Rojas-Palma A., Flores J.D. Dynamical complexities in the Leslie-Gower predator-prey model as consequences of the Allee effect on prey. Applied Mathematical Modelling. 2011; 35:366-381.

25

González-Olivares E., Arancibia-Ibarra C., Rojas A, González-Yañez B. Dynamics of a Leslie-Gower predation model considering a generalist predator and the hyperbolic functional response. Mathematical Biosciences and Engineering. 2019; 16(6):7995-8024.

26

Holling C.S. The components of predation as revealed by a study of small-mammal predation of the European pine sawfly. Canadian Entomologist. 1959; 91:293-320.

27

Korobeinikov A.A. Lyapunov function for Leslie-Gower predator-prey models. Applied Mathematical Letters. 2001; 14:697-699.

28

Leslie P.H. Some further notes on the use of matrices in population mathematics. Biometrika. 1948; 35:213-245.

29

Leslie P.H, and Gower J.C. The properties of a stochastic model for the predator-prey type of interaction between two species. Biometrka. 1960; 47:219-234.

30

May R.M. Stability and complexity in model ecosystems (2nd edition), Princeton University Press. 2001.

31

Murray J.D. Mathematical Biology. Springer - Verlag New-York 1989.

32

Perko L. Differential Equations and Dynamical Systems. Springer-Verlag, 1991.

33

Sáez E., González-Olivares E. Dynamics on a predator-prey model. SIAM Journal of Applied Mathematics. 1999; 59:1867-1878.

34

Scudo F.M., Ziegler J.R. The golden age of Theoretical Ecology 1923-1940. Lecture Notes in Biomathematics 22. Springer-Verlag, Berlin 1978.

35

Tanner J.T. The stability and the intrinsic growth rates of prey and predator populations. Ecology. 1975; 56(4):855-867

36

Taylor R.J. Predation. London: Chapman and Hall, 1984.

37

Tintinago-Ruiz P.C., Gallego-Berrío L.M., González-Olivares E. Una clase de modelo de depredación del tipo Leslie-Gower con respuesta funcional racional no monotónica y alimento alternativo para los depredadores. Selecciones Matemáticas. 2019; 06(2):204-216.

38

Turchin P. Complex population dynamics. A theoretical/empirical synthesis. Monographs in Population Biology 35 Princeton University Press 2003.

39

Volterra V. Variazioni e fluttuazioni del numero d'individui in specie animali conviventi. Memorie della R. Accademia dei Lincei, S.VI, IT 1926; II:31-113.

40

Vera-Damián Y., Vidal C., González-Olivares E. Dynamics and bifurcations of a modified Leslie-Gower type model considering a Beddington-DeAngelis functional response. Mathematical Methods in the Applied Sciences. 2019; 42:3179-3210.

Notas al pie

... González-Olivares

Pontificia Universidad Católica de Valparaıso, Chile (rlopezc@ulima.edu.pe).

...ñiga

Pontificia Universidad Católica de Valparaıso, Chile (rlopezc@ulima.edu.pe).