On the interplay of harvesting and various diffusion strategies for spatially heterogeneous populations

TL;DR

This paper investigates how different harvesting strategies affect the competition outcomes of two species in a spatially heterogeneous environment, revealing conditions for coexistence and exclusion.

Contribution

It introduces a model analyzing the combined effects of harvesting and diffusion strategies on species competition in heterogeneous habitats.

Findings

Harvesting can ensure survival of non-harvested species.

Coexistence is possible within certain harvesting rate ranges.

Specific estimates relate harvesting levels to competition outcomes.

Abstract

The paper explores the influence of harvesting (or culling) on the outcome of the competition of two species in a spatially heterogeneous environment. The harvesting effort is assumed to be proportional to the space-dependent intrinsic growth rate. The differences between the two populations are the diffusion strategy and the harvesting intensity. In the absence of harvesting, competing populations may either coexist, or one of them may bring the other to extinction. If the latter is the case, introduction of any level of harvesting to the successful species guarantees survival to its non-harvested competitor. In the former case, there is a strip of "close enough" to each other harvesting rates leading to preservation of the original coexistence. Some estimates are obtained for the relation of the harvesting levels providing either coexistence or competitive exclusion.

Equations

$$\frac{\partial u(t,x)}{\partial t}=\nabla \cdot \left[a(x)\nabla \left(\frac{u(t,x)}{P(x)}\right)\right]+r(x)u(t,x)\left(1-\frac{u(t,x)+v(t,x)}{K(x)}\right)-\alpha r(x)u(t,x),$$

$$\frac{\partial v(t,x)}{\partial t}=\nabla \cdot \left[b(x)\nabla \left(\frac{v(t,x)}{Q(x)}\right)\right]+r(x)v(t,x)\left(1-\frac{u(t,x)+v(t,x)}{K(x)}\right)-\beta r(x)v(t,x),$$

$$t>0,x\in \Omega ,$$

$$\frac{\partial u}{\partial n}-\frac{u}{P}\frac{\partial P}{\partial n}=\frac{\partial v}{\partial n}-\frac{v}{Q}\frac{\partial Q}{\partial n}=0,x\in \partial \Omega ,$$

$$u(0,x)=u_0(x),v(0,x)=v_0(x),x\in \Omega .$$

$$\frac{\partial u(t,x)}{\partial t}=\nabla \cdot \left[a(x)\nabla \left(\frac{u(t,x)}{P(x)}\right)\right]+r_{1}r(x)u(t,x)\left(1-\frac{u(t,x)+v(t,x)}{K_1(x)}\right),$$

$$\frac{\partial v(t,x)}{\partial t}=\nabla \cdot \left[b(x)\nabla \left(\frac{v(t,x)}{Q(x)}\right)\right]+r_{2}r(x)v(t,x)\left(1-\frac{u(t,x)+v(t,x)}{K_2(x)}\right),$$

$$t>0,x\in \Omega ,$$

$$\frac{\partial u}{\partial n}-\frac{u}{P}\frac{\partial P}{\partial n}=\frac{\partial v}{\partial n}-\frac{v}{Q}\frac{\partial Q}{\partial n}=0,x\in \partial \Omega ,$$

$$u(0,x)=u_0(x),v(0,x)=v_0(x),x\in \Omega .$$

$$\nabla \cdot \left[a(x)\nabla \left(\frac{u^{\ast }(x)}{P(x)}\right)\right]+r_{1}r(x)u^{\ast }(x)\left(1-\frac{u^{\ast }(x)}{K_1(x)}\right)=0,x\in \Omega ,\frac{\partial (\frac{u^{\ast }}{P})}{\partial n}=0,x\in \partial \Omega$$

$$\nabla \cdot \left[b(x)\nabla \left(\frac{v^{\ast }(x)}{Q(x)}\right)\right]+r_{2}r(x)v^{\ast }(x)\left(1-\frac{v^{\ast }(x)}{K_2(x)}\right)=0,x\in \Omega ,\frac{\partial (\frac{v^{\ast }}{Q})}{\partial n}=0,x\in \partial \Omega ,$$

$${\alpha }_1\ge {\alpha }^{\ast }=1-\frac{{\int }_{\Omega }r{v}_{\beta }^{\ast }dx}{{\int }_{\Omega }rKdx},$$

$$\nabla \cdot \left[a(x)\nabla \left(\frac{K(x)}{P(x)}\right)\right]\not\equiv 0,\nabla \cdot \left[b(x)\nabla \left(\frac{K(x)}{Q(x)}\right)\right]\not\equiv 0,$$

$$K(x)=\gamma P(x)+\delta Q(x),x\in \Omega .$$

$${\alpha }_2>{\alpha }^{\ast }=1-\frac{{\int }_{\Omega }Pr{v}_{\beta }^{\ast }/Kdx}{{\int }_{\Omega }rPdx},$$

$$\mathrm{SY}:={\int }_{\Omega }\alpha r(x)u_s(x)dx.$$

$$\mathrm{SY}:={\int }_{\Omega }\alpha r(x)u_s(x)dx+{\int }_{\Omega }\beta r(x)v_s(x)dx.$$

$$\mathrm{MSY}:={\int }_{\Omega }\frac{1}{4}r(x)K(x)dx.$$

$$\underset{\Omega }{\int }r(x)K_2(x)dx>\underset{\Omega }{\int }r(x)v^{\ast }(x)dx.$$

$$\frac{\partial u(t,x)}{\partial t}=\nabla \cdot \left[a(x)\nabla \left(\frac{u(t,x)}{P(x)}\right)\right]+r_{1}r(x)u(t,x)\left(1-\frac{u(t,x)+v(t,x)}{cK(x)}\right),$$

$$\frac{\partial v(t,x)}{\partial t}=\nabla \cdot \left[b(x)\nabla \left(\frac{v(t,x)}{Q(x)}\right)\right]+r_{2}r(x)v(t,x)\left(1-\frac{u(t,x)+v(t,x)}{K(x)}\right),$$

$$t>0,x\in \Omega ,$$

$$\frac{\partial u}{\partial n}-\frac{u}{P}\frac{\partial P}{\partial n}=\frac{\partial v}{\partial n}-\frac{v}{Q}\frac{\partial Q}{\partial n}=0,x\in \partial \Omega ,$$

$$u(0,x)=u_0(x),v(0,x)=v_0(x),x\in \Omega .$$

$$\frac{\partial u(t,x)}{\partial t}=\nabla \cdot \left[a(x)\nabla \left(\frac{u(t,x)}{P(x)}\right)\right]+r_{1}r(x)u(t,x)\left(1-\frac{v^{\ast }(x)}{cK(x)}\right),$$

$$\nabla \cdot \left[a\nabla \left(\frac{\psi }{P}\right)\right]+r_{1}r\psi \left(1-\frac{v^{\ast }}{cK(x)}\right)=\sigma \psi ,x\in \Omega ,\frac{\partial (\psi /P)}{\partial n}=0,x\in \partial \Omega .$$

$${\sigma }_1=\underset{\psi \ne 0,\psi \in W^{1,2}}{\sup }\left[-\underset{\Omega }{\int }a|\nabla (\psi /P){|}^{2}dx+\underset{\Omega }{\int }{r}_{1}r\frac{{\psi }^2}{P}\left(1-\frac{v^{\ast }}{cK}\right)dx\right]/\underset{\Omega }{\int }\frac{{\psi }^2}{P}dx.$$

