Optimal harvesting of a competitive n-species stochastic model with delayed diffusions

1.   Introduction

In mathematical modelling, the term diffusion is used to describe the motion of species from one region to another. Influenced by various natural factors, such as geographic, hydrological or climatic conditions and human activities, migrations occur between patches, which affects the population dynamics, for example the persistence and extinction of species ^{[1,2,3,4,5,6,7,8]}. The growth of species population is also affected by competition caused by disputing food, resources, territories and spouses, including intraspecific and interspecific competitions among populations. To see the effects of the diffusion and competition on population dynamics, we propose the following mathematical model with $n$ species, for $\;i, j = 1, 2, \ldots, n,$

where $x_{i}(t)$ is the population size at time $t$ of the $i$th species, positive constants $r_{i}$, $a_{ii}$ are the growth rate and the interspecific competition rate of the $i$th species respectively, $a_{ij}>0\;(j\neq i)$ is the competition rate between species $i$ and $j$, $D_{ij}\geq 0$ is the diffusion coefficient from species $j$ to species $i$, $\alpha_{ij}\geq 0$ indicates the diffusion boundary condition.

Recently, time delays have been widely used in biological and ecological models in order to get more realistic mathematical models, for example ^{[9,10,11,12,13,14,15,16]}. In this paper, we also consider the time delay, which is accounted for the diffusion. For example, birds cannot migrate immediately after they were born, so the time delay here is the time it takes for them to learn to fly before they can migrate, and death can also occur in the process. Then, from (1) we have the model with time delays as follows

where $\tau_{ij}\geq 0$ is the time delay and $d_{j}$ is the death rate of the $j$th species. Let $\tau = \max_{i, j = 1, \ldots, n}\{\tau_{ij}\}$ and $C([-\tau, 0];R_{+}^{n})$ denote the family of all bounded and continuous functions from $[-\tau, 0]$ to $R_{+}^{n}$. We assume model (2) is subject to the following initial condition

Reference ^[17] suggests that the growth rate of organisms is generally affected by environmental fluctuations accounted for the disturbance of ecological environment in nature, consequently parameters in biologic models will exhibit random perturbations ^[18]. Thus, the deterministic models, like (2) are not applicable to capture the essential characters. In the past years, researchers have suggested the use of white noises to capture the main characters of these stochastic fluctuations, see ^{[18,19,20,21,22,23,24,25,26,27]} for example. Denote by $\{B_{i}(t)\}_{t\geq 0}, (i = 1, 2, \ldots, n)$ the independent standard Brownian motions defined on a complete probability space $(\Omega, \{F_{t}\}_{t\in R_{+}}, P)$ with $\sigma_{i}^{2}$ represents the intensity of the environment noises. Then, the growth rate subject to random perturbation can be described by

with which the model (2) reads

We further consider the optimal harvesting problem of model (4). The research on the optimal harvesting of the population is of great significance to the utilization and development of resources, and can also help mankind to get the optimal strategy of harvesting in order to obtain the most long-term benefits ^{[28,29,30,31,32,33,34,35]}. Then, we reach the following model accounted for harvesting:

where $h_{i}\geq 0$ denotes the harvesting effort of the species $i$.

In the rest of the paper, we will devote ourselves to explore the dynamics and the optimal harvesting strategy of model (5). More precisely, in Section 2, we establish necessary conditions for persistence of species in mean and extinction of the species. In Section 3, we investigate conditions of stability, and prove asymptotic stability in distribution of the model, namely, there is a unique probability measure $\rho(\cdot)$ such that for each $\phi \in C([-\tau, 0];R_{+}^{n})$, the transition probability $p(t, \phi, \cdot)$ of $x(t)$ converges weekly to $\rho(\cdot)$ when $t\rightarrow\infty$. In Section 4, by the use of the Hessian matrix and theorems of optimal harvesting due to ^[36], we investigate the optimal harvesting effort and gain the maximum of expectation of sustainable yield (ESY). In Section 5, we numerically illustrate our theoretical results obtained in previous sections, and then conclude our study in Section 6.

2.   Persistence and extinction

For the convenience of the following discussion, we define some notations as follows

and assume that $\sum_{j = 1, j\neq i}^{n}a_{ij}\geq \sum_{j = 1, j\neq i}^{n}D_{ij}e^{-d_{j}\tau_{ij}}$ holds in the rest of the paper.

Following the same argument as in ^[37], we can prove the existence of the positive solution.

