Preventing extinction in Rastrelliger brachysoma using an impulsive mathematical model

1.   Introduction

One of the most important foods for human is fish. It is about $35\%$ of protein humans consumed. The South Pacific jack mackerel (JM) (Trachurus murphyi) is one of the most important fish for consumption in the world ^[6]. The Atlantic Mackerel (Scomber scombrus) is one of the most fish species in the North Atlantic which is widely spawned to produce human's food ^[10]. The India mackerel (Rastrilliger kangurta) ^[4] and Spanish mackerel (Scomberomorus commerson) ^[17] are also the fish production using for human consumption. In the Asia-Pacific region, Indo-Pacific mackerel (Rastrelliger brachysoma) or short mackerel is a fish indigenous in this region. This species spread almost in coasts and islands in the Gulf of Thailand and the Andaman sea ^[5]. It is about 14–20 centimeters long. It can spawn about $20,000$ eggs at one time. Short mackerel is a pelagic fish that is highly economically important in the region because it is cheap and tasty ^[12]. It is therefore widely consumed in this region. So, there is high demand in South East Asia, resulting in numerous fisheries.

In 1995, the Department of Fisheries, Ministry of Agriculture and Cooperatives, Thailand reported that short mackerels spawn two periods which are January to March and June to August each year. Moreover, they spawn more than $20$ meters below sea level. Recently, the Department of Fisheries reports that the number of catching short mackerel drastically decreased every year. In 2014, on average 145,000 tons were caught per year 70,000 tons were caught in 2015, then 31,000 tons, 25,000 tons and 17,700 tons in 2016–2018, respectively ^[18]. The major problem is catching fish in the spawning season. The economic value loss caused by catching one ton of Indo-Pacific mackerel fry is about seven to eight million baht. The government proposed a law to protect fry from fisheries in 2015.

The mathematical model is used for preparing to predict future phenomena in real-world problems. The model allows simulating any possibility. Many mathematical models have been used to describe the dynamic of the fish population. Bueno et al. reviewed the mathematical model for managing the carrying capacity of aquaculture activities in lakes ^[3]. The model of fish population dynamic was proposed by Hallam et al. by using the ordinary differential equations ^[7]. The model studied the changes in lipid and structure components which considered the energy demand and available energy. Khatun et al. proposed the mathematical model to studied renewable fishery management ^[11]. The two new second-order characteristic scheme was proposed to solve an age-structured population model with nonlinear diffusion and reaction ^[15] which used to consider the environmental and spatial region effect. Raymond et al. proposed the model to describe the dynamics of a two-prey one predator system in fishery ^[21].

In addition, the toxicity in the sea is one of the problems which affect the population of short mackerels and their food. Bergami et al. studied the effect of plastic pollution in the marine ecosystem especially nanoplastics ^[2] to the plankton species which is the food of many fishes in the marine. Hashiguchi et al. identified and evaluated the toxicity of palm oil mill which affected the plankton species ^[8]. The antifouling compound zinc pyrithione (ZPT) was studied the effect on the natural planktonic communities by Hjorth et al. ^[9]. They found evidence of the diverse effect of ZPT on marine plankton. Kaur et al. proposed the dynamical model to study the effect of toxicity to the plankton system ^[19]. Randall and Tsui studied the effect of ammonia in the aquatic environment on the central nervous system of fish ^[20].

To prevent the extinction of Indo-Pacific mackerels, we propose a mathematical model to find a way to control the catching of the Indo-Pacific mackerels by using the impulsive model. The novelty of this work is the prey-predator model of Indo-Pacific mackerels (short mackerels) and their foods which composes the toxic environment and impulsive fisheries. The impulsive model is a suitable model for fisheries because the fishermen are allowed to harvest fish on some periods for preventing the extinction of fish. The proposed model can be forecast the population of short mackerel and short mackerel food (plankton and small fishes). We analyze the model to consider the behavior of the model solution. The numerical simulations are used to verify the theoretical results.

The rest of this paper is organized as follows: The modification model is proposed in Section 2, and model analysis is presented in Section 3. In Section 4, numerical simulations are showed, followed by conclusions in Section 5.

2.   Model modification