$$-\underset{\Omega }{\int }a|\nabla (\Psi /P){|}^{2}dx+\underset{\Omega }{\int }{r}_{1}r\frac{{\Psi }^2}{P}\left(1-\frac{v^{\ast }}{K}\right)dx-{\sigma }_1\underset{\Omega }{\int }\frac{{\Psi }^2}{P}dx=0.$$

$$0<\underset{\Omega }{\int }\frac{{\Psi }^2}{P}{r}_{1}r\frac{v^{\ast }}{K}dx=-\underset{\Omega }{\int }a|\nabla (\Psi /P){|}^{2}dx+\underset{\Omega }{\int }{r}_{1}r\frac{{\Psi }^2}{P}dx-{\sigma }_1\underset{\Omega }{\int }\frac{{\Psi }^2}{P}dx.$$

$$c>c_1:=\underset{\Omega }{\int }\frac{{\Psi }^2}{P}{r}_{1}r\frac{v^{\ast }}{K}dx{\left[-\underset{\Omega }{\int }a|\nabla (\Psi /P){|}^{2}dx+\underset{\Omega }{\int }{r}_{1}r(x)\frac{{\Psi }^2}{P}dx\right]}^{-1},$$

$${\sigma }_1\ge \frac{r_1}{M}\underset{\Omega }{\int }r(x)K(x)\left(1-\frac{v^{\ast }(x)}{cK(x)}\right)dx=\frac{r_1}{M}\left[\underset{\Omega }{\int }r(x)K(x)dx-\frac{1}{c}\underset{\Omega }{\int }r(x)v^{\ast }(x)dx\right]>0$$

$$c^{\ast }=\frac{\underset{\Omega }{\int }r(x)v^{\ast }(x)dx}{\underset{\Omega }{\int }r(x)K(x)dx}<1.$$

$$0<m:=\underset{x\in \overline{\Omega }}{\min }K(x),M:=\underset{x\in \overline{\Omega }}{\max }K(x).$$

$$0<b:=\underset{x\in \overline{\Omega }}{\min }{v}^{\ast }(x),B:=\underset{x\in \overline{\Omega }}{\max }{v}^{\ast }(x).$$

$$r(x)\overline{u}\left(1-\frac{\overline{u}+v_s(x)}{cK(x)}\right)\le r(x)\overline{u}\left(1-\frac{\overline{u}}{cK(x)}\right)=r(x)\overline{u}\left(1-\frac{b}{2cM}\right)<0$$

Full text

††journal: Journal of Theoretical Biology

1 Introduction

Harvesting is quite a common cause for extinction of exploited populations. However, the study of one harvested species [31], generally, does not adequately describe the population dynamics in multi-species environment, as it ignores interactions, such as competition, predation or mutualism. If, in the absence of harvesting, two populations coexist, will intervention bring the harvested population to extinction? If originally the resident species was immune to a possible invasion by a competitor, will its being harvested open the habitat to a possible invasion? If two species coexist and compete for the same resources, and both are harvested, which relations between the harvesting rates will preserve coexistence?

Spatial structure, as well as competition for resources with other species, can significantly influence the outcome of the competition. Structural or spatial heterogeneity can be beneficial for exploited species, for example, providing refuge areas, or harmful, for instance, creating age-specific over-exploitation, or crucial disruptions in the spread over habitat, or vulnerability under seasonal migrations.

Studies on harvesting were usually focused on establishing the maximum sustainable yield, investigating at the same time optimality and sustainability of harvesting policies. Usually, they were limited to one population, see, for example, [1, 3, 5, 19, 20, 22, 31]. Special caution is required for the choice of harvesting strategies when the Allee effect is incorporated in population dynamics [1, 20]. More interesting scenarios arise when harvesting is applied to one or more of several interacting populations [14, 15, 18, 19, 27, 32, 33]. In certain cases, a generalist predator was introduced to control invasion and was considered as a source of harvesting [29].

While either homogeneous environments (modeled by a system of ordinary differential equations) or non-selective harvesting can result in competitive exclusion of all but one species [14, 19], consideration of structured environment and different types of harvesting applied to distinct species, can lead to coexistence as well, see, for example, [33]. In certain cases, populations with and without diffusion [29], or in the case of small diffusion coefficients [16] were considered. In the non-harvested case for two competing species, common scenarios are either competitive exclusion or coexistence [13]. If competitors differ by one parameter only, for example, dispersal strategies or carrying capacities, this led to a competitive exclusion to favor a more evolutionarily advantageous property or strategy. Still, for harvested species, the interplay of two or more factors (different spatial distributions, model parameters and harvesting rates) can eventually lead to anyone of the three different scenarios, two of competitive exclusion, and coexistence. In principle, if both species are over-exploited, the fourth scenario of the total extinction is possible.

In the present paper, we investigate a competition of two species choosing different dispersal strategies and being harvested with various harvesting efforts. Evolutionary advantage or disadvantage of fast and slow, random or directed movements are closely connected to heterogeneity of the environment. For spatially homogeneous environments and regular dispersal, a slower diffuser [17] is a competition winner, if no other differences between competitors exist. If the environment is spatially heterogeneous then the species using knowledge on the carrying capacity of the environment in its dispersal strategy can bring its regularly diffusing competitor to extinction [7, 10, 12, 23]. The ideal free distribution describes the spatial structure when any movement would decrease the population fitness. For spatially heterogeneous but stationary environments, this may correspond to the case of the population density coinciding with the carrying capacity of the environment. The idea that an evolutionarily stable strategy corresponds to the case, when the ideal free distribution is a limit solution was justified, for example, in [2, 9]. However, the logistic equation with a regular diffusion (the Fisher equation) does not have an ideal free distribution as a solution, which is extensively discussed in the recent paper [16]. Moreover, faster diffusion leads to a disadvantage in a competition [17]. Several developments, such as introducing advection or considering the modified dispersal term, were aimed to remedy this situation [2, 5, 9, 10, 11, 12, 25]. An interesting investigation in this area in the discrete case was recently reported in [13]. In the continuous case for two (or more) species, cooperation in resources consumption is possible when the densities of the two populations complement each other forming an ideal free pair [2, 6]. For such a pair, the total density exactly matches the carrying capacity, and the diffusion strategies are aligned with these distributions. This can describe the situation when accessibility of resources is spatially different for the two species, i.e. specialization by resource consumption.

For a given type of diffusion, the average population rates depend on the two important parameters, assumed to be space-dependent in heterogeneous environments: the intrinsic growth rate $r$ and the carrying capacity $K$, see [16, 19] and references therein. If a population is exploited, different harvesting strategies can be introduced, such as constant level, proportional (constant effort) and threshold harvesting. For the latter case, the deduction is constant for high enough population rates but either becomes close to proportional [31] or is canceled [19] when the population level is below a prescribed threshold.

Most results on the outcomes of the competition stipulated by various diffusion strategies and other parameters, are based on the theory of monotone dynamical systems, see [21] and its recent generalization in [26], and the eigenvalue technique developed and used in [8, 9, 10, 11, 12].