Lemma 2.1. Given initial value (3), model (5) admits a unique global positive solution $x(t) = \left(x_{1}(t), \ldots, x_{n}(t)\right)^{T}$ almost surely. Furthermore, for each $p>1$, there exists a positive constant $K = K(p)$ such that

To show our main result of this section, we consider the following auxiliary equations

with initial value

By [38,Stochastic Comparison Theorem], we know that for $t\geq -\tau$,

Remark 1. It is easy to see from ^[39] that the explicit solution of (7) is

Similar calculation gives

Then, by using ^[40], we obtain the following.

Lemma 2.2. Let $b_{i}>0$. Then, from (7) we have

Based on Lemma 2, we assume:

Assumption 2.1. $b_{i}>0, \;c_{i}>0, \;i = 1, 2, \ldots, n.$

Remark 2. A result due to Golpalsamy ^[10] and Assumption 2.1 imply that there exists a unique positive solution $(\det(A_{1})/\det (A), \ldots, \det(A_{n})/\det (A))^{T}$ for the following system

in which

and $A_{i}$ is the matrix given by using the $(b_{1}, b_{2}, \ldots, b_{n})^{T}$ to replace the $i$th column of matrix $A$.

Now we are in the position to show our main results.

Theorem 2.1. All species in system (5) are persistent in mean $a.s.$, $i.e.$,

when Assumption 2.1 is satisfied.

Proof. Let $b_{i}>0$, according to (12) that for $i, j = 1, 2, \ldots, n, \;j\neq i$, one has

which together with (9) yields

By using Itô's formula to (5), one can see that

According to (16) together with the property of Brownian motion, we obtain

We next to show that

In view of (9) and (12), we have

Therefore we obtain

From (15) and (12), we get

By using $\lim_{t\rightarrow+\infty}t^{-1}B_{i}(t) = 0$ together with what we have just obtained, yields that for any given $\varepsilon>0$, there exists a $T = T(\omega)$ thus for $t\geq T, \;i, j = 1, 2\ldots, n, \;i\neq j,$

Applying these inequalities to (11), we have

in which $M_{ij}>0$ is a constant. Note that $c_{i} = b_{i}-\sum_{j = 1, j\neq i}^{n}\frac{a_{ij}}{q_{ji}}b_{j}>0$, thereby for large enough $t$, one has that

Hence for sufficiently large $t$, we obtain

Rearranging this inequality shows that

Since $t$ is large enough and $\varepsilon$ is arbitrary, we get (14). This completes the proof of Theorem 2.1.

Corollary 2.1. If there is a $b_{i}<0$, then according to (17), one has $\limsup_{t\rightarrow+\infty}t^{-1}\ln x_{i}(t)\leq b_{i}<0, \;a.s.$ It is to say $\lim_{t\rightarrow+\infty}x_{i}(t) = 0, \;a.s.$, which means that the $i$th species in system (5) will die out.

3.   Stability in distribution

In this section, we study the stability of the model. To this end, we suppose the following holds:

Assumption 3.1. $a_{ii}+\sum_{j = 1, j\neq i}^{n}D_{ij}\alpha_{ij}\geq \sum_{j = 1, j\neq i}^{n}a_{ji}+\sum_{j = 1, j\neq i}^{n}D_{ji}e^{-d_{i}\tau_{ji}}, \;i = 1, 2, \ldots, n.$

Then, we can prove the following.

Theorem 3.1. The system (5) is asymptotically stable in distribution if Assumption 3.1 holds.

Proof. Given two initial values $\phi(\theta), \psi(\theta) \in C\left([-\tau, 0];R_{+}^{n}\right)$ of model (5), the corresponding solutions are $x^{\phi}(t) = (x_{1}^{\phi_{1}}(t), \ldots, x_{n}^{\phi_{n}}(t))^{T}$ and $x^{\psi}(t) = (x_{1}^{\psi_{1}}(t), \ldots, x_{n}^{\psi_{n}}(t))^{T}$ respectively. Let

Applying Itô's formula yields

Therefore

Together with $\mathbb{E}(V(t))\geq 0$, one has

Hence we have $\mathbb{E}\bigg|x_{i}^{\phi_{i}}(s)-x_{i}^{\psi_{i}}(s)\bigg|\in L^{1}[0, \infty), \;i = 1, 2, \ldots, n.$ At the same time, by using (5) we obtain that