Now, an impulsive mathematical model was proposed in which the control of the extinction of the short mackerel by considering biological control, toxic environment, and impulsive fisheries are as follows:

For $t\neq nT$,

For $t = nT$,

where $x_1(t)$ is the density of small fishes (short mackerel food) population at time $t$, $x_2(t)$ is the density of the short mackerel population at time $t$, $T$ is the period of impulsive fisheries, $\gamma$ is the negative effect of fisheries on the density of small fishes population, and $\omega$ is the negative effect of fisheries on the density of short mackerel population with $n\in {\mathbb{Z}_{+}}$, ${\mathbb{Z}_{+}} = \{1, 2, 3, ...\}$, $0 < \gamma < 1$, and $0 < \omega < 1$. All parameters in the model are non-negative.

Equation (2.1) expresses the rate of change of the population density of short mackerel food (small fishes). The small fishes increase by the logistic function with the growth rate $a_1$ and the carrying capacity $k_1$. On the other hand, the density of them decreases due to being hunted by short mackerel with rate $b$, natural death with the rate $d_1$, and toxic death with the rate $u_1$.

Equation (2.2) expresses the rate of change of the population density of the short mackerel (Rastrelliger brachysoma). The short mackerels increase by the logistic function with the growth rate $a_2$ and the carrying capacity $k_3$. On the other hand, the density of them increases due to hunting small fishes with the rate $r$, natural death with the rate $d_2$, and toxic death with the rate $u_2$.

3.   Model analysis

Let

where ${{\mathbb{R}}_{+}} = [0, \infty), {\mathbb{R}}_{+}^{2} = \{X\in {{\mathbb{R}}^{2}}:X = (x_1, x_2), x_1\ge 0, x_2\ge 0\}$. The map defined by the right hand side of system (2.1)–(2.4) is denoted by $F = ({{F}_{1}}, {{F}_{2}})$.

Definition 3.1 (^[1]). The function $V$ defined in (3.1) is said to belong to class ${{V}_{0}}$ if

(a) $V$ is continuous in $(nT, (n+1)T]\times \mathbb{R}_{+}^{2}\to{{\mathbb{R}}_{+}}$ and for each $X\in \mathbb{R}_{+}^{2}, n\in {\mathbb{Z}_{+}},$ $\lim_{(t, Y)\to(nT^+, X)} V(t, Y) = V(nT^+, X)$ exists.

(b) $V$ is locally Lipschitzian in $X$.

Definition 3.2 (^[1]). Suppose $V\in {{V}_{0}}.$ For $(t, X)\in (nT, (n+1)T]\times \mathbb{R}_{+}^{2},$ the upper right derivative of $V(t, X)$ with respect to system (2.1)–(2.4) is defined by

where $F = ({{F}_{1}}, {{F}_{2}})$.

The solution of system (2.1)–(2.4), $X(t) = (x_1(t), x_2(t))$, is assumed to be a piecewise continuous function. It means that, $X(t)$ is continuous on $(nT, (n+1)T]$, $n\in {\mathbb{Z}_{+}}$ and $\lim_{t\to nT^+}X(t) = X(nT^+)$ exists. By the smoothness properties of $F$, the system (2.1)–(2.4) has a unique solution.

Lemma 3.3. Suppose $X(t) = (x_1(t), x_2(t))$ is a solution of the system (2.1)–(2.4) with the initial value $X({{0}^{+}})\ge 0$. Then the solution $X(t)\ge 0$ for all $t\ge 0$.

Proof. The solution $x_1(t)$ with non-negative initial value can be negative when the slope of $x_1(t)$ at $0$ is negative. So, the proof of this Lemma can be expressed as follows:

Case $t\ne nT$.

Consider the Eq (2.1)

Whenever $x_1(t) = 0$, the slope of $x_1(t)$ can be described by $\dfrac{dx_1}{dt} = 0$. This means that $x_1(t)$ cannot be negative. So, $x_1(t)$ is a non-negative solution.

Consider the Eq (2.2)

Whenever $x_2(t) = 0$, the slope of $x_2(t)$ can be described by $\dfrac{dx_2}{dt} = 0$. This means that $x_2(t)$, cannot be negative. So, $x_2(t)$ is a non-negative solution.

Case $t = nT$.