In several cases, models with harvesting can be reduced to non-harvested systems with modified parameters, see [4]. We apply a similar modification which leads to the update of the original space-dependent intrinsic growth rates and carrying capacities. The main focus of the present paper is on the interplay of various diffusion strategies and levels of harvesting for two competing populations. To efficiently handle the model, we introduce the following assumptions:

both competing species have proportional time-independent intrinsic growth rates and the same stationary carrying capacities;

• 2.

one or both species can be subject to space-dependent harvesting, and the spatial harvesting level is proportional to the intrinsic growth rate;

• 3.

diffusion strategies can vary but both populations are isolated, there is no flux through the boundary of the common habitat.

The first assumption that intrinsic growth rates and carrying capacities coincide is quite common for a comparison, whenever one or more other properties are evaluated (in our case, a diffusion strategy and a harvesting level). Harvesting of a percentage of new growth, which is proportional to the current population and the space-dependent intrinsic growth rate, is quite natural in resource management. This means that higher harvesting efforts are applied to more fertile areas with higher harvest expectations. The third assumption that both populations are considered in an isolated area is also general in population dynamics. Evolutionary advantage is evaluated without the interference of emigration and immigration through the boundaries. The system describes an isolated area and can be modeled in laboratory conditions.

Under certain assumptions, an equation describing a harvested species can be reduced to a regular diffusive logistic equation, with an intrinsic growth rate and carrying capacity updated. Combined with non-harvested species, we obtain a system with different diffusion strategies and carrying capacities considered earlier, and distinct intrinsic growth rates, its particular cases were investigated in [7].

Following [7], we can conclude that under any level of harvesting, the advantage of the carrying capacity driven dispersal is no longer sufficient to bring the competing species to extinction and observe, using numerical simulations, how a growing harvesting rate brings exploited population to extinction. However, in the present paper we in addition evaluate harvesting levels for which coexistence and extinction should happen. In particular, explicit sufficient estimates are given for each scenario to occur, unlike [7], where observations on transitions from coexistence to extinction were based on numerical simulations only. Moreover, in the present paper, sufficient coexistence bounds are compared to numerically computed maximal harvesting levels when species coexist.

We also consider the situation when, without harvesting, the two populations sustain. This is attained when the dispersal strategies form an ideal free pair [2, 6] and can describe some specialization (for example, one species foraging primarily in shallow and the other in deep water). In this case, both species can be subject to harvesting. If harvesting efforts are close, there is coexistence. If we fix a harvesting effort for the second species, as harvesting of the first one exceeds this effort and increases, its average density decreases, necessarily bringing it to extinction for intensive enough harvesting. However, for this fixed effort being small enough, lower levels of harvesting and even its absence cannot, generally, bring the second species to extinction. Some estimates when each of the the three scenarios necessarily happens, are presented.

The paper is organized as follows. In Section 2, we state the problem and the main results of the paper. Section 3 illustrates these conclusions with numerical simulations and compares theoretically computed constants to numerical bounds. In Section 4 we discuss conclusions of the present paper and outline some further directions in the study of the influence of dispersal strategies on extinction and survival of exploited populations. Section 5 contains proofs of the main results and all the auxiliary statements required.

2 Preliminaries and Main Results

We consider an isolated spatially distributed competitive system, with harvesting effort proportional to the spatially heterogeneous but time-independent intrinsic growth rate. Both populations exist in a spatial domain, in an isolated environment, there is no flux through the boundary of the domain.

Let $\Omega$ be an open nonempty bounded domain of $\mathbb{R}^{n}$ with $\partial\Omega\in C^{2+p}$, $0<p<1$. We assume that $K(x)$, $K_{i}(x)$, $i=1,2$, $P(x)$ and $Q(x)$ are in the class $C^{1+p}(\overline{\Omega})$, and they are positive for any $x\in\overline{\Omega}$, $r_{i}>0$, $i=1,2$, $r(x)\geq 0$, $x\in\overline{\Omega}$, and $r(x)>0$ in an open nonempty subdomain of $\Omega$. Denote by $u(t,x)$ and $v(t,x)$ population densities, assume that the carrying capacity $K(x)$ and the intrinsic growth rate $r(x)$ are time independent.

The main object of the paper is the system

[TABLE]

For the initial conditions, we assume that $u_{0}(x)\geq 0$, $v_{0}(x)\geq 0$, $x\in\overline{\Omega}$, and both $u_{0}$ and $v_{0}$ are positive in an open nonempty subdomain of $\Omega$.

The type of diffusion in (1), which was discussed in [6], includes the following models as special cases (here we consider the first equation in (1) only):

If $P$ and $a$ are both constant, we obtain the regular diffusion $d\Delta u$, $d>0$.

• 2.

If $a$ is constant, the dispersal term $\Delta(u/P)$ in the particular case when $P=\frac{1}{d}K$ was introduced in [5] and later used in [22, 23], while the constant $P$ and $a=d/K$, $d>0$, leading to $d\nabla\cdot(\frac{1}{K}\nabla v)$, was studied in [23].

• 3.

If $a_{1}=\mu_{1}P$, $r=K$ and $\displaystyle P=e^{\mu_{2}K}$, where $\mu_{i}$, $i=1,2$ are space-independent, we obtain the directed advection $\displaystyle\nabla\cdot\left[\mu\nabla u-\nu u\nabla K\right]$ with $\mu=\mu_{1}$, $\nu=\mu_{1}\mu_{2}$, considered, for example, in [2, 10, 11, 12].

We will also use a modification of (1) without harvesting ($\alpha=\beta=0$) but with different proportional growth rates and carrying capacities

[TABLE]

Further, we explore semi-trivial stationary solutions $(u^{*},0)$ and $(0,v^{*})$, where $u^{*}$ and $v^{*}$ satisfy

[TABLE]

and

[TABLE]

respectively.

Lemma 1**.**

[23*]**

If $P$ is proportional to $K_{1}$ then $u^{*}(x)=K_{1}$ is the only solution of (3). Problem (4) has a unique solution $v^{*}$ which is positive on $\overline{\Omega}$.*

**For (1), a solution can either tend to one of the semi-trivial equilibria, when one of the two species is brought to extinction, or both $u$ and $v$ can persist. In the persistence case, we will say that $u$ and $v$ strongly persist. Numerical simulations illustrate that in this case, for any initial conditions, the solution converges to a certain coexistence equilibrium. Under harvesting of both $u$ and $v$, extinction of both species can occur, but this scenario will be described later. **

The two main results of the present paper are the following.

Theorem 1**.**

Let $P$ be proportional to $K$, while $\displaystyle\nabla\cdot\left[b(x)\nabla\left(\frac{K}{Q}\right)\right]\not\equiv 0$ on $\Omega$, $\alpha,\beta\in[0,1)$. Then, the following three scenarios can occur.

For any $\beta\in[0,1)$ and $\alpha\leq\beta$, all solutions of (1) converge to $((1-\alpha)K,0)$. 2. 2.

For any $\beta<1$, there exists $\alpha_{2}\in(\beta,1)$ such that whenever $\alpha\in(\alpha_{2},1)$, all solutions of (1) converge to $(0,v_{\beta}^{*})$, where $v_{\beta}^{*}$ satisfies (4) with $K_{2}=(1-\beta)K$, $r_{2}=1-\beta$.

