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

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

modelo matemático de una cadena alimenticia depredador-presa: plancton-anchoveta

mathematical model of a predator-prey food chain: plankton-anchovy

Neisser Pino Romero^{∗} Christian Ulises Salazar Fernández^{†}

Received: 21/10/2020; Revised: 19/07/2021; Accepted: 27/09/2021

∗Universidad Peruana Cayetano Heredia, Departmento de Ciencias Exactas, Facultad de Ciencias y Filosofía, Lima, Perú. E-Mail: neisser.pino@upch.pe

†Universidad Nacional de Ingeniería, Departamento académico de métodos cuantitativos, Facultad de Ingeniería Económica, Estadística y Ciencias Sociales, Lima, Perú. E-Mail: csalazarf@uni.edu.pe

70 n. pino – c. salazar

Resumen

En el presente trabajo, se construye un modelo matemático que representa la dinámica que existe entre el fitoplancton (F), el zooplancton (Z) y la anchoveta (A) a partir del modelo de Samares y Anal. Esta cadena alimenticia se presenta en el área marítima peruana donde se encuentra un ecosistema de tres eslabones, además que la cadena alimenticia es un eje principal en el equilibrio ecológico dentro del mar. Se estudiará el comportamiento de las poblaciones (especies marítimas) donde se determinará la estabilidad del modelo (criterio de Routh-Hurwitz) en el largo plazo, donde se analizará con el criterio de Dulac-Bendixon para evaluar la existencia de órbitas periódicas cerradas. Y las respectivas simulaciones computacionales para complementar el estudio e interpretar las situaciones que afectan al ecosistema marítimo.

Palabras clave: ecología matemática; ecuaciones diferenciales ordinarias; cadena alimenticia; estado estacionario; simulación computacional.

Abstract

In the present work, a mathematical model is built that represents the dynamics that exist between phytoplankton (F), zooplankton (Z) and anchovy (A) from the model of Samares and Anal. This food chain occurs in the Peruvian maritime area where there is a three-link ecosystem, in addition to the fact that the food chain is a main axis in the ecological balance within the sea. The behavior of the populations (maritime species) will be studied where the stability of the model (Routh-Hurwitz criterion) in the long term will be determined, where it will be analyzed with the Dulac-Bendixon criterion to evaluate the existence of closed periodic orbits. And the respective computational simulations to complement the study and interpret the situations that affect the maritime ecosystem.

Keywords: mathematical ecology; ordinary differential equations; food chain; stationary state; computational simulation.

Mathematics Subject Classification: 34A34, 34A45, 34C07.

1 Introducción

Se realizará una modificación al modelo desarrollado en el artículo “Dynamics of the interaction of plankton and planktivorous fish with delay” por Samares y Anal [36] como también se utilizará como referencia el artículo de Debasis, Alakes y Guruprasad [22] donde se consideran la dinámica de interacción entre el fitoplancton, el zooplancton y una especie de pez [5, 38], en nuestro caso la anchoveta [16]. La cadena alimenticia considerada en el presente trabajo es de tres eslabones [5, 37], que se representa en la depredación del zooplancton hacia el fitoplancton, y a la vez, el zooplancton es depredado por la anchoveta [13, 17]. Pero una cadena interesante es la que se genera con la anchoveta en cuanto a su alimentación dentro del ecosistema como depredador principal [14, 15]. Con este estudio se puede focalizar en aspectos económicos, ecológicos y geopolíticas de un país en cuanto a su desarrollo industrial [9, 31].

Según el Instituto del Mar del Perú (IMARPE), la anchoveta se alimenta tanto del fitoplancton como el zooplancton [14]; esto ayudará a poder visualizar como este sistema ecológico se mantiene en una coexistencia en el tiempo en un ambiente cerrado sin intervención del hombre [6, 15]. Por otro lado, la anchoveta es una especie fundamental en la fauna marina beneficiosa para el país [3, 23], tanto como alimento de otras especies marítimas que formarían parte de eslabones superiores de la cadena alimenticia [14, 29], como también generadora de actividad pesquera artesanal e industrial. Esta actividad no se incluirá en el análisis del modelo matemático de la cadena alimenticia [2, 5].

La pesca artesanal y empresarial se centran en una actividad económica que se realiza dentro del área marítima peruana [2, 25]. Esta acción pesquera se representa en la mortalidad de la anchoveta que estará considerado dentro del modelo propuesto de manera implícita sin considerar la dinámica que interviene en el ecosistema [3, 31]. La pesquería ayuda a reducir la población de la anchoveta para que no haya una sobrepoblación donde se pueda generar una depredación total del fitoplancton y del zooplancton [5, 29], y la vez es provechosa para el consumo humano y soporte de una beneficiosa actividad económica [16, 38].

Figura 1: Cadena alimenticia entre el fitoplancton, el zooplancton y la anchoveta. Fuente: Página web Significados [37].

Se mostró el diagrama de la cadena alimenticia en la Figura 1, donde de manera directa se genera el equilibrio ecológico en el ecosistema dentro de los eslabones considerados sin una intervención externa que pueda ocasionar su ruptura [13, 17, 26]. Pero la interacción humana se encuentra de manera implícita en el modelo, de manera específica en la tasa de mortalidad de la anchoveta [22, 36] debido que esta actividad tiene una dinámica propia y se debe tener mayores hipótesis para poder expresarlo de manera explícita de como impacta en el ecosistema en el tiempo [3, 6]. La reproducción de las especies marinas tienen un crecimiento rápido, de forma específica para el fitoplancton (crecimiento logístico) como también para el zooplancton que su crecimiento depende de la depredación que realiza al fitoplancton, y esta dinámica se desarrolla dentro del ecosistema considerado como eslabones menores de la anchoveta [14, 23].

La mayor importancia de la reproducción de la anchoveta, es porque representa un soporte alimenticio de otros peces, aves y mamíferos [14] dentro del ecosistema marítimo donde se incorpore una cadena alimenticia con un mayores eslabones [29]; esta consideración no se representará en el modelo matemático propuesto [9, 36]. Además que es una de las especies pelágicas de mayor importancia debido a los grandes volúmenes de captura con destino a una vasta producción empresarial (consumo externo) y artesanal (consumo interno) [2, 15]. La actividad pesquera se centra en periodos normales y es capturada en la franja costera, dentro de las 60 millas náuticas y a profundidades menores de 100 metros como lo expresa el IMARPE [15, 38]. Donde también la eficiencia y equidad en la pesca se debe realizar de una forma controlada y supervisada de tal forma que la captura y esfuezo sea lo más adecuada posible para beneficio de la actividad económica con una vigilancia ecológica [6, 31].

Por lo cual, la coexistencia entre la anchoveta y el plancton (fitoplancton y zooplancton) es fundamental para el equilibrio biológico entre las especies marinas en los eslabones inferiores de un ecosistema más grande del estudiado en este trabajo [27, 36], y de forma directa con la actividad humana debido que genera una alta producción pesquera tanto por la misma anchoveta como también de otras especies que dependen de la anchoveta [6, 14, 31].

Además, el cuidado de la anchoveta de que mantenga sus volumen población beneficia a la producción pesquera del ser humano porque la anchoveta es el alimento de una abundante variedad de especies pertenecientes a la cadena alimenticia siendo un eje importante [3, 9], pero en el presente trabajo no se considerará en la construcción del modelo matemático debido que se quiere cimentar un estudio solido con tres eslabones y la dinámica considerada de depredación para la construcción del modelo [5, 17, 38].

2 Modelo matemático

Sean F, Z, A ∈ C^{1} ([0; +∞)), es decir, F, Z, A : ℝ_{0}^{+} → ℝ_{0}^{+} son diferenciables en todo el dominio ℝ_{0}^{+}, [18]. Sean las variables no negativas (poblaciones ecológicas) del fitoplancton (F), zooplancton (Z) y de la anchoveta (A) y los parámetros (tasas ecológicas), con lo cual se podrá formular el modelo matemático, y además estarán definidas con su respectivo valor positivo que se utilizarán en la Sección 3, cual corresponde a la Simulación Computacional del sistema (1) para generar interpretaciones de sentido ecólogico [4, 21, 33].

El modelo describe la dinámica de la cadena alimenticia que surge entre el fitoplancton con el zooplancton y como la anchoveta se alimentan de estas dos especies siendo un detalle adicional presentado por Samares y Anal [36], y teniendo como referencia sobre la evolución de las poblaciones en el tiempo de acuerdo a los reportes registrados por el IMARPE [22, 14]. El enfoque matemático que se formulará en el modelo intenta describir la dinámica que se presentan entre las especies aunque se tendrá un detalle en cuanto a la evolución de la dinámica de depredación que será constante en todo el tiempo considerado por el modelo que evolucionará a tasas constantes [6, 27, 40].

La reproducción del fitoplancton (F) se modelará con un crecimiento logístico mientras que la reproducción del zooplancton (Z) y la anchoveta (A) se reproducirán de acuerdo a la depredación que realicen [22, 29]. Se utilizará una respuesta funcional Holling tipo II en cada interacción dinámica que se presentan entre las especies [26, 27], porque en su dinámica de depredación se presenta una tasa de ingesta desacelerada, que se deriva de la suposición de que el consumidor está limitado por su capacidad para procesar alimentos [26, 27]. Con este detalle se tendrá una distinción en las tasas de depredación y de beneficio por la depredación que se suscita en el ecosistema [36, 40].

Además de una tasa de mortalidad por causas naturales en cada especie de plancton (fitoplancton y zooplancton); para la especie de la anchoveta se considera una tasa de mortalidad por causas naturales o por el efecto de la pesca realizada tanto por la pesca artesanal o la pesca industrial que se realiza en el mar peruano [3, 9]; este parámetro es importante debido que la pesca puede generar un desequilibrio ecológico [23, 25, 38]. Como no se considera la dinámica que ocasiona la pesca de manera directa y explícita, el sistema tiene una debilidad en el modelaje ecológico presentado pero permite tener un bosquejo de como evoluciona en el corto plazo [6, 31]; en este sentido se debería plantearse hipótesis adicionales que representen de forma más cerca a la problemática de la pesca artesanal e industrial dentro del ecosistema lo cual abriría un estudio más amplio con mejores resultados a mediano plazo [9, 31].