Consider the Eq (2.3)

Since $x_1 \ge 0$ and condition $0 < \gamma < 1$, we have that $x_1(t^+)\ge 0.$

Consider the Eq (2.4)

Since $x_2 \ge 0$ and condition $0 < \omega < 1$, we have that $x_2(t^+)\ge 0.$

Lemma 3.4. The solution $X(t) = (x_1(t), x_2(t))$ has upper bound i.e. $x_1(t)\le M$ and $x_2(t)\le M$, for sufficiently large $t$, provided that

Proof. We let $V(t) = x_1(t)+x_2(t)$,

and

Consider $t\ne nT$, we choose $\hat{c} > 0$ and

Then,

Hence $D^{+}V +\hat{c}V\le M_0.$

Consider $t = nT,$

By Lemma 2.2 of Liu et al. ^[16] we obtain that, for $t\in (nT, (n+1)T]$,

Since there exists $M > 0$ such that $x_1(t)\le M$ and $x_2(t)\le M$ for sufficiently large $t$ thus $V(t)$ is uniformly ultimately bounded.

3.1.   The reduced impulsive system

The reduced impulsive system of system (2.1)–(2.4) when the population density of small fishes is zero ($x_1 = 0$) is:

where $A \equiv \dfrac{a_2}{k_3}+u_2 > 0$ and $B \equiv a_2-d_2$.

We obtain $B > 0$ if

The solution of Eq (3.3) is

where $c$ is arbitrary constant.

By the conditions (3.4), (3.6) and $x_2$ is increasing function, a periodic solution of system (3.3)–(3.5) is

with $\tilde{x}_2(0^+) = \dfrac{B(1-\omega-e^{-BT})}{A(1-e^{-BT})} > 0.$

Therefore, the system (3.3)–(3.5) has the positive solution

Lemma 3.5. The periodic solution $\tilde{x}_2(t)$ of system $(3.3)–(3.5)$ exists and $x_2(t) \to \tilde{x}_2(t)$ as $t \to \infty$ for all solution $x_2(t)$ of system $(3.3)–(3.5)$. Hence,

is a periodic solution of the original system $(2.1)–(2.4)$ at the zero density of small fishes for $t \in (nT, (n+1)T]$ and

Theorem 3.6. Suppose

and

where

and

Then the solution $\left(0, \tilde{x_2}(t)\right)$ of the system $(2.1)–(2.4)$ is locally asymptotically stable.

Proof. Here, we focus on a small perturbation $\left(v_1(t), v_2(t)\right)$ from the point $\left(0, \tilde{x}_2(t)\right)$:

Then,

where $\Phi(t)$ satisfies

with the identity matrix $\Phi(0) = I$.

Therefore, the matrix of fundamental solution is

Note that the terms (*) and (**) are not necessary to be calculated because the further analysis does not require these terms.

Linearization of Eqs (2.3)–(2.4) yields

The eigenvalues of the following matrix $W,$

are

Since $0 < \gamma < 1, 0 < \omega < 1,$ and the conditions (3.6), (3.10)–(3.12) hold, it implies that

Therefore, $|\lambda_1| < 1$ and $|\lambda_2| < 1.$ We can conclude that the solution $\left(0, \tilde{x_2}(t)\right)$ of the system (2.1)–(2.4) is locally asymptotically stable by Floquet thery. The proof is completed.

3.2.   Permanence of system

Definition 3.7 (^[1]). The system $(2.1)–(2.4)$ is permanent if the solution is bounded. That is, there are positive constants $\bar{a}$ and $\bar{b}$ and a finite time ${t_0}$ such that for all solution with the positive initial values $x_1(0^+) > 0$ and $x_2(0^+) > 0$

for all $t > {t_0}$.

Theorem 3.8. The system $(2.1)–(2.4)$ is permanent if inequalities $(3.2)$, $(3.6)$, $(3.11)$, $(3.12)$ hold and the following conditions:

and

are satisfied.

Proof. By Lemma 3.4, there is a constant $M > 0$ such that the solution $x_1(t)\le M$ and $x_2(t)\le M$ for sufficiently large $t$.

Since $\dfrac{{r}{b}x_1x_2}{1+{k_2}x_1}\ge0$, Eq (2.2) implies that