For any $\beta\in[0,1)$, there exists $\alpha_{1}\in(\beta,1)$ satisfying

[TABLE]

**where $v_{\beta}^{*}$ is a solution of (4) with $K_{2}=(1-\beta)K$, $r_{2}=1-\beta$, such that for any $\alpha\in(\beta,\alpha_{1})$, all solutions of (1) strongly persist.

---

Remark 1**.**

Theorem 1 outlines the two stages of the influence of harvesting in the case when a harvested species chooses a better dispersal strategy:

Any level of harvesting, whatever small it is, alleviates this advantage and first leads to coexistence. 2. 2.

At a certain level of harvesting, depending on the parameters of equations (carrying capacity, intrinsic growth rate and dispersal $Q$), competitive exclusion of a species being harvested is inevitable.

If both populations are harvested, a more exploited species cannot be the only survivor. Both the lower and the upper bounds of the critical value for a harvesting effort bringing the influenced species to extinction are evaluated in the process of the proof.

The next result refers to the case when $P,Q$ form an ideal free pair, i.e. neither $(K,0)$ nor $(0,K)$ is a solution of (1) with $\alpha=\beta=0$ as

[TABLE]

but some positive multiple of $K$ is in the convex hull of $P$ and $Q$, i.e. there exist $\gamma>0$ and $\delta>0$ such that

[TABLE]

Equality (7) yields that $(\gamma P,\delta Q)$ is a solution of (1).

Theorem 2**.**

Let $P$, $Q$ and $K$ satisfy (6), (7) hold for some $\gamma>0$, $\delta>0$ and $\alpha,\beta\in[0,1)$.

**For any $\beta\in(0,1)$, there exist $\alpha_{1}\in[0,\beta)$ and $\alpha_{2}\in(\beta,1)$ such that for any $\alpha\in(\alpha_{1},\alpha_{2})$, all solutions of (1) strongly persist.

---

For any $\beta\in[0,1)$, there exists $\alpha_{3}\in(\beta,1)$ such that whenever $\alpha\in(\alpha_{3},1)$, all solutions of (1) converge to $(0,v_{\beta}^{*})$, where $v_{\beta}^{*}$ satisfies (4) with $K_{2}=(1-\beta)K$, $r_{2}=1-\beta$. 3. 3.

There exists $\beta_{1}\in(0,1)$ such that for any $\beta\in(\beta_{1},1)$ there is an $\alpha_{0}(\beta)\in[0,\beta)$ such that for any $\alpha\in[0,\alpha_{0})$, all solutions of (1) converge to $(u_{\alpha}^{*},0)$, where $u_{\alpha}^{*}$ satisfies (3) with $K_{1}=(1-\alpha)K$, $r_{1}=1-\alpha$.

Remark 2**.**

If two species form an ideal free pair then harvesting of any of them at a low enough rate does not change the fact of coexistence. However, if we fix a harvesting rate for one species high enough, the absence or low harvesting for its competitor can bring it to extinction. On the other hand, if the competitor is harvested intensively enough, the chosen species becomes the only survivor.

Remark 3**.**

As in Theorem 1, we can present an estimate such that, say, for a fixed $\beta$ and

[TABLE]

all solutions of (1) strongly persist.**

For completeness, we state that over-exploitation brings both species to extinction.

Theorem 3**.**

If $0\leq\beta<1\leq\alpha$ then all solutions of (1) with nontrivial in $v$ initial conditions converge to $(0,v_{\beta}^{*})$, and for $0\leq\alpha<1\leq\beta$, all solutions of (1) with nontrivial in $u$ initial conditions converge to $(u_{\alpha}^{*},0)$. If $\alpha\geq 1$ and $\beta\geq 1$, the zero solution is a global attractor.

Finally, we determine the optimal harvesting policy assuming that we can control dispersal and either one or both of two similar species are appropriate for harvesting. This result generalizes the conclusion for a one-species model in [5]. We recall that, assuming an attractive stationary solution $(u_{s},v_{s})$ of (1), we compute Sustainable Yield as

[TABLE]

If both species are harvested,

[TABLE]

Theorem 4**.**

Maximum Sustainable Yield (MSY) for (1) is attained for $P=K$, $\alpha=0.5$, $\beta>0.5$ and any $Q$ (or $\beta\geq 0.5$ and $\displaystyle\nabla\cdot\left[b(x)\nabla\left(\frac{K}{Q}\right)\right]\not\equiv 0$), if only the first species is harvested (and the second one is subject to culling). If both species are harvested, there is a variety of strategies $P,Q$ such that $K$ is in their convex hull (we choose $P$ and $Q$ such that $P+Q=K$) and $\alpha=\beta=0.5$ leading to MSY. In both cases,

[TABLE]

3 Numerical Simulations

Example 1**.**

Consider (1) with $P(x)=K(x)=2+\cos(\pi x)$, $a(x)=b(x)=Q(x)=1$ for all $x\in(0,4)$, $t>0$. Let $r(x)\equiv 1.1$, $\alpha,\beta\in(0,1)$, and the initial conditions be $u_{0}=2.1=v_{0}$. The system becomes a well-known model of carrying capacity driven diffusion competing with a regularly diffusing population. We consider the solution at the time $T=2000$ for which it is sufficiently close to a steady state. To illustrate the results of Theorem 1, Parts 2 and 3, we consider solutions at $t=T$ as functions of harvesting effort $\alpha\in[0,1)$ for a fixed $\beta\in[0,1)$. We use the functions $K(x)$, $r(x)$, and the stationary solution to the second equation $v_{\beta}^{*}(x)$ to compute $\alpha^{*}$ for each fixed value of the harvesting effort $\beta$ as in (5). In Figure 1 we present a spatial distribution of solutions at $t=T$. Without harvesting, a regularly diffusing population goes extinct, while adding harvesting of $u$ at small rates leads to coexistence.

Figure 2* illustrates the relation between $\alpha^{*}$ computed in (5) for various levels of $\beta$ (which is just a lower estimate of $\alpha$ for which coexistence necessarily occurs) and the maximal value $\alpha^{**}$ of $\alpha_{1}$ when coexistence still occurs computed numerically.*

Figure 3* represents dependence of the average population density on the harvesting rates $\alpha$ and $\beta$.*

Example 2**.**

Let us consider competition model (1), where the carrying capacity is a combination of Gaussian functions and a positive constant $K(x)=10e^{-12.5\pi^{2}(x-2)^{2}}-e^{-50\pi^{2}(x-2)^{2}}+1$ with $P(x)\equiv K(x)$, $a(x)=b(x)=Q(x)=1$ for all $x\in(0,4)$, $t>0$. Let $r(x)\equiv 1.1$, $\alpha,\beta\in(0,1)$, and set the initial conditions as $u_{0}=2.1=v_{0}$. Similarly to the previous example, without harvesting, there is a competitive exclusion of a regularly diffusing population, while harvesting of $u$ at small levels leads to coexistence, as shown in Figure 4.

We provide graphs of average stationary solutions as functions of $\alpha$ for fixed values $\beta$ (Figure 5), using (5) to compute the upper bound $\alpha^{*}$.

Figure 6* presents dependency of the average population density on the harvesting rates $\alpha$ and $\beta$.*

Example 3**.**