El modelo presentado por Samares y Anal [36] es la base del modelo presentado en (1) con las tasas ecológicas (parámetros) consideradas por los autores. Pero con la introducción de la depredación del fitoplancton por la anchoveta donde se utiliza la respuesta funcional Holling tipo II [36, 40]. El enfoque del modelo intentará describir como la anchoveta es un depredador del plancton (fitoplancton y zooplancton), pero a la vez, la anchoveta es un eslabón de una cadena alimenticia más grande; y es por esto que la anchoveta es importante para otras especies de la cadena alimenticia donde ella se convierte en presa [3, 31]; este detalle no se considera en el análisis del presente trabajo [22, 15].

Con las hipótesis mencionadas, se puede formular el siguiente modelo matemático que tiene base del modelo presentado en [36]. Donde se considera la dinámica de interacción en el ecosistema marítimo de nuestras poblaciones (fitoplancton, zooplancton y anchoveta) en un sistema no cerrado [25, 40]:

\{ \begin{aligned}\frac{dF}{dt} & = & rF (1 − \frac{F}{K} ) − \frac{α_{1}F ·Z}{k_{1}+F} − \frac{α_{2}F ·A}{k_{2}+F} − μ_{1}F & , & F (0) = F_{o} ≥ 0 \\ \frac{dZ}{dt} & = & \frac{β_{1}F ·Z}{k_{1}+F} − \frac{α_{3}Z·A}{k_{3}+Z} − μ_{2}Z & , & Z(0) = Z_{o} ≥ 0 \\ \frac{dA}{dt} & = & \frac{β_{2}F ·A}{k_{2}+F} + \frac{β_{3}Z·A}{k_{3}+Z} − μ_{3}A & , & A(0) = A_{o} ≥ 0 \end{aligned} ; r, K, k_{1}, k_{2}, k_{3}, α_{1}, α_{2}, α_{3}, β_{1}, β_{2}, β_{3}, μ_{1}, μ_{2}, μ_{3} > 0 ; α_{1} > β_{1} , α_{2} > β_{2} , α_{3} > β_{3}

(1)

De la ecuación (1) donde representa el modelo matemático que se estudiará para analizar la dinámica entre el depredador y la presa de nuestra cadena alimenticia de tres eslabones del ecosistema marítimo [15, 27, 36].

2.1 Existencia y positividad de soluciones

Primero, se mostrará la existencia de las soluciones y luego la positividad de las mismas del sistema (1) que modela especies ecológicas, y no pueden ser negativas [27]. Con estas justificaciones, se formalizará las hipótesis del fenómeno ecológico entre las tres especies que habitan en el mar (cadena alimenticia de tres eslabones) para la buena definición del modelo matemático. Por otro lado, el modelo propuesto expresado en (1), es no homogéneo [8, 12, 21].

Sean las funciones asociadas al modelo (1) de la siguiente forma:

\begin{aligned}\frac{dF}{dt} = G_{1}(F, Z, A) & = & rF (1 − \frac{F}{K} ) − \frac{α_{1}F ·Z}{k_{1}+F} − \frac{α_{2}F ·A}{k_{2}+F} − μ_{1}F \\ \frac{dZ}{dt} = G_{2}(F, Z, A) & = & \frac{β_{1}F ·Z}{k_{1}+F} − \frac{α_{3}Z·A}{k_{3}+Z} − μ_{2}Z \\ \frac{dA}{dt} = G_{3}(F, Z, A) & = & \frac{β_{2}F ·A}{k_{2}+F} + \frac{β_{3}Z·A}{k_{3}+Z} − μ_{3}A \end{aligned}  (2)

De la ecuación (2) se determinan sus derivadas parciales como sigue:

• \frac{dG_{1}}{dF} = r−\frac{2rF}{K} −\frac{α_{1}k_{1}Z}{(k_{1}+F)^{2}} −\frac{α_{2}k_{2}A}{(k_{2}+F)^{2}} −μ_{1} , \frac{dG_{1}}{dZ} = −\frac{α_{1}F}{k_{1}+F} , \frac{dG_{1}}{dA} = −\frac{α_{2}F}{k_{2}+F}

• \frac{dG_{2}}{dF} = \frac{β_{1}k_{1}Z}{(k_{1}+F)^{2}} , \frac{dG_{2}}{dZ} = \frac{β_{1}F}{k_{1}+F} − \frac{α_{3}k_{3}A}{(k_{3}+Z)^{2}} − μ_{2} , \frac{dG_{2}}{dA} = −\frac{α_{3}Z}{k_{3}+Z}

• \frac{dG_{3}}{dF} = \frac{β_{2}k_{2}A}{(k_{2}+F)^{2}} , \frac{dG_{3}}{dZ} = \frac{β_{3}k_{3}A}{(k_{3}+Z)^{2}} , \frac{dG_{3}}{dA} = \frac{β_{2}F}{k_{2}+F} + \frac{β_{3}Z}{k_{3}+Z} − μ_{3}

Con estas derivadas parciales asociadas a la ecuación (2), se podrá satisfacer las condiciones del Teorema 1 [21, 33].

Teorema 1 Sea la función G : ℝ_{+}^{3} → ℝ_{+}^{3}, donde G(x) = (G_{1}(x), G_{2}(x), G_{3}(x)) con x = (F, Z, A) son funciones continuas y existen \frac{∂G_{j}}{∂x_{j}} continuas en ℝ_{+}^{3} ∀j = \overline{1; 3}. Entonces, la función G es localmente Lipschitz continua en ℝ_{+}^{3}.

A partir del Teorema 1, se sigue de la siguiente forma para que el modelo (1) satisfaga una hipótesis del Teorema 2, y así garantizar la existencia y positividad de soluciones del sistema [18, 21].

• G_{1}(0 ; Z > 0 ; A > 0) = 0 ≥ 0
• G_{2}(F > 0 ; 0 ; A > 0) = 0 ≥ 0
• G_{3}(F > 0 ; Z > 0 ; 0) = 0 ≥ 0

Teorema 2 Sea G : ℝ_{+}^{3} → ℝ_{+}^{3} localmente lipschitz continua en ∀j = \overline{1; 3} y satisface G_{j}(x) ≥ 0 para cualquier x ∈ ℝ_{+}^{3}, x_{j} = 0. Entonces para cada x_{0} ∈ ℝ_{+}^{3}, existe una única solución de x′ = G(x) con x(0) = x_{0} en ℝ_{+}^{3} donde está definida en algún intervalo (0; b] con b ∈ (0; +∞].

2.2 Solución acotada del modelo

Ahora, se formulará el Teorema 3 que nos ayudará a mostrar que el modelo (1) alcanza una solución acotada en el largo plazo. Con este teorema se complementará a los dos teoremas anteriores, Teorema 1 y Teorema 2, donde se mostró la existencia, la unicidad y la positividad de las soluciones del modelo [8, 36].

Teorema 3 Todas las soluciones del Sistema (1) están acotadas en el largo plazo.

Demostración. Definimos W (t) = F (t) + Z(t) + A(t) ≥ 0. De la primera ecuación del sistema (1), se tiene

\frac{dF}{dt} = rF (1 − \frac{F}{K} ) − \frac{α_{1}F · Z}{k_{1} + F} − \frac{α_{2}F · A}{k_{2} + F} − μ_{1}F ≤ rF (1 − \frac{F}{K} ) .

Las soluciones positivas con las condiciones iniciales positivas están definidas debajo de la solución de la siguiente ecuación:

\frac{dF}{dt} = rF (1 − \frac{F}{K} ) .  (3)

Con lo cual, de la ecuación (3), se puede obtener su solución analítica en la condición inicial de (1). Considerando una constante de integración (C) en la solución exacta, [35].

F (t) = \frac{K}{1 + Ce^{−(rK)t}} .  (4)

Luego, todas las soluciones de F (t) tienden a K cuando t ⟶ +∞. Es decir,

lim_{t⟶+∞} sup F (t) = K  (5)

De la función W (t), se puede expresar W′(t) = F′(t) + Z′(t) + A′(t). Por lo cual, se obtiene lo siguiente

\begin{aligned}\frac{dW}{dt} & = & rF (1 − \frac{F}{K} ) + \frac{(β_{1}−α_{1})F Z}{k_{1}+F} + \frac{(β_{2}−α_{2})F A}{k_{2}+F}+ \\ & & \frac{(β_{3}−α_{3})ZA}{k_{3}+Z} − μ_{1}F − μ_{2}Z − μ_{3}A \end{aligned} .  (6)

Se definirá lo siguiente, μ_{o} = min{μ_{1}, μ_{2}, μ_{3}}. Además, el término de interacción entre presa y depredador (β_{j} − α_{j} , ∀j = \overline{1; 3}) es negativa porque la tasa de depredación es mayor a la tasa de beneficio por la depredación [26, 27]. Luego, se obtiene la siguiente ecuación

W′(t) ≤ rF − μ_{o}W (t).  (7)

De una forma similar que se realizó para la primera ecuación del sistema (1) para obtener la solución acotada, se construirá la ecuación (8). Además, se tiene la condición inicial positiva W (0) = F (0) + Z(0) + A(0) ≥ 0 [36]:

W′(t) = rF − μ_{o}W (t) , W (0) ≥ 0  (8)

Con lo cual, es fácil obtener la solución exacta a la ecuación (8):

W (t) = \frac{rK}{μ_{o}} + (W (0) − \frac{rK}{μ_{o}} ) · e^{−μ_{o}t}.  (9)

Así, todas las soluciones de W (t) tienden a \frac{rK}{μ_{o}} cuando t ⟶ +∞. Es decir,

lim_{t⟶+∞} sup W (t) = \frac{rK}{μ_{o}} .  (10)

Por consiguiente, se obtiene el siguiente resultado

lim_{t⟶+∞} sup\{F (t) + Z(t) + A(t)\} ≤ \frac{rK}{μ_{o}} .  (11)

El resultado obtenido en la ecuación (11) nos indica que la solución del sistema (1) está acotada en el largo plazo. Es decir, la sumatoria de todas las poblaciones en el largo plazo no llegan a sobrepasar la cota \frac{rK}{μ_{o}}. Por consiguiente, se obtiene la región factible del sistema [27, 36]:

ℭ = \{(F, Z, A) ∈ ℝ_{+}^{3} : 0 ≤ F ; 0 ≤ Z ; 0 ≤ A; F + Z + A ≤ \frac{rK}{μ_{o}} \}. (12)

Finalmente, se ha deteminado que nuestro sistema (1) tiene solución única y es positiva (Teorema 2), además que está acotada en el largo plazo (Teorema 3) satisfaciendo las condiciones para un modelo matemático aplicado a la ecología.

2.3 Análisis cualitativo

Después de haber presentado la existencia y unicidad, y además de positividad de soluciones, y las soluciones del sistema (1) están acotadas en el largo plazo. Se procederá a realizar el análisis cualitativo del modelo (1) que nos permite saber como es su comportamiento a un largo plazo de tiempo. Por lo cual se determinará los puntos críticos o puntos de equilibrio en el estado estacionario del modelo matemático [7, 33, 36].

2.3.1 Puntos de equilibrio

Los puntos de equilibrio (puntos críticos) lo determinaremos igualando a cero cada ecuación del sistema de ecuaciones (1) [21, 27]:

\{ \begin{matrix}rF (1 − \frac{F}{K} ) − \frac{α_{1}F ·Z}{k_{1}+F} − \frac{α_{2}F ·A}{k_{2}+F} − μ_{1}F & = & 0 \\ \frac{β_{1}F ·Z}{k_{1}+F} − \frac{α_{3}Z·A}{k_{3}+Z} − μ_{2}Z & = & 0 \\ \frac{β_{2}F ·A}{k_{2}+F} + \frac{β_{3}Z·A}{k_{3}+Z} − μ_{3}A & = & 0. \end{matrix}  (13)

Al resolver el sistema (13), se obtienen cuatro puntos críticos explícitos asociados al sistema (1). Además se considerará la existencia de un punto de coexistencia (e_{5}), siendo el estado estacionario buscado por la ecología porque es ahí donde se presenta el equilibrio biológico en el ecosistema. Aunque se tenga en consideración de forma teórica de este estado estacionario [26, 27, 36]:

\begin{aligned}e_{1} & = & ( 0 ; 0 ; 0 ) \\ e_{2} & = & (\frac{K(r−μ_{1})}{r} ; 0 ; 0 ) \\ e_{3} & = & (\frac{k_{2}μ_{3}}{β_{2}−μ_{3}} ; 0 ; \frac{k_{2}β_{2}(K(r−μ_{1})(β_{2}−μ_{3})−k_{2}rμ_{3})}{Kα_{2}(β_{2}−μ_{3})^{2}} ) \\ e_{4} & = & (\frac{k_{1}μ_{2}}{β_{1}−μ_{2}} ; \frac{k_{1}β_{1}(K(r−μ_{1})(β_{1}−μ_{2})−k_{1}rμ_{3})}{Kα_{1}(β_{1}−μ_{2})^{2}} ; 0 ) \\ e_{5} & = & ( F^{∗} ; Z^{∗} ; A^{∗} ) . \end{aligned}  (14)

De la ecuación (14), se tienen los puntos de equilibrio que permiten analizar la estabilidad local asociada al sistema (1) en ese estado estacionario. Los puntos críticos deben ser no negativos debido que las funciones representan poblaciones de las especies ecológicas, como se mostró en el Teorema 2 [8, 32, 36]:

• Extinción total (e_{1}) : Existencia trivial.

• Supervivencia del fitoplancton (e_{2}) : Existe si, r > μ_{1}.

• Extinción del zooplancton (e_{3}) :

Existe si r > μ_{1}, β_{2} > μ_{3}, K(r − μ_{1})(β_{2} − μ_{3}) > k_{2}rμ_{3}.

• Extinción de la anchoveta (e_{4}) :

Existe si, r > μ_{1}, β_{1} > μ_{2}, K(r − μ_{1})(β_{1} > μ_{2}) − k_{1}rμ_{3}.

• Coexistencia total (e_{5}) :

Existe si, F^{∗} > 0 , Z^{∗} > 0 , A^{∗} > 0.

Una obsevación importante sobre el punto de coexistencia (e_{5}) de nuestro sistema (1) es la hipótesis de su existencia como una expresión natural que se presenta en el ecosistema por el equilibrio biológico [27, 36]. A pesar que no se ha podido determinar el punto de manera explícita el punto de equilibro [17, 26].

2.3.2 Estabilidad del modelo

Después de haber obtenido los puntos críticos (14) asociados al sistema. Ahora, se procederá a determinar la matriz Jacobiana (J(F, Z, A)) asociada al modelo matemático (1). Esta matriz representa un sistema lineal asociado al modelo no lineal en un punto de equilibrio donde se encuentra el estado estacionario asociado al sistema no lineal [8, 33].

Utilizaremos las siguientes variables auxiliares para poder expresar mejor la matriz jacobiana asociada al sistema (1) de una forma más accesible [21]:

• q_{1} = r − \frac{2rF}{K} − \frac{α_{1}k_{1}Z}{(k_{1}+F)^{2}} − \frac{α_{2}k_{2}A}{(k_{2}+F)^{2}} − μ_{1}.

• q_{2} = \frac{β_{1}F}{k_{1}+F} − \frac{α_{3}k_{3}A}{(k_{3}+Z)^{2}} − μ_{2}.

• q_{3} = \frac{β_{2}F}{k_{2}+F} + \frac{β_{3}Z}{k_{3}+Z} − μ_{3}.

• n_{1} = \frac{α_{1}F}{k_{1}+F} > 0  ;  n_{2} = \frac{α_{2}F}{k_{2}+F} > 0.
• n_{3} = \frac{β_{1}k_{1}Z}{(k_{1}+F)^{2}} > 0  ;  n_{4} = \frac{α_{3}Z}{k_{3}+Z} > 0.
• n_{5} = \frac{β_{2}k_{2}A}{(k_{1}+F)^{2}} > 0  ;  n_{6} = \frac{β_{3}k_{3}A}{(k_{3}+Z)^{2}} > 0.

La matriz jacobiana asociada al sistema (1) sería el siguiente:

J(F, Z, A) = [ \begin{matrix}q_{1} & −n_{1} & −n_{2} \\ n_{3} & q_{2} & −n_{4} \\ n_{5} & n_{6} & q_{3} \end{matrix} ] .  (15)

De la matriz jacobiana (15) se puede obtener el polinomio característico asociado al sistema (1). Se recurrirá variables auxiliares para expresar una forma más sencilla [32]:

• Q_{1} = −(q_{1} + q_{2} + q_{3}).
• Q_{2} = q_{1}q_{2} + q_{1}q_{3} + q_{2}q_{3} + n_{1}n_{3} + n_{2}n_{5} + n_{4}n_{6}.
• Q_{3} = n_{2}n_{3}n_{6} − (q_{1}q_{2}q_{3} + q_{1}n_{4}n_{6} + q_{3}n_{1}n_{3} + q_{2}n_{2}n_{5} + n_{1}n_{4}n_{5}).

Con estas variables, el polinomio característico queda de la siguiente forma

p(λ) = λ^{3} + Q_{1}λ^{2} + Q_{2}λ + Q_{3}.  (16)

Para determinar la estabilidad local del sistema (1) se utilizará el criterio de Routh-Hurwitz, donde se debe garantizar que cada constante del polinomio característico (16) debe ser positivo:

Q_{1} > 0 ⟺ q_{1} + q_{2} + q_{3} < 0

Q_{2} > 0 ⟺ q_{1}q_{2} + q_{1}q_{3} + q_{2}q_{3} > 0

Q_{3} > 0 ⟺ n_{2}n_{3}n_{6} > q_{1}q_{2}q_{3} + q_{1}n_{4}n_{6} + q_{3}n_{1}n_{3} + q_{2}n_{2}n_{5} + n_{1}n_{4}n_{5}.

(17)

Después de garantizar la primera condición de Routh-Hurwitz expresado en (17), se debe garantizar los menores principales positivos [1, 32]:

\begin{matrix}△_{1} = Q_{1} > 0 (trivial) \\ △_{2} = | \begin{matrix}Q_{1} & 1 \\ Q_{3} & Q_{2} \end{matrix} | = Q_{1}Q_{2} − Q_{3} > 0. \end{matrix}  (18)

Por consiguiente, se satisfacen las condiciones del criterio de Routh-Hurwitz expresadas en las ecuaciones en (17) y (18) para obtener la estabilidad local en el punto de equilibrio asociado al sistema (1) [1, 33].

Teorema 4 Sea el modelo matemático expresado en (1), y posee su matriz jacobiana asociada (15), y además su polinomio característico (16) en el punto de equilibrio (e_{j}) para todo j = \overline{1, 5}. Si Q_{1}, Q_{2}, Q_{3} > 0, además Q_{1}Q_{2} − Q_{3} > 0, entonces el sistema (1) es localmente estable.

Una observación delicada para determinar la estabilidad de cada estado estacionario asociado al punto de equilibrio, son las condiciones de las variables auxiliares [5, 27]. El modelo se centra en la existencia del punto de coexistencia (e_{5}) debido que se busca este estado estacionario para beneficio de las tres especies y la intervención implícita del hombre como también se podría considerar otras especies depredadoras de la anchoveta [2, 25, 38]. Es verdad, que la posibilidad de la coexistencia en un ecosistema abierto como se presenta en el mar peruano viene a ser algo bastante difícil de modelar, por consiguiente, esta problemática no se considera en nuestro análisis, pero quedaría para un futuro trabajo donde se incorpore estas hipótesis para construir un modelo más amplio que ayude a la preservación del equilibrio biológico [15, 27, 31].

Con respecto al punto de equilibrio de extinción total (e_{1}), tiene que ser obligatoriamente inestable debido que no se debe llegar a este escenario trágico para el ecosistema [5, 14, 17], además que este punto es inestable de acuerdo al Teorema 4 [27, 41]. Además, la extinción de la anchoveta que se presenta en el punto de equilibrio (e_{2}), lo cual no ayudaría al equilibrio ecológico del ecosistema marino, y a la producción en la actividad pesquera del país en el tiempo. De forma adicional, este punto de equilibrio es inestable de acuerdo al Teorema 4 que garantiza que no llega a ese estado [23, 27].

Por otro lado, el punto de equilibrio (e_{3}) es un punto donde se extingue el zooplancton pero aún se mantiene la supervivencia de la anchoveta mediante el consumo del fitoplancton [22, 9]; y con respecto al punto (e_{4}) es un estado donde se extingue la anchoveta pero se mantienen las otras dos especies que afecta de manera directa a la producción pesquera y al equilibrio ecológico (depredadores de la anchoveta) de eslabones superiores que no se consideran en el modelo propuesto [27, 29]. Finalmente, el punto de coexistencia (e_{5}) es el punto que siempre se debe buscar y conocer la sensibilidad de los parámetros donde se mantenga en el largo plazo [32, 41], a pesar que sea una situación que es difícil de conseguir por factores externos que en muchos casos la intervención del ser humano causa el desequilibrio ecológico entres las especies marítimas [14, 38].

2.4 Criterio Dulac-Bendixon

Se sabe que un campo X al multiplicarse por una función positiva f, el campo resultante fX tienen el mismo diagrama fase que X. Así, por ejemplo, si uno de ellos tiene órbitas periódicas, el otro también. Como se sabe, la divergencia de un campo X = (f, g) se define como div(X(x, y)) = \frac{∂f}{∂x}(x, y) + \frac{∂g}{∂x}(x, y). Si X tiene una órbita periódica γ y si div(X) ≢ 0 en la región acotada por γ, entonces esta región intersecta tanto div(X) > 0 como div(X) < 0 [4, 24, 42].

A veces ocurre que div(X) cambia de signo en una región en U, pero existe una función f, positiva en U tal que div(fX) no cambia de signo, implicando que X no tiene órbitas periódicas contenidas en U. A esta función f se le llama función de Dulac de X en U [4, 19, 42].

Para poder utilizar el criterio de Dulac-Bendixon, se tendrá que reducir el modelo (1) de tres poblaciones a dos poblaciones (sistema bidimensional) donde se obtenga el plano fase respectivo [30, 41]. Se ha obtenido que las soluciones están acotadas en el largo plazo (Teorema 3). En este sentido, la región factible expresada en (2.4) donde la población total está acotada por \frac{rK}{μ_{o}}. En el largo plazo, el sistema (1)está sujeto a la restricción F +Z+A = \frac{rK}{μ_{o}}, y utilizando esta consideración se tiene A = \frac{rK}{μ_{o}} − F − Z en el modelo. Con esta transformación, se puede eliminar la población de las anchovetas (A) de las ecuaciones [27, 30].

Esta sustitución nos brinda un modelo más simple (fitoplancton-zooplancton):

\{ \begin{matrix}\frac{dF}{dt} & = & rF (1 − \frac{F}{K} ) − \frac{α_{1}F ·Z}{k_{1}+F} − ^{\frac{α_{2}F ( \frac{rK}{μo} −F −Z)}{k_{2}+F}} − μ_{1}F. \\ \frac{dZ}{dt} & = & \frac{β_{1}F ·Z}{k_{1}+F} − ^{\frac{α_{3}Z( \frac{rK}{μo} −F −Z)}{k_{3}+Z}} − μ_{2}Z. \end{matrix}  (19)

El plano fase F − Z asociado al sistema (19) está dado por

ℭ_{F Z} = \{(F, Z) ∈ ℝ_{+}^{2} : 0 ≤ F ; 0 ≤ Z ; F + Z ≤ \frac{rK}{μ_{o}} \} .

Proposición 1 Si \frac{α_{1}K}{μ_{o}(k_{1}+F)^{2}} < \frac{1}{KZ} , 1 + k_{2} < Z y \frac{rK}{μ_{o}} < k_{3} + F, entonces el multiplicador

Φ(F, Z) = \frac{1}{F Z}

es una función de Dulac del sistema (19) en el interior de ℭ_{F Z}.

Demostración. Sea la función Φ = \frac{1}{F Z} para todo F, Z > 0.

• U_{1} = ΦF′ = \frac{r}{Z} − \frac{rF}{KZ} − \frac{α_{1}rK}{μ_{o}(k_{1}+F)} − \frac{α_{2}}{Z(k_{2}+F)} + \frac{α_{2}F}{Z(k_{2}+F)} + \frac{α_{2}}{k_{2}+F} − \frac{μ_{1}}{Z}

• U_{2} = ΦZ′ = \frac{β_{1}}{k_{1}+F} − \frac{α_{3}rK}{μ_{o}F (k_{3}+Z)} + \frac{α_{3}}{k_{3}+Z} + \frac{α_{3}Z}{k_{3}+Z} − \frac{μ_{2}}{F}

Luego, realizamos las derivadas parciales para después obtener la divergencia asociada al sistema.

• \frac{∂U_{1}}{∂F} = −\frac{r}{KZ} − \frac{α_{1}rK}{μ_{o}(k_{1}+F)^{2}} + \frac{α_{2}}{Z(k_{2}+F)^{2}} + \frac{α_{2}k_{2}}{Z(k_{2}+F)^{2}} − \frac{α_{2}}{(k_{2}+F)^{2}}

• \frac{∂U_{2}}{∂Z} = \frac{α_{3}rK}{μ_{o}F (k_{3}+Z)^{2}} − \frac{α_{3}}{(k_{3}+Z)^{2}} − \frac{α_{3}k_{3}}{F (k_{3}+Z)^{2}}

Luego, determinamos la divergencia con el multiplicador Φ(F, Z), y con las hipótesis consideradas en la Proposición 1 se tiene lo siguiente:

div(U_{1}, U_{2}) = [−\frac{r}{KZ} + \frac{α_{1}rK}{μ_{o}(k_{1} + F)^{2}} ] + \frac{α_{2} [−Z + k_{2} + 1]}{Z(k_{2} + F)^{2}}

+\frac{α_{3} [−F − k_{3} + \frac{rK}{μ_{o}} ]}{Z(k_{3} + Z)^{2}} < 0.

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

Por consiguiente, el modelo expresado en (19) admite una función de Dulac (Φ = \frac{1}{F Z}) para todo F, Z > 0. Por lo tanto, el subsistema no posee órbitas periódicas en el interior de la región ℭ_{F Z} [24, 42].

Este plano fase es interesante debido que son las poblaciones del primer y segundo eslabón de la cadena alimenticia considerada en el modelo [26]. Además que estas poblaciones son presas de la anchoveta (tercer eslabón), y nos permite tener un bosquejo de cómo podría afectar el ecosistema marítimo. Asimismo, de la actividad pesquera vigilada y monitoreada por el Instituto del Mar del Perú (IMARPE) [14, 27].

De forma similar que se construyó las hipótesis necesarias para formular la Proposición 1. Se puede introducir la consideración de F = \frac{rK}{μ_{o}} − Z − A en el modelo. Con lo cual, se puede eliminar la población del fitoplancton (F) de las ecuaciones [27, 30]. Esta sustitución nos brinda un modelo más simple (zooplancton-anchoveta) para analizar si el modelo posee órbitas periódicas:

\{ \begin{matrix}\frac{dZ}{dt} & = & \frac{β_{1}F ·Z}{k_{1}+F} − ^{\frac{α_{3}Z( \frac{rK}{μo} −F −Z)}{k_{3}+Z}} − μ_{2}Z. \\ \frac{dA}{dt} & = & \frac{β_{1}F ·Z}{k_{1}+F} − ^{\frac{α_{3}Z( \frac{rK}{μo} −F −Z)}{k_{3}+Z}} − μ_{3}A. \end{matrix}  (20)

El plano fase Z − A asociado al sistema (20) está dado por

ℭ_{ZA} = \{(Z, A) ∈ ℝ_{+}^{2} : 0 ≤ Z ; 0 ≤ A; Z + A ≤ \frac{rK}{μ_{o}} \} .

Proposición 2 Si \frac{α_{3}}{(k_{3}+Z)^{2}} < _{\frac{β_{1}k_{1}}{A(k_{1}+ \frac{rK}{μo} −Z−A)^{2}}} + _{\frac{β_{2}k_{2}}{Z(k_{2}+ \frac{rK}{μo} −Z−A)^{2}}} , entonces el

multiplicador

Φ(Z, A) = \frac{1}{ZA}

es una función de Dulac del sistema (20) en el interior de ℭ_{ZA}.

Demostración. Sea la función Φ = \frac{1}{ZA} para todo Z, A > 0.

• R_{1} = ΦZ′ = \frac{β_{1}( \frac{rK}{μo} −Z−A)}{A(k_{1}+ \frac{rK}{μo} −Z−A)} − \frac{α_{3}}{k_{3}+Z} − \frac{μ_{2}}{A}.

• R_{2} = ΦA′ = \frac{β_{2}( \frac{rK}{μo} −Z−A)}{Z(k_{2}+ \frac{rK}{μo} −Z−A)} + \frac{β_{3}}{k_{3}+Z} − \frac{μ_{3}}{Z}.

Luego, realizamos las derivadas parciales para después obtener la divergencia asociada al sistema:

• \frac{∂R_{1}}{∂Z} = −_{\frac{β_{1}k_{1}}{A(k_{1}+ \frac{rK}{μo} −Z−A)}} + \frac{α_{3}}{(k_{3}+Z)^{2}} .

• \frac{∂R_{2}}{∂A} = −_{\frac{β_{2}k_{2}}{Z(k_{2}+ \frac{rK}{μo} −Z−A)}}.

Luego, determinamos la divergencia con el multiplicador Φ(Z, A), y con las hipótesis consideradas en la Proposición 2 se tiene lo siguiente:

div(R_{1}, R_{2}) = \frac{α_{3}}{(k_{3} + Z)^{2}} − [\frac{β_{1}k_{1}}{A (k_{1} + \frac{rK}{μ_{o}} − Z − A)}+

+\frac{β_{2}k_{2}}{Z (k_{2} + \frac{rK}{μ_{o}} − Z − A)}] < 0.

Por consiguiente, el modelo expresado en (20) admite una función de Dulac (Φ = \frac{1}{ZA}) para todo Z, A > 0. Por lo tanto, el subsistema no posee órbitas periódicas en el interior de la región ℭ_{ZA} [24, 42].

Este plano fase es importante debido que son las poblaciones del segundo y tercer eslabón de la cadena alimenticia considerada en el modelo [15, 27]. Además que la anchoveta tiene como alimento al zooplancton en su eslabón próximo y es predominante conocer su comportamiento en el tiempo porque es una población que es presa de eslabones superiores (especies marítimas) pero también de la población humana (actividad pesquera) [14, 16].

Una consideración de la investigación presentado por Samares y Anal [36] guiado del trabajo Mukhopadhyay y Bhattacharyya [28] permite considerar la existencia de la bifurcación de Hopf en la coexistencia de las especies marítimas. Por otro lado, la consideración del criterio de Dulac-Bendixon nos permite conocer que no poseen órbitas cerradas para los subsistemas (19) y (20). En el presente trabajo no se ha realizado el estudio de la existencia de la bifurcación de Hopf para el estado de coexistencia de las tres especies marítimas [16, 27], pero esta consideración sería importante analizar para bosquejar con qué parámetros se podría obtenerlaen el tiempo [29, 41].

3 Simulación numérica

Dentro del modelamiento matemático aplicado a la Ecología, se ha realizado el Análisis Cualitativo (sección 2) que se relacionará con el Análisis Numérico (soluciones aproximadas) [5, 21]. Esto nos ayuda a visualizar que el análisis matemático con las simulaciones computacionales se complementan mutuamente para mostrar la estabilidad local del modelo matemático (1), donde se puede considerar interpretaciones ecológicas en un corto, mediano y largo plazo bajo las hipótesis consideradas [22, 13, 26]. Para lo cual, se presentará una tabla con los valores a considerar para la respectiva simulación computacional y poder obtener curvas poblacionales en el tiempo [14, 36, 29]. Donde las tasas ecológicas del modelo (1) tendrán su unidad dimensional respectiva como las poblaciones ecológicas (toneladas) [27, 40].

En el Instituto del Mar del Perú (IMARPE) se puede visualizar diversos índices y valores referenciales que representan la dinámica de este ecosistema marítimo [16, 29], que nos permite una posible proyección de la evolución de cada población de acuerdo a la recolección de datos diarios, mensuales, semestrales y anuales [2, 15, 38]. Además que hay un parámetro sensible que es la mortalidad de la anchoveta debido a la acción depredadora de otras especies marítimas dentro del mar (eslabones superiores de la cadena alimenticia) como también la actividad pesquera que se realiza a nivel industrial y artesanal [6, 9]; la pesquería sin medida ni un control guiado y supervisado puede ocasionar un desequilibrio ecológico o hasta la posible extinción de la especie de la anchoveta [5, 23, 25] y con esta trágica situación conllevaría un problema crítico para otras especies como para la actividad pesquera [3, 22, 16]. Estas consideraciones no están en el modelo propuesto debido que tiene una dinámica que debe ser estudiada de forma más completa; dicho estudio no se realizará en este trabajo de investigación pero se intenta considerar la existencia de estos supuestos para un trabajo futuro que permita tener esta dinámica dentro del modelo [3, 31, 41].

Las tasas ecológicas expresadas en la Tabla 1, se han recopilado del artículo de Samares y Anal [36] en cuanto a los parámetros ecológicos de la respuesta funcional Holling tipo II; esto se puede utilizar debido que son datos de la dinámica entre las especies (anchoveta y plancton). De forma particular, se han extraído de los datos provistos por el IMARPE en cuanto a la reproducción, mortalidad y depredación que se estudia en el mar peruano, y se presentan en los reportes sobre las actividades marítimas registradas [14, 16]; se ha omitido la consideración del impacto realizado por la pesquería como dinámica propia [6, 9]. Siempre se ha mantenido una inspección diaria en temporadas críticas causados por factores externos que puedan afectar el ecosistema, esto quiere decir una consideración detallada de las actividades pesqueras [14, 15]. Esto se debe que el ecosistema marítimo tiene una importancia para el país por intereses ecológicos, económicos y geopolíticos con los demás países vecinos [6, 31, 29].

Como se había mencionado, la pesquería es una actividad económica importante para el país tanto por la elaboración de productos derivados de la anchoveta como un eslabón de una cadena mucho mayor según lo manifiesta el Ministerio de Producción del Perú [25], y la Sociedad Nacional de Pesquería [38]. Y por este motivo, se debe tener una vigilancia predominante en la anchoveta por ser el alimento principal de otras especies marítimas que conforman una cadena alimenticia de mayor eslabón que forma parte también de la vigilancia que realiza el IMARPE para mantener el equilibrio biológico en el mar peruano [14, 16], y la interacción de la acción pesquera sin control y supervisión [3, 2, 23].

Esto nos permitiría poder bosquejar la dinámica entre las tres especies al mediano y largo plazo bajo la perspectiva que no hay cambios en las tasas ecológicas (tasas constantes), pero en el corto plazo ya brinda como está evolucionando las especies [17, 29]; también se podría considerar cuando cambian en el tiempo en otro trabajo de investigación [27, 36]. De manera adicional, el trabajo de Medina et al. (2015) nos presenta indicadores de como se alimenta la anchoveta en el norte de Chile que forma parte de la frontera con el sur de Perú, donde se podrá proyectar el comportamiento de la depredación [16] y cómo se puede generar un sistema más amplio donde se tenga una mayor cadena alimenticia que representaría de como se relacionan las especies marítimas en un amplio espectro biológico [14, 23]. Estas perspectivas permiten tener un estudio complementario entre dos naciones que comparten un recurso marítimo importante de actividad económica y manteniendo en mayor medida el equilibrio biológico para salvaguardar todas las especies que comprenden el ecosistema [15, 29, 38].

Tabla 1: Tabla de valores numéricos de las poblaciones iniciales.

Especie  Descripción  Valor Dim.
F  Población inicial del fitoplancton 1.30 T n
Z  Población inicial del zooplancton 0.70 T n
A  Población inicial del anchovetas 0.05 T n

Las unidades de la poblaciones ecológicas (fitoplancton, zooplancton y anchoveta) están medidas en toneladas (Tn), por mientras que las tasas ecológicas están consideradas por unidad de tiempo (días) de acuerdo a la vigilancia que realiza el IMARPE como parte de su labor institucional para servicio de la sociedad de acuerdo a los intereses del Estado [14, 16]. También se considera valores referenciales a la literatura citada [3, 28, 36].

Tabla 2: Tabla de valores numéricos de los parámetros (constantes).

Param.  Descripción  Valor Dim.
r  Tasa de reproducción intrínseca del fitoplancton 1.20 1/T
k Capacidad de carga del fitoplancton 0.60 T n

α_{1} Tasa de ingestión del zooplancton por consumir fitoplancton

2.00 1/T

α_{2} Tasa de ingestión de la anchoveta por consumir fitoplancton

0.70 1/T

α_{3} Tasa de ingestión de la anchoveta por consumir zooplancton

0.90 1/T

β_{1} Tasa de beneficio del zooplancton por consumir fitoplancton

1.55 1/T

β_{2} Tasa de beneficio de la anchoveta por consumir fitoplancton

0.50 1/T

β_{3} Tasa de beneficio de la anchoveta por consumir zooplancton

0.70 1/T

k_{1} Constante media de saturación del fitoplanctonzooplancton

0.80 T n

k_{2} Constante media de saturación del fitoplanctonanchoveta

0.80 T n

k_{3} Constante media de saturación del zooplanctonanchoveta

0.80 T n
μ_{1}  Tasa de mortalidad del fitoplancton  0.02 1/T
μ_{2}  Tasa de mortalidad del zooplancton  0.03 1/T
μ_{3}  Tasa de mortalidad de la anchoveta  0.25 1/T

De aquí, se podrá realizar la simulación computacional del modelo (1) con los valores de las Tablas 1 y 2, y deducir interpretaciones ecológicas de la dinámica que se presenta en la cadena alimenticia [17, 29, 32]. Se determinará un tiempo de simulación de 300 días, para conocer el comportamiento de las poblaciones, su respectiva evolución cuando interactúan entre sí [2, 25]. Por consiguiente, se realizará su interpretación ecológica respectiva para poder evaluar la situación que se suscita en el ecosistema marítimo peruano [22, 36].

3.1 Simulación computacional

El modelo presentado en (1) donde se ha garantizado la existencia, unicidad, y positividad de soluciones (Teoremas 1 y 2) [18, 21], y posee soluciones acotadas en el largo plazo (Teorema 3) [28, 36]; y estabilidad en el punto de equilibrio asociado (Teorema 4) [22, 33]. Además, el sistema (1) se ha dividido en subsistemas (19) y (20) para determinar si poseen órbitas periódicas en el largo plazo (Proposiciones 1 y 2) [30, 41].

Por ende, al poseer los valores numéricos de las tasas ecológicas asociadas al modelo podremos analizar y verificar la estabilidad del modelo y contrastar con la simulación computacional de acuerdo con las Tablas 1 y 2 [7, 16, 36]:

\{ \begin{matrix}\frac{dF}{dt} = 1.2F (1 − \frac{F}{0.60}) − \frac{2F ·Z}{0.80+F} − \frac{0.70F ·A}{0.80+F} − 0.02F, & F (0) = 1.30 \\ \frac{dZ}{dt} = \frac{1.55F ·Z}{0.80+F} − \frac{0.90Z·A}{0.80+F} − 0.03Z, & Z(0) = 0.70 \\ \frac{dA}{dt} = \frac{0.50F ·A}{0.80+F} + \frac{0.70Z·A}{0.80+Z} − 0.25A, & A(0) = 0.05. \end{matrix}

(21)

De acuerdo con el modelo propuesto (1), donde se ha introducido los valores de las Tablas 1 y 2, se obtuvo el sistema (21) para poder analizar su comportamiento numérico con interpretaciones ecológicas [27, 36]. De aquí, se puede obtener los puntos de equilibrio no negativos en el estado estacionario para poder verificar la situación ecológica que se tiene; y se pueda complementar con la simulación computacional en un tiempo de simulación determinado [33, 36].

Por consiguiente, se ha obtenido cuatro puntos de equilibrio, y además se ha determinado su condición en su estado estacionario respectivo con el apoyo del Teorema 4 (criterio de Routh-Hurwtiz) para nuestro modelo [1, 7].

3.1.1 Punto de extinción total

Este punto de equilibrio, w_{1} = (0; 0; 0), representa la situación menos deseada para el ecosistema marino, y donde la intervención humana puede ser el principal gestor de generación de este escenario crítico [25]. De acuerdo a nuestro modelo (21), el estado estacionario asociado a este punto es Inestable [15, 26].

3.1.2 Punto de la supervivencia del fitoplancton

Este punto de equilibrio, w_{2} = (0.5899; 0; 0) representa la situación que no se desea en el ecosistema marino, aún cuando la especie del fitoplancton se mantiene en el tiempo. La intervención humana puede ser el principal gestor de generación de este escenario crítico. De acuerdo a nuestro modelo (21), el estado estacionario es Inestable [5, 23, 25].

3.1.3 Punto de la extinción de la anchoveta

Este punto de equilibrio, w_{3} = (0.01579; 0.46843; 0) representa la situación que no se desea en el ecosistema marino, aún cuando la especie del fitoplancton y la especie del zooplancton se mantienen en el tiempo. La intervención humana puede ser el principal gestor de generación de este escenario crítico debido a la actividad pesquera sin ninguna consideración. De acuerdo a nuestro modelo (21), el estado estacionario es Inestable [5, 25, 38].

3.1.4 Punto de coexistencia

Este punto de equilibrio, w_{4} = (0.3297; 0.1397; 0.4410) representa la situación que siempre debe permanecer en el ecosistema marino. La intervención humana puede ser el principal gestor de romper esta coexistencia debido a un abuso de la actividad pesquera sin ninguna consideración de la preservación de la especie. De acuerdo a nuestro modelo (21), el estado estacionario donde se encuentra la coexistencia es Estable [7, 23, 32].

En resumen, se puede expresar la condición de cada punto de equilibrio del modelo (21) en la siguiente tabla de acuerdo a los valores numéricos [1, 22, 36].

Tabla 3: Tabla de estabilidad del modelo en su punto de equilibrio.

Punto de Equilibrio  Estado Estacionario
w_{1} = ( 0 ; 0 ; 0 )  Inestable
w_{2} = ( 0.5899 ; 0 ; 0 )  Inestable
w_{3} = ( 0.01579 ; 0.46843 ; 0 )  Inestable
w_{4} = ( 0.3297 ; 0.1397 ; 0.4410 )  Estable

Hay una consideración que se tiene con el modelo presentado (1) donde se obtuvo las condiciones de existencias de los puntos de equilibrio (14) donde se tiene cinco puntos, pero al introducir los valores de los parámetros se obtiene cuatro puntos de equilibrio asociados al sistema (21) debido que el punto de equilibrio (e_{3}) no satisface la condición de existencia positiva [27, 36].

Por consiguiente, se realizará la simulación computacional donde podremos visualizar el comportamiento de las especies (fitoplancton, zooplancton y anchoveta) en el tiempo. Con la simulación en el tiempo determinado, se podrá analizar el comportamiento en el corto plazo, mediano y largo plazo (análisis computacional) de acuerdo a las tasas ecológicas constantes consideradas en las tablas 1 y 2 [27, 36]. Además, se busca complementar el análisis cualitativo realizado en la sección 2, y también el desarrollo en la subsección 3.1 con los valores numéricos de los parámetros (constantes) y las poblaciones iniciales del modelo matemático de la cadena alimenticia de tercer eslabón [7, 22, 36].

Para la presente simulación, se recurrió al método numérico de Runge-Kutta de orden cuatro con un margen de error absoluto (10^{−6}) para tener una buena aproximación de la solución del sistema (21). Este método que se realiza por el tiempo discreto ayuda enormemente para la solución numérica en cuanto al costo computacional y el tiempo de compilación [32, 33].

Figura 2: Simulación computacional del modelo de cadena alimenticia.

Fuente: Elaboración propia.

La simulación computacional representa al modelo matemático (1) que comprende la dinámica que existe entre las especies marinas dentro del ecosistema y están expresadas en la Figura 1; con lo cual, la Figura 2 busca el escenario propuesto de acuerdo a los parámetros indicados en la Tabla 2 donde se puede observar la coexistencia de las especies, como se mostró en la Tabla 3, donde el estado estacionario de estabilidad se daba en el punto de equilibrio de la coexistencia, w_{4} = ( 0.3297 ; 0.1397 ; 0.4410 ), y la simulación concuerda con el valor estacionario de cada especie marina [27, 33].

Siempre debe concordar el resultado del análisis cualitativo (sección 2) con la solución numérica que provee la simulación computacional. Lo cual, en la simulación computacional presenta el estado de coexistencia de las tres especies (equilibrio biológico) que siempre se busca que se mantenga pero lastimosamente siempre está en peligro de una depredación descontrolada [14, 29].

La consideración de la mortalidad de las anchovetas (μ_{3}) se encuentra en una situación delicada donde la intervención del hombre puede ayudar a que se mantenga el equilibrio ecológico en favor del plancton, como también generar la extinción de la especie [2, 25]; por eso siempre se debe considerar esta tasa como un parámetro sensible que busque la coexistencia de las tres especies marítimas de manera que las cuatro poblaciones se vean beneficiadas (la población humana es la cuarta población que interviene en la cadena alimenticia que no aparece en el modelo matemático (1) de manera explícita) [5, 14, 38].

Figura 3: Simulación computacional del modelo (puntos de equilibrios).

Fuente: Elaboración propia.

En la Figura 3 se puede analizar el comportamiento de las tres poblaciones a través del tiempo de simulación donde el punto negro marca la población inicial de las especies marítimas (1.3, 0.7, 0.05). Y los puntos rojos son los puntos de equilibrio inestables (w_{1}, w_{2}; w_{3}) del sistema (21), y el punto verde es el punto de equilibrio estable (w_{4}) en el largo plazo [5, 22]. En el corto plazo se puede visualizar un comportamiento oscilatorio pero luego se estabiliza, y las poblaciones se mantienen constantes mientras no haya una pertubación externa que pueda afectar el ecosistema [27, 36]. Por lo tanto, las Figuras 2 y 3 presentan de una forma más detallada el comportamiento de las poblaciones dentro del ecosistema tanto en su plano a través del tiempo como para la visualización del punto estacionario estable [16, 23, 41]. Este estado en un ecosistema donde intervienen otras especies y la actividad pesquera no tienden a generar un estado estable pero esta simulación nos brinda información de como se puede bosquejar el comportamiento en el tiempo de nuestro modelo [6, 17, 40].

La anchoveta juega un rol fundamental e importante en la compleja red alimenticia del ecosistema de Humboldt, dado que cientos de diferentes especies de peces, mamíferos y aves marinas dependen de ella para su alimentación [29, 38]. Sin lugar a duda, siempre se debe considerar la manera que se pueda generar el aprovechamiento de los recursos naturales, y que el ecosistema se mantenga en un comportamiento adecuado de producción para la población humana, sin generar caos ni crisis en las demás especies de manera específica el fitoplancton, el zooplancton y la anchoveta [5, 14, 29]. Y este lineamiento viene ser la idea central de la investigación de como ir construyendo modelos que permitan conocer el comportamiento que se puede generar ante intervenciones de la pesquería en el ecosistema marítimo [3, 9, 15].

Figura 4: Plano fase de las especies marítimas. Fuente: Elaboración propia.

El gráfico de la Figura 4 nos ayuda a visualizar el comportamiento entre dos especies es el plano fase que permite analizar la evolución que realizan ambas en el tiempo sin tener el eje temporal. La Figura 4 presentará el plano fase de cada especie considerada en el modelo (21), de manera que nos ayude a ver la competencia que se suscita en la convivencia entre la anchoveta y el plancton (fitoplancton y zooplancton) [14, 36]. Por lo cual, un análisis más estructurado de cómo interactúan las especies dentro del tiempo de simulación considerado nos permite ver la evolución entre un versus de especie a especie. Y es, en la Figura 4 donde se puede observar la estabilidad local (punto atractor) desde el punto de vista del plano fase mostrado [29, 36, 38]. Teniendo las simulaciones computacionales en las Figuras 2 y 4 se puede observar como se presenta un comportamiento oscilatorio en las poblaciones sin llegar a la extinción, además que en un tiempo de mediano plazo las poblaciones se llegan a estabilizar generándose la coexistencia marítima que se obtiene de acuerdo al punto de equilibrio (w_{4}). Además de las Proposiciones 1 y 2 complementan bien a la Figura 4 debido que no se muestra la existencia de órbitas periódicas en el largo plazo [30, 41].

De acuerdo con la Figura 2, se presenta un comportamiento oscilatorio entre las especies marítimas donde la población de la anchoveta comienza a crecer lentamente hasta que se superpone a las otras dos especies (plancton) donde todas las especies se mantienen en un estado de coexistencia estable en el tiempo donde el teorema 4 respalda estos resultados numéricos para poder contrastar e interpretar de diferentes perspectivas de análisis [18, 27]. De forma adicional, se puede apreciar en la Figura 4 una proyección mediante regresión lineal como previsión en el comportamiento generado por la depredación tanto animal (explícita) como humana (implícita), esto quiere decir, la interacción de la presa con la depredadora en cada nivel de la cadena alimenticia [14, 29].

En la Figura 5, presentamos dos simulaciones que se han realizado con los valores de las Tablas 1 y 2, donde solamente se ha perturbado la tasa de mortalidad de la anchoveta (μ_{3} = 0.25), y como se puede observar, el nivel de perturbación en el equilibrio ecológico que puede ocasionar leves cambios en su tasa de mortalidad (causas naturales o actividad pesquera) [25, 29]. Para un valor de mortalidad de μ_{3} = 0.20, se obtiene un equilibrio ecológico en el tiempo pero con una menor población de zooplancton (casi a extinguirse) que en el escenario presentado en la Figura 2 donde tenía una tasa de mortalidad inicial; lo cual no ayuda al equilibrio ecólogico para que se mantenga las especies en un nivel poblaciones adecuado de manera que otras especies se vean afectadas por la ausencia del zooplancton y la depredación se centre únicamente en el fitoplancton lo cual en el tiempo podría extinguirse también [15, 17].

Por otro lado, si se incrementa levemente la tasa de mortalidad a un valor

Figura 5: Simulación computacional del modelo de cadena alimenticia.

Fuente: Elaboración propia.

de μ_{3} = 0.265, se obtiene una considerable reducción de la anchoveta en el corto y mediano plazo, lo cual no es bueno para el ecosistema marítimo, pero también para la actividad pesquera que no tendrá esta especie para seguir pescando, cuando la misma actividad pesquera es una de las principales causas de producir la extinción de esta especie [15, 17]. Por consiguiente, la pesca es una actividad sensible que puede ayudar a las especies como también puede eliminar a las mismas; esto conlleva a una conciencia ecológica para mantener el ecosistema por el bien de todas las especies (especies marítimas que depredan a la anchoveta (eslabones superiores de la cadena alimenticia) como la actividad pesquera del hombre) del mar [2, 25, 38]. Pero de acuerdo al modelo propuesto en (21) en el largo plazo la población de la anchoveta retoma su densidad donde comienza a crecer y estabilizarse [6, 36]. Esta simulación viene a ser algo alejada de la realidad debido que la consideración de tasas constantes no aborda la evolución de las poblaciones, como se reporta en los informes del IMARPE, de forma más cercana que se manifiesta en el mar donde interactúan especies marítimas y el hombre (actividad pesquera) [15, 29, 40].

Un comportamiento que llama la atención el lado derecho de la Figura 5 es una evolución oscilatoria del zooplancton en mayor densidad que el fitoplancton debido a la ausencia de la anchoveta [6]. Esto se debe a comportamientos de ausencia del depredador mayor genera causa que la respuesta funcional Holling Tipo II nos permita analizar sobre la cantidad de presas que puede consumir el depredador restante [27, 31]. Además, la actividad pesquera se ve colapsada en el mediano plazo debido a un abrupto realizado en el corto plazo que conlleva a un desequilibrio biológico tanto para la cadena alimenticia considerada de tres niveles, sino también para una cadena mucho mayor [15, 29, 38].

En este sentido, el IMARPE custodia las temporadas de vedas y controles respectivos para que no suceda este escenario que sería catastrófico para todo el ecosistema marítimo como el sistema económico que se desprende de la actividad pesquera [16, 29]. Y por lo cual, se necesita proyecciones de la evolución de las poblaciones en diferentes situaciones críticas que pueda perjudicar el ecosistema marítima tanto para el desarrollo natural de las especies marítimas y la intervención del hombre sin control [14, 15].

4 Discusiones y conclusiones

El modelo matemático de la cadena alimenticia de la anchoveta que se ha presentado en la sección 1, y se ha estudiado en la sección 2, donde se ha considerado las Ecuaciones Diferenciales Ordinarias (1) con una respuesta funcional Holling Tipo II que nos ayuda a describir la dinámica que existe entre las tres especies, dentro de un ecosistema marítimo que permita analizar la evolución en el tiempo para ver la mejor manera de preservar la coexistencia de las tres especies [7, 36]. Además se ha realizado la estabilidad del sistema en el estado estacionario de la coexistencia (w_{4}) que es donde se encuentra el equilibrio ecológico teniendo en cuenta la intervención del hombre de manera implícita como una de las causas principales de la mortalidad de la anchoveta [15, 23, 27]. La consideración de la función de Dulac (criterio de Dulac-Bendixon) en la sección 2.4 donde el sistema (1) se desprendió dos subsistemas (19) y (20) donde se probó que no poseían órbitas cerradas dentro de sus respectivas regiones [30, 41].

Además, en la sección 3, de manera específica en la subsección se centra en las simulaciones computacionales que nos permiten analizar el comportamiento de las poblaciones en el tiempo donde el enfoque es mantener un nivel adecuado que la coexistencia mantenga en el equilibrio ecológico donde la intervención del hombre no cause un rompimiento de este equilibrio porque él mismo crearía la ruptura de su propia producción pesquera [15, 25, 38].

También, se ha presentado variaciones, en la tasa ecológica sensible donde el hombre interviene, siendo la tasa de mortalidad de la anchoveta (μ_{3} = 0.25), y al disminuirse esta tasa en cinco unidades porcentuales (μ_{3} = 0.20) causa de una decrecimiento en la población del zooplancton; y al perturbar un leve aumento porcentual (μ_{3} = 0.265) genera la extinción de la anchoveta en el tiempo [38]. La impotancia del análisis cualitativo que se complementa con la simulación numérica representa una gran ayuda a comprender la sensibilidad de los parámetros que indica los teoremas de estabilidad (Teorema 4), para que obtener la situación ecológica más adecuada [5, 12, 17, 36].

Se tiene que considerar muchos detalles en los parámetros para el modelo (1) que nos permita poder realizar simulaciones más cercanas a la realidad que se vive en el especio marítimo perteneciente al país donde los estudios, en muchos casos, no alcanza para concretar la verdadera situación de esta cadena alimenticia (fitoplancton, zooplancton y anchoveta) [17, 36]; donde la cadena alimenticia se encuenta en el tercer eslabón, más aún cuando la fauna marina es mucho más extensa de la considerada en el modelo, y los parámetros no necesariamente deben ser constantes sino que evolucien en el tiempo [5, 22], y la más importante de consideración debería ser una expresión matemática para representar la dinámica de la pesquería y cómo afecta al ecosistema marítimo donde se podría perturbar el modelo de tal manera que representa la dinámica de la depredación [5, 15, 25].

Finalmente, una de las discusiones más importante sería una construcción de la cadena alimenticia más completa para poder analizar la sensibilidad de los parámetros de mortalidad tiene dos acciones de generarse por la depredación de las especies marítimas y por la depredación del ser humano [2, 15, 25]; más aún conocer las condiciones para garantizar la coexistencia de las poblaciones [7, 28]. La actividad pesquera es una de las razones que el equilibrio ecológico puede perderse cuando se realiza sin ninguna vigilancia y control [14, 29]. Por lo cual, siempre se debe estudiar los impactos que se genera en toda la cadena alimenticia para que sea provechoso para todas las especies y no se rompa el equilibrio ecológico de la fauna marina en provecho de todo el ecosistema [23, 38].

Por lo cual, podemos obtener las siguientes conclusiones.

1. El modelo matemático presentado en (1), es una extensión del modelo desarrollado en [36], donde se focaliza la dinámica de alimentación que realiza la anchoveta (tercer eslabón), siendo el depredador mayor del ecosistema estudiado, con respecto al fitoplancton (primer eslabón) y al zooplancton (segundo eslabón) como se muestra en la Figura 1 [29]. Esto permite modelar de manera más real la cadena alimenticia generada en el mar peruano y analizar la dinámica depredadora de la anchoveta [14, 23].

2. Los teoremas de acotación de soluciones (Teorema 3) y de estabilidad (Teorema 4) ayudan a visualizar el estado estacionario que se genera de acuerdo las restricciones para garantizar la no negatividad de los puntos de equilibrio [17, 36]. Así, se pueden buscar las condiciones sensibles para lograr la coexistencia de las especies desde el enfoque matemático [13, 24].

3. Las proposiciones 1 y 2 presentan a los subsistemas (19) y (20) no poseen órbitas periódicas cerradas respectivamente en sus regiones [30, 41]. Esto nos permite conocer un bosquejo del comportamiento del sistema (1) en el tiempo. No obstante, el sistema podría admitir una bifurcación de Hopf en la coexistencia de las especies como lo indican Samares y Anal [36].

4. La simulación computacional permitió conocer el comportamiento de las especies en el tiempo donde en el inicio se presentó crecimiento oscilatorio pero luego se llega a estabilizar las tres especies, sobretodo la población de la anchoveta que se mantiene por encima de las otras dos y coexisten en el tiempo [15, 17, 29].

5. La Figura 4 presenta el comportamiento en el plano fase de cada especie para visualizar el comportamiento oscilatorio y luego converger a un punto como complemento a la Figura 2, donde se analiza la evolución en el tiempo de las tres especies donde se analiza como llegan al estado estacionario de coexistencia [23, 36]. Este estado estacionario se obtiene ante la sensibilidad de los parámetros que se tiene en la recolección de los datos por parte del IMARPE [5, 15].

6. La mortalidad de la anchoveta dentro de la cadena alimenticia presentada viene a ser un parámetro muy sensible [5] debido a la actividad pesquera que se incluye de manera implícita dentro del modelo (1), siendo una debilidad del modelo pero puede ser el inicio de la construcción de un modelo más robusto y que considere la dinámica de la pesquería de forma explícita [14, 36].

7. La pesquería puede ser la acción para romper el equilibrio ecológico si no se tiene intervenciones adecuadas para preservar la dinámica de depredación marítima [5, 16], y producción pesquera como se presenta en la Figura 5, donde la acción de pescar menos no beneficia al zooplancton y pescar un poco más puede generar la extinción de la anchoveta [25, 38].

Agradecimientos

Nos gustaría agradecer a cada árbitro anónimo por la cuidadosa lectura del manuscrito, y por brindarnos comentarios constructivos que ayudaron a mejorar el mismo. Además, de forma personal quisiéramos agradecer por el apoyo, guía y sugerencias a nuestra compañera Alejandra Melchor Barrionuevo, y a la profesora Roxana López Cruz por sus sugerencias académicas para el desarrollo de esta investigación.

Financiamiento

La presente investigación no ha recibido ningún tipo de financiamiento para su realización.

Referencias

[1] B. Aguirre-Hernández, C.A. Loredo-Villalobos, E.C. Díaz-González, E. Campos-Cantón, Estabilidad de sistemas por medio de polinomios Hurwitz, Revista de Matemática: Teoría y Aplicaciones 24(2017), no. 1, 61–77. Doi: 10.15517/RMTA.V24I1.27751

[2] H. Aquije-Ballon, L. Ortíz-Chávez, Relación de la pesca industrial de anchoveta con la pesca artesanal: evidencia del caso peruano, Facultad de Ciencias Económicas y Empresariales, Universidad de Piura, 2020. Recuperado de https:// pirhua. udep. edu.

pe/bitstream/ handle/ 11042/ 3110/ ECO-L_004. pdf?
sequence=1&  isAllowed=y

[3] E.M. Aranguren-Lizarzaburu, K. Arteaga-Zarate, F.R. Chavez-Mapelli, Impacto ecónomico de la implementación de los límites máximos de captura por embarcación en la industria pesquera del Perú entre los años 2005 hasta el 2014, Facultad de Hotelería y Negocios, Universidad Peruana de Ciencias Aplicadas, 2016. Recuperado de https://

repositorioacademico.   upc. edu. pe/bitstream/ handle/
10757/ 622221/ chavez_fm.pdf? sequence=5&  isAllowed=y

[4] R. Bamón, Ecuaciones Diferenciales y Sistemas Dinámicos, Pro Mathematica 5(1991), no. 9-10, 145-168, 1991. Recuperado de http://

revistas. pucp. edu. pe/index. php/ promathematica/
article/ download/ 3419/ 3266/

[5] R.C. Barret, J.P. Caulkins, A.J. Yates, D.L. Elliot, Population dynamics of the Peruvian anchovy, Mathematical Modelling 6(1985), no. 6, 525–548. Doi: 10.1016/0270-0255(85)90052-1 [6] M. Bouchon-Corrales, S. Cahuín-Villanueva, E. Díaz-Acuña, M. Ñiquen-

Carranza, Captura y esfuerzo pesquero de la pesquería de anchoveta peruana (Engraulis Ringens), Boletín Instituto del Mar del Perú 19(2000), no.

1-2, 109–115. Recuperado de https:// revistas. imarpe. gob.

pe/index. php/ boletin/ article/ view/ 159/ 148

[7] A. Boussouar, S. Le Bihan, O. Arino, P. Prouzet, Mathematical model and numerical simulations of the migration and growth of Biscay Bay anchovy early larval stages, Oceanologica Acta 24(2001), no. 5, 489-504. Doi: 10.1016/S0399-1784(01)01167-7

[8] D. Driver-Rodney, Ordinary and Delay Differential Equations, Springer, New York NY, 1977. Doi: 10.1007/978-1-4684-9467-9

[9] E. Galarza, N. Collado, Los derechos de pesca: el caso de la pesquería de anchoveta peruana, Apuntes, Revista de Ciencias Sociales 40(2013),

no. 73, 7–42. Recuperado de https:// dialnet. unirioja. es/

servlet/ articulo? codigo=4786988

[10] E. Galarza, J. Kámiche-Zegarra, Pesca artesanal: oportunidades para el desarrollo regional, Universidad del Pacífico, Repositorio Institucional.

Recuperado de http:// hdl. handle. net/ 11354/ 1013

[11] E. González-Olivera, J. Mena-Lorca, Análisis cualitativo de un modelo de pesquerías de acceso abierto, Investigaciones Marinas 22(1994), 3–11. Doi: 10.4067/S0717-71781994002200001

[12] S.A. Gourley, K. Yang-Kuang, A stage structured predator-prey model and its dependence on maturation delay and death rate, J. Math. Biol. 49(2004), 188–200. Doi: 10.1007/s00285-004-0278-2

[13] A. Guillermo, Tópicos de Biología Matemática. Notas de Clase, Instituto de Balseiro, Universidad Nacional de Cuyo-CNEA, 2018. Recuperado de

http:// ricabib. cab. cnea. gov. ar/680/ 2/biologicos_
2.pdf

[14] Instituto del Mar del Perú, Biología anchoveta, Reportes del Mar Peruano, Ministerio de Producción, Perú, 2020 [citado 1 septiembre 2020]. Recupe-

rado en http:// www. imarpe. gob. pe/imarpe/ lista__. php?

id_seccion=I0170010402000000000000

[15] Instituto del Mar del Perú, Pesquería anchoveta, Reportes del Mar Peruano, Ministerio de Producción, Perú, 2020 [citado 1 septiembre 2020].

Recuperado en http:// www. imarpe. pe/imarpe/ lista__.

php? id_seccion=I0170010412000000000000

[16] Instituto del Mar del Perú. Situación del stock norte-centro de la anchoveta peruana (Engraulis ringens) al mes de mayo de 2020 y perspectivas de explotación para la primera temporada de pesca del año, Ministerio de Producción, Perú, 2020. Recuperado en http:// www. imarpe. gob.

pe/imarpe/ descarga_informe. php? archivo=107

[17] V.F. Krapivin, The estimation of the Peruvian current ecosystem by a mathematical model of biosphere, Ecological Modelling 91(1996), no. 1-3, 1–14. Doi: 10.1016/0304-3800(95)00155-7

[18] Y. Kuang, Basic properties of mathematical population models, Department of Mathematics and Statistics, Arizona State University, 2002.

Recuperado de http:// citeseerx. ist. psu. edu/ viewdoc/

download? doi=10. 1.1.9.3382& rep=rep1& type=pdf

[19] M. Li, Dulac criteria for autonomous systems having an invariant affine manifold, Journal of Mathematical Analysis and Applications 199(1996), 374–390. Doi: 10.1006/JMAA.1996.0147

[20] J.A. Lobón-Cerviá, Dinámica de poblaciones de peces en ríos: pesca eléctrica y métodos de capturas sucesivas en la estima de abundancias (Vol. 3), Museo Nacional de Ciencias Naturales, Edi-

torial CSIC, 1991. in: https:// bibliotheques.  mnhn. fr/

medias/ detailstatic.  aspx? INSTANCE=EXPLOITATION&
RSC_BASE=HORIZON&  RSC_DOCID=423056

[21] R. López-Cruz, Structured SI Epidemic Models with Applications to HIV Epidemic, Ph.D. Diss. Arizona State University; 27–45, 2006. In:

https:// math. la.asu. edu/ ~kuang/ paper/ Roxana. pdf

[22] D. Manna, A. Maiti, G.P. Samanta, A Michaelis-Menten type food chain model with strong Allee effect on the prey, Applied Mathematics and Computation 311(2020), no. 5, 390–409. Doi: 10.1016/j.amc.2017.05.040

[23] M. Medina, L. Herrera, J. Castillo, J. Jaque, N. Pizarro, Alimentación de la anchoveta (Engraulis ringens) en el norte de Chile (18^{o}25’-25^{o}40’S) en diciembre de 2010. Lat. Am. J. Aquat. Res. 43(2015), no. 1, 46–

58. Recuperado de https:// www. redalyc. org/ articulo. oa?
id=175036979005

[24] E. Mendoza, O. Osuna, G. Villavicencio-Pulido, Funciones de Dulac para modelos matemáticos de la ecología, Rev. de Matemática: Teoría y Aplicaciones 27(2020), no. 2, 367–382. Doi: https:// doi. org/ 10.

15517/ rmta. v27i2. 34430

[25] Ministerio de Producción, Reporte mensual de la evolución del sector pesca – Junio 2020, Estudios Económicos, Estadística Pesquera.

Recuperado de http:// ogeiee. produce. gob. pe/index. php/

shortcode/ estadistica-oee/  estadistica-pesquera

[26] F.R. Momo, A.F. Capurro, Ecología Matemática, Principios y Aplicaciones, Ediciones Corporativas, 1ra. Ed., 2006.

[27] J.D. Murray, Mathematical Biology I. An Introduction, 3rd Ed., Interdisciplinary Applied Mathematics 17, Springer, New York NY, 2002. Doi: 10.1007/b98868

[28] B. Mukhopadhyay, R. Bhattacharyya, A stage-structured food chain model with stage dependent predation: Existence of codimension one and codimension two bifurcations, Nonlinear Analysis: Real World Applications 12(2011), no. 6, 3056–3072. Doi: 10.1016/j.nonrwa.2011.05.007

[29] OCEANA, La anchoveta y el Niño, Reporte, Lima-Perú, 2020 [citado 1 septiembre 2020]. Recuperado de https:// peru. oceana. org/ es/

la-anchoveta-y-el-nino

[30] O. Osuna, C. Vargas-De-León, Construction of Dulac functions for mathematical models in population biology, International Journal of Biomathematics 8(2020), no. 3, 1–20. Doi: 10.1142/S1793524515500357

[31] C.E. Paredes, Eficiencia y equidad en la pesca peruana: la reforma y los derechos de pesca. Reporte del Consorcio de Investigación Económica y Social (CIES), Instituto del Perú, Universidad de San Martín de Porres,

2012. Recuperado de https:// cies. org. pe/sites/ default/

files/ investigaciones/  eficiencia_y_equidad_en_la_
reforma_de_la_pesca_en_el_peru. pdf

[32] J.R. Pérez-Nuñez, R. López-Cruz, Análisis matemático de una cadena alimenticia - presa - depredador - control - biológico, Selecciones Matemáticas 4(2020), no. 1, 112–123. Doi: 10.17268/sel.mat.2017.01.11 [33] N. Pino-Romero, Análisis y simulación numérica de un modelomatemático

SI con retardo discreto para las enfermedades de transmisión sexual, Tesis de Maestría en Matemática Aplicada, Facultad de Ciencias Matemáticas, Escuela de Posgrado, UNMSM, 2017. Recuperado de http://

cybertesis.  unmsm. edu. pe/handle/ 20.500. 12672/ 5735

[34] N. Pino-Romero, R. López-Cruz, Teoremas de estabilidad en un modelo matemático SI con dinámica vital estructurado por sexo para el estado libre de infección desarrollado mediante las ecuaciones diferenciales ordinarias y las ecuaciones diferenciales con retardo aplicado a la salud pública del Perú, Selecciones Matemáticas 4(2017), no. 2, 202–210. Doi: 10.17268/sel.mat.2017.02.07

[35] N. Pino-Romero, C.U. Salazar-Fernández, R. López-Cruz, Solución uniformemente acotada y estabilidad asintótica del punto libre de infección de un modelo matemático SI con dinámica vital (crecimiento logístico) mediante las ecuaciones diferenciales con retardo, Selecciones Matemáticas 6(2019), no. 1, 66–76. Doi: 10.17268/sel.mat.2019.01.09

[36] S. Pal, A. Chatterjee, Dynamics of the interaction of plankton and planktivorous fish with delay, Cogent Mathematics, 2(2015), no. 1, article 1074337. Doi: 10.1080/23311835.2015.1074337

[37] Significados [Internet], Significado de cadena alimenticia acuática, Ciencia [citado 2 de septiembre 2020]. Recuperado de https:// www.

significados.  com/ cadena-alimenticia-acuatica/

[38] Sociedad Nacional de Pesquería, Industria pesquera: Contribución a la economía peruana, Lima-Perú, 2020 [citado 1 septiembre 2020]. Recupe-

rado de https:// www. snp. org. pe/relevancia-economica/

[39] D.A. Valencia-Talavera, Ecología trófica de la cachema Cynoscion analis (Jenyns, 1842) procedente de la pesca artesanal en las zonas de Santa Rosa, Chimbote y Callao, Tesis de Licenciatura, Facultad de Ciencias Veterinarias y Biológicas, Universidad Científica del Sur,

2020. Recuperado de https:// repositorio.  cientifica.  edu.

pe/handle/ 20.500. 12805/ 1148

[40] J. Vandermeer, Omnivory and the stability of food webs, Journal of Theoretical Biology 238(2006), 497–504. Doi: 10.1016/j.jtbi.2005.06.006

[41] V. Vásquez-Hipólito, Modelación matemática de la interacción plantapolinizador y bifurcación hacia atrás en un modelo epidemiológico con fecundidad aumentada, Tesis de Doctorado en Matemáticas Aplicadas, Centro de Investigación en Matemáticas (CIMAT), 2016. Recu-

perado de https:// cimat. repositorioinstitucional.    mx/

jspui/ bitstream/ 1008/ 394/ 2/TE%20565. pdf

[42] G. Villaseñor-Aguilar, Sobre las órbitas periódicas de ecuaciones diferenciales en el plano, Tesis de Maestría, Instituto de Física y Matemáticas, Universidad Michoacana de San Nicolás de Hidalgo, 2010. Recuperado de

http:// bibliotecavirtual.   dgb. umich. mx:8083/ xmlui/
bitstream/ handle/ DGB_UMICH/ 1130/ IFM-M-2010-0001.
pdf? sequence=1