then we obtain

for some positive $\epsilon$ and large enough $t$.

Hence,

for large enough $t.$

Hence, the remaining proof is the system has a lower bound. That is, there exists a positive constant ${m_2}$ such that $x(t) > {m_2}$. To obtain the result, we let

for some ${m_3} > 0$.

Next, there are two steps as follows:

Step 1. Prove by contradiction that there is $t_1$ such that $x_1(t_1)\ge{m_3}$. Suppose that $x_1(t) < {m_3}$ for all $t > 0$. From Eqs (2.2) and (2.4), we get

Let us consider the comparison system

and

Hence,

is a positive periodic solution of system (3.15)–(3.17) with

The system (3.15)–(3.17) has a positive solution

with $t\in (nT, (n+1)T]$ and $\dfrac{1}{Y(t)}\to\dfrac{1}{\tilde{Y}(t)} \quad$ as $t\to \infty$, where

By the comparison theorem in ^[14], we obtain that

Now, we consider (2.1)

Since $x_2(t)\leq Y(t)$, there is a ${T^*} > 0$ such that,

for a sufficiently small ${\epsilon_1} > 0$.

Therefore,

for $t\ne nT, t\ge {T^*}$ and

for $t = nT, t\ge {T^*}.$

Letting $N\in {\mathbb{Z}_+}$ and $NT\ge{T^*}$, and integrating over $(nT, (n+1)T], n\ge N,$ we get

where

Consider

For sufficiently small ${\epsilon_1} > 0$,

Since $\hat{M_1} < {a_1}-{d_1}$ and (3.11) is satisfied, a small positive ${m_3}$ is chosen so that $\ln\eta > 0$.

We get

We observe that $x_1((n+k)T)\ge x_1(nT)\eta^{k}\to \infty$ as $k\to \infty$ which contradicts the boundedness of $x_1(t).$ Therefore, there exists ${t_1} > 0$ such that $x_1(t_1)\ge{m_3}.$

Step 2. The proof is completed if $x_1(t)\ge{m_3}$ for all $t > {t_1}$. Otherwise, $x_1(t) < {m_3}$ for some $t > {t_1}$. Setting $t^* = \inf_{t > t_1}\{x_1(t) < m_3\}$. There are two cases as follows:

Case 1: ${t^*} = {n_1}T$ for some ${n_1}\in {Z_+}$. That is $x_1(t)\ge{m_3}$ for $t\in (t_1, t^*]$ and $x_1(t^*) = {m_3}$ by continuity of the solution $x_1(t)$.

Since $x_1(t) < M$ and $m_1 < x_2(t) < M$ for some positive $M$ and $m_1$ with sufficiently large $t$, we can choose $\bar{M} > 0$ and $\bar{m}_1 > 0$ so that

and

such that

Then, we choose $n_2, n_3\in{\mathbb{Z}_+}$ such that

and

where

Define ${T}' = {n_2}T+{n_3}T$. There exists $t_2\in (t^*, t^*+T']$ so that $x_1(t_2) > {m_3}$.

Otherwise, considering Eq (3.19) with

we obtain

for $t\in(nT, (n+1)T]$ where $n_1\le n\le n_1+n_2+n_3$.

For ${n_2}T\le t-{t^*}\le T'$, we have

Since the condition (3.14), we get

Then,

According to Step 1, we obtain

From Eq (2.1), we have

Integrating inequality (3.27) over $[t^*, t^{*}+{n_2}T]$, we have

hence,