Consider the case where $P$ and $Q$ form an ideal free pair, concretely, $P=1.1+0.5\cos(\pi x)$, $Q(x)=0.9+0.5\cos(\pi x)$, $K(x)=2+\cos(\pi x)=P(x)+Q(x)$, $x\in(0,4)$. The initial conditions are $u_{0}=2.1=v_{0}$. We set $r(x)\equiv 1.1$ and consider $\alpha,\beta\in(0,1)$. Note that $P$ and $Q$ are linearly independent. According to [6, Theorem 1] the coexistence solution will be $(P(x),Q(x))$ in the absence of harvesting. To illustrate the results of Theorem 2, Part 2, for each fixed $\beta$ we compute the estimate $\alpha^{*}$ as in (8). We compare population densities at time $t=T=2000$ (see Figure 7) with different (zero and non-zero) values of harvesting, observing the achieved coexistence state.

In Figure 8 the case of solutions at $t=T=2000$ as functions of harvesting $\alpha,\beta$ is illustrated. We construct the dependence of these functions on the value of $\alpha$ for several fixed $\beta\in[0,1)$ indicating the value of $\alpha^{*}$ computed by (8) and $\alpha^{**}$ as the maximal $\alpha_{1}$ when coexistence occurs computed numerically.

In Figure 9, we compare average population densities $u$ and $v$ at time $t=T=2000$ as functions of both $\alpha$ and $\beta$.

Example 4**.**

Figure 10* illustrates how the values of $\alpha$ and $\beta$ that provide coexistence are distributed for an ideal free pair $P(x)=1+10e^{-12.5\pi^{2}(x-2)^{2}}$, $Q(x)=e^{-50\pi^{2}(x-2)^{2}}+2$, $K=10e^{-12.5\pi^{2}(x-2)^{2}}+e^{-50\pi^{2}(x-2)^{2}}+3=P+Q$, $x\in(0,4)$. Similarly, in Figure 11 we show distribution of coexistence promoting rates $\alpha$ and $\beta$ for the case when $P(x)=K(x)=10e^{-12.5\pi^{2}(x-2)^{2}}+e^{-50\pi^{2}(x-2)^{2}}+3$, $Q(x)=e^{-50\pi^{2}(x-2)^{2}}+2$, $x\in(0,4)$. Note that in Figure 10 there is a region of values near $(0,0)$ for which coexistence is preserved, while in the case depicted on Figure 11, only increasing $\alpha$ allows to avoid competitive exclusion ($\beta=0$).*

4 Discussion

The results obtained can shed some light on the species invasion management. If a more efficient diffuser is an invader, all other parameters being the same, without harvesting or for the same or higher level of harvesting of the resident population, its habitat will be completely populated by invaders bringing the original species to extinction. However, any prevalence in harvesting of the invader can guarantee survival for the resident population; moreover, some estimates were presented when the invader goes extinct. If a resident is a more efficient diffuser, its habitat can be invaded only with more intensive harvesting of the resident population. An estimate on harvesting rate is given when both species should coexist, in spite of the fact that the resident is exploited more than the invader. An additional estimate provides the level of harvesting guaranteeing that the habitat is completely invaded bringing the resident to extinction. If invader-resident form an ideal free pair, for close levels of harvesting coexistence is guaranteed, while harvesting levels were evaluated bringing the invader (or the resident) to extinction.

The estimation for the harvesting effort leading to coexistence is important, as it may give an indication which harvesting levels lead to sustainability of the exploited population. However, the proposed estimate, though approaching the switch value as the harvesting effort for the competitor increases, is still far from being sharp. Also, if culling is applied to bring one of the species to extinction, the minimal required reduction rate is overestimated.

The paper also explores the situation when, in the absence of harvesting, the two species coexist, based on partial specialization (expressed in the dispersal strategies) and resources sharing. In this case, once the exploitation rate of one of the populations is fixed, we obtain a limitation on the harvesting level of the other population, such that the two species still coexist. Usually, this assumes a close level of harvesting for both species. If there is an ideal free pair, harvesting of one of the species has to be high enough for a possibility of the harvested species to be brought to extinction. For substantial efforts, even some smaller harvesting rates for the other species can be allowed, still leaving this population as the only survivor in the competition. Also, as a harvested system is reduced to a model with modified growth rates and carrying capacities, we explored the influence of these parameters, changed proportionally, on the competition outcome.

The general approach to harvesting problems allows to extend research of the present paper in the following direction. Here we considered external influence in the form of harvesting only. More powerful control would involve stocking along with harvesting. For example, if a more efficient diffuser tries to invade a habitat of a less efficient resident, stocking of the resident species can prevent it from extinction. If a resident is in additional disadvantage being a preferred prey for a local predator, an additional stocking effort is required to avoid its extinction. If stocking is allowed in addition to harvesting, for the ideal free pair and any (even zero) harvesting effort, the population can be brought to extinction by a competitor having an excess of external resources. These results, we believe, can be justified in the same way as Theorems 1 and 2. **While numerical simulations illustrate global stability of a unique coexistence equilibrium, whenever both semi-trivial equilibria is unstable, we do not state its uniqueness and global attractivity. Establishing this fact is still an open problem. **

As a combination of the present paper and [7], it is possible to incorporate harvesting in the system where one of the species more efficiently uses resources than the other one expressed in a higher carrying capacity. If the harvesting effort is constant, after appropriate substitutions, Theorems 1, 2, 3 and 4 would lead to the description of possible outcomes. We expect that the two species would coexist for a harvesting effort of the efficient consumer in a certain interval. For smaller efforts, a less efficient consumer goes extinct, while for higher efforts it becomes the only survivor in the competition.

It would be interesting to extend the present analysis to the case when harvesting is semi-constant [31], being constant everywhere but for small population levels, where it is close to proportional to the population size. Such an assumption allows us to avoid the situation when solutions become negative. Certainly, to apply the same results, we have to suppose that, though independent of the population size, the harvesting function is spatially heterogeneous, being proportional to both the growth rate and the carrying capacity. The study of significant and remnant harvested populations (introduced in [31]) will strongly depend on the bounds of the carrying capacity.

Note that the results of the present paper are based on the essential assumption that the carrying capacity is positive on the whole domain. However, both in literature and practical applications, the carrying capacity can become zero at some locations, especially close to the boundary of the domain. Zero or even negative growth rates can be considered, see, for example, [24]. For the type of diffusion discussed in the present paper, the carrying capacity driven diffusion, existence of areas separated by the strips with the zero carrying capacities can lead to separated non-mixing patches, while the other species would freely spread in the whole domain. It would be interesting to explore how a complicated structure of the domain, in particular, quickly changing $r$ and $K$, can influence survival. In Example 2, high variation of carrying capacities promotes survival of the species with a carrying capacity driven diffusion.

Unlike [23], in the present paper both species are subject to the logistic type of growth. However, the case when there is a strong Allee effect, i.e. there is a minimum population size guaranteeing survival, even without harvesting, is still an open problem. The influence of the growth law on harvesting outcome is outlined in [22]: higher per capita growth rates for smaller population densities are crucial. It would also be interesting to explore competition of the two species with different growth laws.

5 Auxiliary Results and Proofs

We start with a result on a stationary solution weighted average.

Lemma 2**.**