That is to say $\mathbb{E}(x_{i}(t))$ is continuously differentiable with respect of $t$. Computing by (5) leads to

in which $K>0$ is a constant. It implies that $\mathbb{E}(x_{i}(t))$ is uniformly continuous. Using ^[41], we get

Denote $p(t, \phi, {\rm d} y)$ as the transition probability density of the process $x(t)$ and $P(t, \phi, A)$ represents the probability of event $x(t)\in A$. By (6) and [42,Chebyshev's inequality], we can obtain that the family of $p(t, \phi, {\rm d} y)$ is tight. Now define $\Gamma\left(C\left([-\tau, 0];R^n_+\right)\right)$ as the probability measures on $C\left([-\tau, 0];R^n_+\right)$. For arbitrary two measures $P_1, P_2\in \Gamma$, we define the metric

where

Since $\{p(t, \phi, {\rm d} y)\}$ is tight, then according to (19) we know that for any $\varepsilon>0$, there is a $T>0$ satisfies that for $t\geq T, s>0$,

Therefore $\{p(t, \xi, \cdot)\}$ is Cauchy in $\Gamma$ with metric $d_\mathbb{L}$, in which $\xi\in C\left([-\tau, 0];R^n_+\right)$ is arbitrary given. Hence there exists a unique $\kappa(\cdot)\in \Gamma(C([-\tau, 0];R^n_+))$ such that $\lim\limits_{t\rightarrow\infty}d_\mathbb{L}(p(t, \xi, \cdot), \kappa(\cdot)) = 0$. At the same time, it follows from (19) that

Consequently,

This completes the proof of Theorem 3.1.

4.   Optimal harvesting

In this section, we consider the optimal harvesting problem of system (5). Our purpose is to find the optimal harvesting effort $H^{*} = (h_{1}^{*}, \ldots, h_{n}^{*})$ such that:

(ⅰ) $Y(H) = \lim_{t\rightarrow+\infty}\sum_{i = 1}^{n}\mathbb{E}(h_{i}x_{i}(t))$ is maximum;

(ⅱ) Every $x_{i}(i = 1, 2, \ldots, n)$ is persistent in the mean.

Before we give our main results, we define

in which $G = (r_{1}-0.5\sigma_{1}^{2}, r_{2}-0.5\sigma_{2}^{2}, \ldots, r_{n}-0.5\sigma_{n}^{2})^{T}$ and $I$ is the unit matrix, and make an assumption:

Assumption 4.1. $A^{-1}+(A^{-1})^{T}$ is positive definite,

Theorem 4.1. Suppose Assumptions 3.1 and 4.1 hold, and If these following inequalities

are satisfied. Then, for system (5) the optimal harvesting effort is

and the maximum of ESY is

Proof. Denote $W = \{H = (h_{1}, \ldots, h_{n})^{T}\in R^{n}\mid b_{i}>0, \;c_{i}>0, \;h_{i}>0, \;i = 1, \ldots, n\}$. Easily we can see that for any $H\in W$, (14) is satisfied. Note that $\Theta\in W$, then $W$ is not empty. According to (14), we have that for every $H\in W$,

Applying Theorem 4.1, there is a unique invariant measure $\rho(\cdot)$ for model (5). By [43,Corollary 3.4.3], we obtain that $\rho(\cdot)$ is strong mixing. Meanwhile, it is ergodic according to [43,Theorem 3.2.6]. It means

Let $\mu(x)$ represent the stationary probability density of system (5), then we have

Since the invariant measure of model (9) is unique, one has

In other words,

Assume that $\Theta = (\theta_{1}, \theta_{2}, \ldots, \theta_{n})^{T}$ is the solution of the following equation

Thus, $\Theta = [A(A^{-1})^\mathrm{T}+I]^{-1}G$. By using of the Hessian matrix (see ^[44,45]),

is negative defined, then $\Theta$ is the unique extreme point of $Y(H)$. That is to say, if $\Theta\in W$ and under the condition of (21), the optimal harvesting effort is $H^* = \Theta$ and $Y^*$ is the maximum value of ESY. This completes the proof of Theorem 4.1.

5.   Numerical simulations

To see our analytical results more clearly, we shall give some numerical simulations in this section. Without loss of generality, we consider the following system

which is the case when $n = 2$ in (5), with initial value

where $r_{i}>0, \;a_{ij}>0, \;\tau_{i}\geq 0, \;i, j = 1, 2$.

Firstly, we discuss the persistence in mean of $x_{1}$ and $x_{2}$. For that, we take the parameter values as follows:

The initial values are $x_{1}(\theta) = 0.5+0.01\sin\theta$, $x_{2}(\theta) = 0.5+0.02\sin\theta$, $\theta\in[-\tau, 0]$. Simple calculations show that $b_{1} = 0.7988>0, \;b_{2} = 0.7488>0, \;c_{1} = 0.6662>0, \;c_{2} = 0.6556>0$ implying Assumption 2.1 is satisfied. Then by Theorem 2.1, we can obtain that in (29)

Applying the Milstein numerical method in ^[47], we then obtained the numerical solution of system (29), see Figure 1. It shows that $x_{1}$ and $x_{2}$ respectively asymptotical approach to 0.2268 and 0.5964 time averagely. And this agrees well with our results in Theorem 2.1. Then we research the distributions of $x_{1}$ and $x_{2}$ under the same conditions. Obviously, we have $a_{11}+D_{12}\alpha_{12}\geq a_{21}+D_{21}e^{-d_{1}\tau_{21}}, \;a_{22}+D_{21}\alpha_{21}\geq a_{12}+D_{12}e^{-d_{2}\tau_{12}}$, it is to say Assumption 3.1 is satisfied. Thus by Theorem 3.1, system (29) is asymptotically stable in distribution as suggested by Figure 2.

Lastly, we consider the optimal harvesting strategy of system (29). It is easy to see that the Assumption 2.1 and Assumption 3.1 are satisfied. Furthermore, we have

in which $I = \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right)$. Since condition (21) is satisfied, by Theorem 4.1, the optimal harvesting effort is