It is in contradiction to the definition of $m_3$. Therefore, there exists $t_2\in (t^*, t^{*}+T']$ so that $x_1(t_2) > {m_3}$.

Now, we define $\tilde{t} = \inf_{t > t^*}\{x_1(t) > m_3\}$. This means that, $x_1(t) < {m_3}$ for $t\in(t^*, \tilde{t})$ and $x_1(\tilde{t}) = {m_3}$ by the continuity of $x_1(t)$. Then, we choose $p\in {\mathbb{Z}_+}$ so that $p\le{n_2}+{n_3}$ and $t^{*}+pT\ge\tilde{t}$, and suppose $t\in (t^{*}+(p-1)T, t^{*}+pT]$. By inequality (3.27), we get

Since ${\eta_1} < 0$ and $p\le {n_2}+{n_3}$.

Let

Thus, $x_1(t)\ge\bar{m}_2$ for $t\in (t^*, \tilde{t})$. Similarly, we use

$\tilde{t}$ instead of $t^*$. Then, we will obtain $x_1(t)\ge\bar{m}_2$ for all sufficiently large $t$.

Case 2: $t^*\ne nT$ for all $n\in {\mathbb{Z}_+}$. That is $x_1(t)\ge {m_3}$ for $t\in [t_1, t^*)$ and $x_1(t^*) = {m_3}$. We assume that $t^* \in ({\bar{n}_1}T, ({\bar{n}_1}+1)T)$, for some ${\bar{n}_1}\in {\mathbb{Z}_+}$. We can consider this in two subcases.

Case 2.1: $x_1(t)\le{m_3}$ for all $t\in (t^*, ({\bar{n}_1}+1)T]$. Suppose that there is ${t'_2}\in [({\bar{n}_1}+1)T, ({\bar{n}_1}+1)T+T']$ so that $x_1(t'_2) > {m_3}$. Otherwise, considering Eq (3.19) with

For $t\in (nT, (n+1)T]$, ${\bar{n}_1}+1\le n\le {\bar{n}_1}+1+{n_2}+{n_3}$, we obtain

In a similar way to Case 1, for ${n_2}T\le t-{t^*}$, we get

Thus,

Since ${n_2}T\le({\bar{n}_1}+1+{n_2})T-{t^*}$, we get

Then,

It is in contradiction to the definition of $m_3$. Thus, there exists $t'_2 \in [({\bar{n}_1}+1)T, ({\bar{n}_1}+1)T+T']$ so that $x_1(t'_2) > {m_3}$.

Now, we define $\bar{t} = \inf_{t > t^*}\{x_1(t) > m_3\}$. Thus, $x_1(t)\le {m_3}$ for $t\in [t^*, \bar{t})$, and $x_1(\bar{t}) = {m_3}$. We choose $p'\in {Z_+}$ such that $p'\le {n_2}+{n_3}+1$ and suppose $t\in ({\bar{n}_1}T+(p'-1)T, {\bar{n}_1}T+p'T]$. From inequality (3.27), we get

Since ${\eta_1} < 0$ and $t-{t^*}\le p'T$. Then,

Let

Therefore, $x_1(t)\ge {m_2}$ for $t\in (t^*, \bar{t})$. We do the same way by using $\bar{t}$ instead of $t^*$. Then, we will obtain $x_1(t)\ge {m_2}$ for all sufficiently large $t$.

Case 2.2: There exists $t''\in (t^*, ({\bar{n}_1}+1)T]$ so that $x_1(t'') > {m_3}$. Define $\underline{t} = \inf_{t > t^*}\{x_1(t) > m_3\}$. Hence, $x_1(t) < {m_3}$ for $t\in [t^*, \underline{t})$, and $x_1(\underline{t}) = {m_3}$. For $t\in [t^*, \underline{t})$, inequality (3.27) holds, we get

since $t < {\bar{n}_1}T+T < {t^*}+T$. For $t > \underline{t}$, we can do the same way since $x_1(\underline{t})\ge {m_3}$. Since ${m_2} < {\bar{m}_2} < {m_3}$, we can conclude that $x_1(t)\ge {m_2}$ for $t\ge {t_1}$. The proof is completed.

3.3.   The positive periodic solution

Now, we carry out the conditions to guarantee the positive periodic solution of the system (2.1)–(2.4) near the periodic solution $(0, \tilde{x}_2)$. For more convenience, we change the variables, and then the new system is shown as follows:

for $t \neq nT$ with

Let

According to Lakmeche and Arini ^[13],

and

where

Now, we can compute

where $\tau_0$ is the root of $d'_0 = 0$. Note that $d'_0 > 0$ if $T < T_{max}$ and $d'_0 < 0$ if $T > T_{max}$.

Note that $a'_0 > 0$ if $T > T_{min}$.

Note that $G^* < 0$ if

and $H^* > 0$ if

Thus, ${G^*}{H^*} < 0$, and by Lakmeche and Arini ^[13], we obtain the following theorem.

Theorem 3.9. If all conditions $(3.2)$, $(3.6)$, $(3.11)$, $(3.12)$, $(3.14)$, $(3.32)$, $(3.33)$ and $T > T_{max} > T_{min}$ hold, then the system $(3.28)–(3.31)$ has a positive periodic solution which is supercritical.

4.   Numerical results and interpretation

In this section, the numerical results of the system (2.1)–(2.4) are carried out by using ode15 package in MATLAB to confirm the analysis of solutions.

The numerical simulations of the (2.1)–(2.4) are computed by using the parameters and initial conditions given in Table 1.

The remaining parameters $\gamma = 0.9, \omega = 0.1, T = 0.1$ are used to simulate the solutions as shown in Figures 1. All parameters are satisfied the conditions in Theorem 3.6. The trend of solution is close to a limit cycle $(0, \tilde{x}_2)$ as proved.

The initial situation starts with 10 units in both of population density of small fishes $x_1(0)$ and short mackerel $x_2(0)$. After that, the population of small fishes continues to decrease and tends to zero due to the short period and high quantity of fisheries. However, the population of shot mackerel decreases in the beginning then it tends to a closed orbit between $0.2593$–$0.2881$ because of the low rate of toxic death and high capacity of hunting with the same period of fisheries.

The computer simulations of the system (2.1)–(2.4) with setting parameters in Table 1 and $\gamma = 0.5, \omega = 0.2, T = 0.5$ are shown in Figure 2. For this case, we have that all parameters are satisfied the conditions in Theorem 3.8. The solution is permanent as proved.

The initial situation starts with 10 units in both of population density of small fishes $x_1(0)$ and short mackerel $x_2(0)$. After that, both small fishes and shot mackerel decrease until tending to a range of $0.0987$–$0.1974$ and $0.2347$–$0.2934$, respectively. The fisherman catches enough of them for the economy and they remain alive in the system. The long period of fisheries, low-frequency fisheries, and the appropriate number of fishing are the most of factors in persistence.

The computer simulations of the system (2.1)–(2.4) with the parametric values $\gamma = 0.7, \omega = 0.3, T = 1.5$ and the remaining parameters as in Table 1 are shown in Figure 3. All parameters in this case are satisfied the conditions in Theorem 3.9. The system has a positive periodic solution as proved.

The initial situation start with $10$ units both of population density of small fishes $x_1(0)$ and short mackerel $x_2(0)$. Then both small fishes and shot mackerel decrease until tending to oscillatory in a narrow range of 0.0592–0.1974 and 0.2054–0.2934, respectively. The decreased population of small fishes and shot mackerel occurs when the period of fisheries starts while the population of them back hits the peak when fisherman do not allow to fishing like a periodic behavior. The suitable amount and period of fisheries like this case to prevent the extinction of small fishes and short mackerel and to prevent the damage of economic are expected situation.

The computer simulations of the system (2.1)–(2.4) with the parametric values in Table 1 and $\gamma = 0.5, \omega = 0.2$ which focused on changes in the period of impulsive fisheries $(T)$ and the negative effect of fisheries on the density of short mackerel population $(\omega)$ with $\gamma = 0.1, T = 0.2$ are shown in Figure 4 and Figure 5 respectively. The results indicated that the densities of the short mackerel were in periodic fashions. Moreover, the different values of $T$ and $\omega$ provided the different highest densities of $x_2(t)$.

5.   Discussion and conclusions

We propose the modification mathematical model to forecast the population density dynamic of Indo-Pacific mackerel (Rastrelliger brachysoma) or short mackerel. The proposed model is utilized to control the population of short mackerel by considering the decrease population affected by catching, toxic, and natural death. The suitable period $T$ and the quantity affecting the decreasing rate of small fishes population $\gamma$ and the decreasing rate of short mackerel population $\omega$ are essential things to maintain short mackerel population and small fish population without extinction. The numerical results show the periodic behavior of the density population. Moreover, there are enough short mackerel for fishermen and humans for a long time.

Acknowledgments

This research project was financially supported by Thailand Science Research and Innovation (TSRI) 2021.

Conflict of interest

All authors declare no conflicts of interest in this paper.