[6*, Lemma 3]**

If $v^{\ast}$ is a positive solution of (4) then*

[TABLE]

As a partial case of (2), consider a competition model with proportional, space-dependent carrying capacities and growth rates with $r_{1}>0$, $r_{2}>0$

[TABLE]

Since (11) is a monotone dynamical system [8], as well as all the above systems, we will further use a modification of [21, Theorem B, p. 4087] in the form earlier described in [23, Theorem 17, pp.1203-1204] in detail.

Lemma 3**.**

*Let $(0,0)$ be a repelling equilibrium of (11). Then exactly one of the following holds: *

(a)

There exists a positive coexistence equilibrium of (11).

(b)

$(u,v)\rightarrow(u^{\ast},0)$* as $t\rightarrow\infty$ for $(u_{0},v_{0})$ from a prescribed set.*

(c)

$(u,v)\rightarrow(0,v^{\ast})$* as $t\rightarrow\infty$ for $(u_{0},v_{0})$ from a prescribed set.*

Moreover, if (b) or (c) holds then for every non-trivial initial conditions either $(u,v)\rightarrow(u^{\ast},0)$ or $(u,v)\rightarrow(0,v^{\ast})$ as $t\rightarrow\infty$.

Remark 4**.**

Note that by Lemma 3, once the trivial equilibrium is a repeller, there is no coexistence equilibrium and one of semi-trivial equailibrium solutions is unstable, the other semi-trivial equilibrium is globally asymptotically stable.

If zero is a repeller and both semi-trivial equailibrium solutions are unstable, all solutions converge to a coexistence equilibrium.

We start with auxiliary statements used in the proof of Theorem 1, in particular, we use special cases of the results in [7].

Lemma 4**.**

[7*]**

If $P$ is proportional to $K$, $\nabla\cdot[b\nabla(K/Q)]\not\equiv 0$ on $\Omega$ and $c\in[1,\infty)$ then the semi-trivial equilibrium $(cK,0)$ of system (11) is globally asymptotically stable.*

Lemma 5**.**

[7*]**

If $c\in(0,1)$ then the semi-trivial equilibrium $(cK,0)$ of (11) is unstable.*

Lemma 6**.**

[7*]**

Let $K_{1}$ and $K_{2}$ be positive on $\overline{\Omega}$, $r$ be nonnegative in $\overline{\Omega}$ and positive in an open nonempty subdomain of $\Omega$, $r_{1}>0$, $r_{2}>0$. Then the trivial equilibrium $(0,0)$ of (2) is a repeller.*

**However, in the case $c\in(0,1)$, the paper [7] does not give an answer whether the two species would eventually coexist, or only the second one survives while the first one goes extinct. In [7], numerical simulation illustrated that for $c<1$ close enough to one, there is coexistence, while for $c$ sufficiently small, the semi-trivial equilibrium $(0,v^{*})$ is globally asymptotically stable. Here we prove existence of such $c^{*}$ that coexistence is inevitable for any $c\in(c^{*},1)$. **

Lemma 7**.**

Let $P$ be proportional to $K$, while $\nabla\cdot\left[b(x)\nabla\left(\frac{K(x)}{Q(x)}\right)\right]\not\equiv 0$ on $\Omega$. Then there exists $c_{1}\in(0,1)$ such that for $c\in(c_{1},1)$, the semi-trivial equilibrium $(0,v^{*})$ of (11), where $v^{*}$ satisfies (4) with $K_{2}=K$, is unstable.

Proof.

For the linearization of the first equation in (11) around the equilibrium $(0,v^{*})$

[TABLE]

consider the associated eigenvalue problem

[TABLE]

According to [8], the principal eigenvalue $\sigma_{1}$ of (12) can be computed as

[TABLE]

Let $c=1$, $\sigma_{1}$ be the principal eigenvalue, and $\Psi$ be the principal eigenfunction, then

[TABLE]

We have

[TABLE]

Note that $v^{*}$ does not depend on $c$, and for

[TABLE]

the principal eigenvalue of (12) is positive, since the substitution of $\psi=\Psi$ in the brackets of (13) gives a positive number.

Further, we construct a simpler estimate of $c$ guaranteeing coexistence, not involving $\Psi$.

Choosing $\psi(x)=K(x)$, denoting $\displaystyle M:=\int\limits_{\Omega}K(x)\,dx,$ and recalling that $P$ is proportional to $K$, we observe by Lemma 2 that the principal eigenvalue satisfies

[TABLE]

for $c\in(c^{*},1)$, where

[TABLE]

Therefore, for $c\in(c^{*},1)$, $\sigma_{1}$ is positive, and thus the semi-trivial steady state $(0,v^{*})$ of (11) is unstable. ∎

All the above results on local instability of semi-trivial equilibrium solutions and the fact that the trivial equilibrium is a repeller, by Lemma 3, imply the following result.

Proposition 1**.**

*Let $P$ be proportional to $K$, while $\displaystyle\nabla\cdot\left[b(x)\nabla\left(\frac{K}{Q}\right)\right]\not\equiv 0$ on $\Omega$. ***There exists $c_{1}\in(0,1)$ such that $c_{1}\leq c^{*}$, where $c^{*}$ is defined in (14), and for $c\in(c_{1},1)$ all solutions of (11) strongly persist. ****

Lemma 8**.**

Let $P$ be proportional to $K_{1}=K$, let $\displaystyle\nabla\cdot\left[b(x)\nabla\left(\frac{K}{Q}\right)\right]\not\equiv 0$ on $\Omega$. There exists $c_{2}\in(0,1)$ such that whenever $c\in(0,c_{2})$, system (11) has no coexistence equilibrium.

Proof.

Everywhere we assumed $K(x)\in C^{1+p}(\overline{\Omega})$, $K(x)>0$ for all $x\in\overline{\Omega}$. These assumptions on $K(x)$ imply existence of positive lower and upper bounds for $K$

[TABLE]

By [22, Appendix Proposition 4] and [8, Propositions 3.2 and 3.3], there exist similar positive bounds for a positive stationary solution $v^{*}(x)$ of (4)

[TABLE]

We have $b\leq M$; assuming the contrary, we obtain that $v^{*}(x)>K(x)$ for any $x$, thus integration of (4) with $K_{2}\equiv K$ over $\Omega$ leads to a contradiction that the integral of a negative on $\Omega$ function equals zero. Thus $b/(2M)<1$; let us choose $c\in\left(0,\frac{b}{2M}\right)$. Let $(u_{s},v_{s})$ be a coexistence solution of (11). Substituting $\overline{u}=\frac{bK(x)}{2M}$ into the first equation in (11), assuming $v_{s}\geq 0$ and taking into account that the diffusion term vanishes as $\overline{u}/P$, or $\overline{u}/K$, is constant, we obtain

[TABLE]

for any $x\in\Omega$ since $c<\frac{b}{2M}$. Thus $\overline{u}=\frac{bK}{2M}$ is an upper solution of the equation, see [30, Chapter 8, Definition 1.2]. We can choose zero as a lower solution and immediately get by (15), (16),

[TABLE]

Choosing $\underline{v}=\frac{v^{*}}{2}$ in the second equation of (11), we can see that it is a lower solution (see [30, Chapter 8, Definition 1.2]), since $v^{*}$ is a solution of (4) with $K_{2}\equiv K$:

[TABLE]

Therefore, $v_{s}\geq\frac{v^{*}}{2}$. Finally, using this inequality together with the fact that $u_{s}\geq 0$ in the first equation of (11), we get, due to $c<\frac{B}{2M}$,

[TABLE]

where the strict inequality is valid in the first and the last cases, unless $u_{s}(x)\equiv 0$. Integrating both sides of the inequality above, using the boundary conditions and assuming $u_{s}\not\equiv 0$ leads to

[TABLE]

which is a contradiction. Hence, there is no coexistence equilibrium. ∎

Now, Lemma 3 leads to the global stability of $(0,v^{*})$, as stated below.

Proposition 2**.**

Let $P$ be proportional to $K$, while $\displaystyle\nabla\cdot\left[b(x)\nabla\left(\frac{K}{Q}\right)\right]\not\equiv 0$ on $\Omega$. There exists $c_{2}\in(0,1)$ such that for $c\in(0,c_{2})$ all solutions of (11) converge to the semi-trivial equilibrium $(0,v^{*})$.

Now, we can proceed to the proof of the main result.

Proof.

(the proof of Theorem 1) Let $\beta\in[0,1)$. Then (1) becomes (11) with

[TABLE]

1. If $\alpha\leq\beta$, we immediately get $c\geq 1$. By Lemma 4, the semi-trivial equilibrium $(u^{\ast},0)=((1-\alpha)K,0)$ is globally asymptotically stable.

2. Fix $\beta\in[0,1)$. By Proposition 1, there is $\displaystyle c_{1}=\frac{1-\alpha_{1}}{1-\beta}\in(0,1)$ such that for any $c\in(c_{1},1)$, all solutions of (11) strongly persist, where $c_{1}\leq c^{\ast}$, $c^{\ast}$ is defined in (14). Since $\alpha_{1}=1-c_{1}(1-\beta)$, taking into account that $K$ is substituted by $(1-\beta)K$ and $v^{\ast}=v_{\beta}^{\ast}$ is also $\beta$-dependent, the bound $\alpha^{\ast}=1-c^{\ast}(1-\beta)$ has form (5). For $\alpha\in(\beta,\alpha_{1})$, which corresponds to $c\in(c_{1},1)$, by Proposition 1, all solutions strongly persist.

3. Again, for a fixed $\beta\in[0,1)$, system (1) becomes (11) with (18). By Proposition 2, there exists a $\displaystyle c_{2}=\frac{1-\alpha_{2}}{1-\beta}\in(0,1)$ such that for any $c\in(0,c_{2})$, all solutions of (11) and thus all solutions of (1) converge to $(0,v_{\beta}^{\ast})$. Choosing $\displaystyle\alpha_{2}=1-c_{2}(1-\beta)\in(0,1)$, we get that for $\alpha\in(\alpha_{2},1)$, all solutions of (1) converge to $(0,v_{\beta}^{\ast})$, which concludes the proof. ∎

Remark 5**.**

If we consider $\beta\geq 1$ then, from the second equation in (1), $\displaystyle\lim_{t\to\infty}v(t)=0$, independently of $u(t)\geq 0$. The first equation implies, for any $\alpha\in[0,1)$, $\displaystyle\lim_{t\to\infty}u(t)=(1-\alpha)K$, thus the semi-trivial equilibrium $((1-\alpha)K,0)$ is globally asymptotically stable.

Next, we proceed to the case of the ideal free pair when (6) is satisfied and for some $\gamma>0$ and $\delta>0$, identity (7) holds.

Let us present auxiliary statements used in the proof of Theorem 2. We further use the following lemma obtained in [28] and recently discussed in [16].

Lemma 9**.**

[28*]**

If $u^{*}/K_{1}\not\equiv 1$ then the solution $u^{*}$ of (3) satisfies*

[TABLE]

Similarly, if $v^{*}/K_{2}\not\equiv 1$ then for the solution $v^{*}$ of (4) the following inequality holds

[TABLE]

Remark 6**.**

If in (3) we assume constant $a$ and $P$ and $r=K$, we get the stationary solution of the Neumann problem for Fisher’s equation

[TABLE]

for which (19) leads to an important conclusion

[TABLE]

that the average solution levels are higher than that of the carrying capacity, while (10) leads to the opposite result for space-independent $r$.

Lemma 10**.**

Let (6) and (7) hold for some $\gamma>0$, $\delta>0$, and $c\in(0,1]$. Then the semi-trivial equilibrium $(u^{*},0)$ of (11) is unstable.

Proof.

Integrating (3) with $K_{1}=cK$ and using the boundary conditions for $u^{*}$, we get

[TABLE]

Here a substitution of $u^{*}$ in (3) with $K_{1}=cK$, where (6) holds, leads to the conclusion that $u^{*}-cK\not\equiv 0$. Thus the first term in the right-hand side above is negative. By (19) in Lemma 9, the second term is also negative. Since the sum equals to zero,

[TABLE]

Moreover, as $c\in(0,1]$ and

[TABLE]

with $(1-c)/c\geq 0$, the first term in the right-hand side of the equality above is positive

[TABLE]

Next, consider the eigenvalue problem for the linearization of the second equation in (11) around the equilibrium $(u^{*},0)$

[TABLE]

According to [8], the principal eigenvalue $\sigma_{1}$ of (22) can be computed as

[TABLE]

It is not less than the value computed for $\psi=Q$ which is

[TABLE]

by (21). Since the principal eigenvalue is positive, the equilibrium $(u^{*},0)$ is unstable. ∎

Lemma 11**.**

Let (6) and (7) hold for some $\gamma>0$, $\delta>0$. Then there exists $c_{1}\in(0,1)$ such that for $c\in(c_{1},1)$, the semi-trivial equilibrium $(0,v^{*})$ of (11) is unstable.

Proof.

Integrating (4) over $\Omega$, taking into account the boundary conditions, gives

[TABLE]

where the first and the third terms are negative.

Thus the second term is positive, and $\displaystyle\int_{\Omega}rP\left(1-\frac{v^{*}}{K}\right)~{}dx>0.$ Moreover, denoting

[TABLE]

we obtain that

[TABLE]

Consider eigenvalue problem (12) for the linearization of the first equation in (11). Its principal eigenvalue $\sigma_{1}$ is given by (13), and it is not less than the value computed for $\psi=P$, i.e.

[TABLE]

for $c\in(c_{1},1]$, and $c_{1}\leq c^{*}$, by (24), where $c^{*}$ is defined in (23). Since the principal eigenvalue is positive, the equilibrium $(0,v^{*})$ for $c\in(c^{*},1)$ is unstable. ∎

Lemma 12**.**

Let (6) and (7) hold for some $\gamma>0$, $\delta>0$. There exists $c_{2}\in(0,1)$ such that whenever $c\in(0,c_{2})$, system (11) has no coexistence equilibrium.

Proof.

As we assumed $K,P\in C^{1+p}(\overline{\Omega})$, $K(x)>0$, $P(x)>0$ for all $x\in\overline{\Omega}$, there is an upper bound for $K$ in (15), $M:=\max_{x\in\overline{\Omega}}K(x)$. Similar smoothness estimates are valid for $P(x)$ and the ratio $K/P$, and the bounds are positive