on the other hand, the maximum of ESY is

By using the Monte Carlo method (see ^[48]) and the parameters in Table 1, we can obtain Figure 3, showing our results in Theorem 4.1.

Next, we consider a case of three species.

We use the following parameter values:

The initial values are $x_{1}(\theta) = 0.5+0.01\sin\theta$, $x_{2}(\theta) = 0.5+0.02\sin\theta$, $x_{3}(\theta) = 0.5+0.001\sin\theta$, $\theta\in[-\tau, 0]$. Easily we get that $b_{1} = 0.1.5536>0, \;b_{2} = 0.7881>0, \;b_{3} = 0.2681>0, \;c_{1} = 1.4502>0, \;c_{2} = 0.0552>0, \;c_{3} = 0.0229>0.$ Thus Assumption 2.1 is hold. By Theorem 2.1, we have for (30)

The numerical results of Theorem 2.1 when $n = 3$ are shown in Figure 4.

The stable distribution for $n = 3$ are shown in Figure 5.

To numerical illustrate the optimal harvesting effort of (30), we set

which yield $H^{*} = \Theta = (1.1052, 0.5537, 0.1663)^{T}$, and the maximum of ESY is $Y^{*} = 0.2263$, see Figure-6.

6.   Conclusions and discussions

In this paper, a stochastic n-species competitive model with delayed diffusions and harvesting has been considered. We studied the persistence in mean of every population, which is biologically significant because it shows that all populations can coexist in the community. Since the model (5) does not have a positive equilibrium point and its solution can not approach a positive value, we considered its asymptotically stable distribution. By using ergodic method, we obtained the optimal harvesting policy and the maximum harvesting yield of system (5). We have also done some numerical simulations of the situations for $n = 2$ and $n = 3$ in model (5) to illustrate our theoretical results as it is very useful whether in terms of mathematics or biology to visualize our conclusions.

Our studies showed some interesting results

(a) Both environmental disturbance and diffused time delay can effect the persistence and optimal harvesting effort of system (5)..

(b) Environmental noises have no effect on asymptotic stability in distribution of system (5), but the time delays have.

There are other meaningful aspects that can be studied further since our paper only consider the effects of white noises on population growth rate. In future, for example, we can consider the situation when white noises also have influences over harvesting (see ^[45]) and non-autonomous system (see ^[46]); the time delay will also be reflected in competition (see ^[49]). Furthermore, we can consider something more complex models such as the ones with regime-switching (see ^[50,51]) or Lévy jumps (see ^[14,42]).

Acknowledgments

This work was supported by the Research Fund for the Taishan Scholar Project of Shandong Province of China, and the SDUST Research Fund (2014TDJH102).

Conflict of interest

The authors declare that there is no conflict of interest regarding the publication of this paper.