[TABLE]

Also, as mentioned in the proof of Lemma 8, the lower bound for a positive stationary solution $v^{*}(x)$ of (4) with $K_{2}=K$ is positive $0<b:=\min_{x\in\overline{\Omega}}v^{*}(x)$, see (16). Again, $b\leq M$, $b/(2M)<1$ and, by the definition of $g$ and $G$, $\displaystyle g\leq\frac{K(x)}{P(x)}\leq G$. We choose

[TABLE]

Then $\displaystyle c<\frac{b}{2M}$. Let $\displaystyle\overline{u}(x)=\frac{bgP(x)}{2M}$. Substituting $\overline{u}$ into the first equation in (11), assuming $v_{s}\geq 0$ and noticing that the diffusion term vanishes as $\overline{u}/P$ is constant, we obtain

[TABLE]

for any $x\in\Omega$. Thus $\overline{u}=\frac{bgP(x)}{2M}$ is an upper solution, choosing zero as a lower solution, we get

[TABLE]

as in (17). Choosing $\underline{v}=\frac{v^{*}}{2}$ in the second equation of (11), we can see that it is a lower solution, see [30, Chapter 8, Definition 1.2], since $v^{*}$ is a solution of (4) with $K_{2}\equiv K$:

[TABLE]

Thus, $v_{s}\geq\frac{v^{*}}{2}$. Finally, using this inequality together with the fact that $u_{s}\geq 0$ in the first equation of (11), we get

[TABLE]

and the strict inequality is valid unless $u_{s}\equiv 0$. Integrating both sides of the inequality above using boundary conditions and assuming $u_{s}\not\equiv 0$ leads to

[TABLE]

The contradiction proves that there is no coexistence equilibrium. ∎

Similarly to Lemma 11, we establish stability for $c\in[1,\infty)$.

Lemma 13**.**

Let (6) and (7) hold for some $\gamma>0$, $\delta>0$. Then there exists $c_{1}\in(1,\infty)$ such that for $c\in[1,c_{1})$, the semi-trivial equilibrium $(u^{*},0)$ of (11) is unstable.

By Lemma 11, for certain $c$, $(0,v^{*})$ is unstable, while by Lemma 13, for another set of $c$, $(u^{*},0)$ is unstable. Whenever $c$ satisfies both conditions, Lemma 3 implies coexistence conditions.

Proposition 3**.**

Let (6) and (7) hold for some $\gamma>0$, $\delta>0$. There exist $c_{1}\in(0,1)$ and $c_{2}\in(1,\infty)$ such that whenever $c\in(c_{1},c_{2})$, all solutions of (11) strongly persist.

Next, competitive exclusion can be considered for $c>1$ as well.

Lemma 14**.**

Let (6) and (7) hold for some $\gamma>0$, $\delta>0$. There exists $c_{3}\in(1,\infty)$ such that whenever $c\in(c_{3},\infty)$, system (11) has no coexistence equilibrium.

Based on Lemmata 12 and 14, we establish stability of the semi-trivial equilibria.

Proposition 4**.**

Let (6) and (7) hold for some $\gamma>0$, $\delta>0$. There exist $c_{1}\in(0,1)$ and $c_{2}\in(1,\infty)$ such that for $c\in(0,c_{1})$ all solutions of (11) converge to the semi-trivial equilibrium $(0,v^{*})$, while for $c\in(c_{2},\infty)$ all solutions of (11) converge to the semi-trivial equilibrium $(u^{*},0)$.

At this stage, we can proceed to the proof of Theorem 2.

Proof.

(the proof of Theorem 2) We recall that (1) has form (11) with notations (18).

1. Fix $\beta\in(0,1)$. By Proposition 3, there are $c_{1}\in(0,1)$ and $c_{2}\in(1,\infty)$ such that for $c\in(c_{1},c_{2})$, all solutions strongly persist. Denote $\alpha_{1}=1-\max\{1-c_{2}(1-\beta),0\}$, $\alpha_{2}=1-c_{1}(1-\beta)$. Then $0\leq\alpha_{1}<\beta<\alpha_{2}<1$ and, by Proposition 3, for $\alpha\in(\alpha_{1},\alpha_{2})$, all solutions of (1) strongly persist.

2. For any fixed $\beta\in[0,1)$, by Proposition 4, there is a $c_{1}\in(0,1)$ such that for $c\in(0,c_{1})$, all solutions of (11) converge to $(0,v_{\beta}^{*})$. Denote $\alpha_{3}=1-c_{1}(1-\beta)\in(\beta,1)$. Then for $\alpha\in(\alpha_{3},1)$ we get $c\in(0,c_{1})$ and thus all solutions of (1) converge to $(0,v_{\beta}^{*})$, where $v_{\beta}^{*}$ is a solution of (4) with $K_{2}=(1-\beta)K$, $r_{2}=1-\beta$.

3. Denote $\displaystyle\beta_{1}=1-\frac{1}{c_{2}}$, where $c_{2}>1$ is defined in Proposition 4, such that for all $c\in(c_{2},\infty)$, all solutions of (11) converge to $(u_{\alpha}^{*},0)$. Next, let $\beta\in(\beta_{1},1)$ be fixed, and $\alpha_{0}=1-(1-\beta)c_{2}$. Then $c_{2}>1$ and $c\in(c_{2},\infty)$ correspond to $0<\alpha_{0}<\beta<1$ and $\alpha\in[0,\beta)$, respectively. Thus for any $\alpha\in[0,\alpha_{0})$, all solutions of (1) converge to $(u_{\alpha}^{*},0)$, where $u_{\alpha}^{*}$ satisfies (3) with $K_{1}=(1-\alpha)K$, $r_{1}=1-\alpha$. This concludes the proof. ∎

Theorem 3 can be justified, similarly to [22], using differential inequalities.

Proof.

(the proof of Theorem 4) We recall that for a fixed $K$, the maximum of $f(x)=x(1-x/K)$ is attained at $x=K/2$ and equals $\frac{1}{4}K$. Integrating the first equation in (1) over $\Omega$ and taking into account the boundary conditions, we get

[TABLE]

For either $P=K$, $\alpha=0.5$, $\beta>0.5$ and any $Q$, or $\beta\geq 0.5$ and $\displaystyle\nabla\cdot\left[b(x)\nabla\left(\frac{K}{Q}\right)\right]\not\equiv 0$, the stationary solution $(0.5K(x),0)$ leads to (9).

Further, in the case of the ideal free pair, MSY in (9) is attained for appropriate $P$ and $Q$, $P+Q=K$, $\alpha=\beta=0.5$. Similarly to the above case

[TABLE]

which concludes the proof. ∎

Remark 7**.**

Let us note that the maximum $\displaystyle\int_{\Omega}u(x)~{}dx$ may be attained for $\displaystyle\nabla\cdot\left[a(x)\nabla\left(\frac{K}{P}\right)\right]\not\equiv 0$, see (20) in Remark 6.

Acknowledgment

**The authors are grateful to the referees whose valuable comments significantly contributed to the paper presentation. ** The authors were partially supported by the NSERC research grant RGPIN-2015-05976. **The second author was also supported by Pacific Institute for the Mathematical Sciences with PIMS Student Training Acceleration Award. **