A time-periodic reaction-diffusion population model of a species subjected to harvesting on the boundary of a protected habitat

Abstract

The objective in this paper is to determine analytically a maximally sustainable pro-rata density-dependent harvest rate of a hypothetical biological species on the spatial boundary of a protected habitat (i.e. no harvesting of the species is allowed in the interior of the spatial domain). This is achieved by analysing an initial-boundary value model involving a reaction-diffusion equation in which the reaction term is a concave function of the population density, depends periodically on time, and varies spatially. The model is equipped with Robin boundary conditions modelling a continuum of potential pro-rata density-dependent harvest rates on the spatial boundary. A set of necessary and sufficient criteria for the existence of a unique positive periodic attractor of solutions to the model is established by employing the theory of comparison principles and monotone iteration schemes. Thereafter, a long-time asymptotic analysis of the population density is undertaken by invoking classical results from the theory of eigenproblems. In this way, the aforementioned necessary and sufficient conditions are rendered more practical in terms of an upper bound on the maximal allowable pro-rata density-dependent harvest rate. Moreover, important properties of the time-periodic attractor of solutions are established to guarantee the existence of a density pro-rata harvest rate which maximises the total harvest per unit time at equilibrium.

1 Introduction

It is well known that as the global population increases, so too does its demand for natural resources—both on land and in the oceans. As a result, the stock of marine fish have, for instance, been in sharp decline since the early 1900s as a direct result of over-harvesting these populations [11, 17]. In recent years, marine ecologists and fishery industrialists have debated the necessity of enforcing harvesting quotas or no-harvesting protection zones to protect marine species in a bid to circumvent the collapse of fishing industries [2, 9]. In an attempt to contribute objectively to the aforementioned debate, mathematical models are often derived which are capable of describing the spatio-temporal evolution of biological species.

There are two commonly adopted approaches in such a modelling endeavour, namely the derivation of abstract models and the construction of practical models [13]. Abstract models are derived with a conscious effort directed at keeping assumptions to a minimum, resulting in an apparently simple mathematical model, allowing for an elegant mathematical treatment of the model [18]. In contrast, practical models are based on realistic assumptions and usually involve inelegant parametrisations of the interrelationships of a large number of variables. Simulations may be employed to capture the long-term spatio-temporal behaviours of solutions to such models in special cases [18]. In this paper, we focus squarely on an abstract model, more specifically—a reaction-diffusion model.

There have been a number of studies in which reaction-diffusion models were employed to determine the effects of harvesting on the prediction of persistence or extinction of a biological species [5,6,7,8, 12, 14]. The objectives of these studies have included finding maximally allowable harvest rates distributed over the interior of the species habitat and determining the size of a protection zone within the habitat such that unrestricted harvesting may occur in the interior of the habitat outside of this protection zone. There have also been studies aimed at finding a spatial harvesting distribution in the interior of the habitat which maximises the total catch of the species over time at equilibrium. In contrast, there have been comparatively few studies in which harvesting is allowed only on the boundary of an otherwise protected zone. We consider this type of harvesting scenario in the current paper. We show that the creation of a protection zone per se may not be enough to guarantee survival of the species if it is harvested too heavily on the spatial boundaries of this protection zone.

We consider a single biological population exploiting natural resources for the purpose of survival in a bounded, m-dimension isotropic spatial domain $\Omega$ with a smooth boundary $\partial \Omega$ (typically, in applications, $m=2$ or 3). Let $u({\varvec{x}},t)$ denote the population density of the population at position ${\varvec{x}} = [x_1, \ldots , x_m]^T \in \Omega$ at time $t \in [0,\infty )$. Then the spatio-temporal evolution of the population density u may be modelled abstractly by the initial-boundary value problem

$$\begin{aligned}&\text{ D.E. } \quad {\displaystyle \frac{\partial u}{\partial t} = \Delta u + g(u,\varvec{x},t)}, \quad {\varvec{x}}\in \Omega \setminus \partial \Omega , \quad t \ge 0, \end{aligned}$$

(1)

$$\begin{aligned}&\text{ B.C. }\quad {\displaystyle \frac{\partial u}{\partial \varvec{\eta }}+ \kappa u = 0}, \qquad \qquad \quad ~{\varvec{x}} \in \partial \Omega , \quad ~ \qquad t \ge 0, \end{aligned}$$

(2)

$$\begin{aligned}&\text{ I.C. } \quad u(\varvec{x},0) = f(\varvec{x})>0, \qquad ~{\varvec{x}}\in \Omega , \end{aligned}$$

(3)

where $\Delta$ is the Laplace operator, $\kappa$ is a positive constant called the pro-rata boundary harvest rate, $\partial /\partial \varvec{\eta }$ is the directional derivative of u in the direction of ${\varvec{\eta }}$, an outward-pointing unit normal vector on the boundary surface $\partial \Omega$ of $\Omega$, and $f(\varvec{x})$ is a prescribed initial condition with a finite p-norm. In this paper we assume the reaction term $g(u,\varvec{x},t)$ to be Lipschitz continuous in u, measurable in $\varvec{x}$ and periodic in t. In particular, it is assumed that there exists a positive real number T such that the function g is T-periodic in t, in the sense that

$$\begin{aligned} g(u,\varvec{x},t)=g(u,\varvec{x},t+kT)\text { for all }u\ge 0, t\ge 0\text { and }\varvec{x}\in \Omega ,\text { for any }k\in {\mathbb {N}}. \end{aligned}$$

(4)

It is also assumed that

$$\begin{aligned} g(u,\varvec{x},t)\text { and }\frac{\partial g}{\partial u}(u,\varvec{x},t) \end{aligned}$$

are Lipschitz continuous and of class $C^\alpha (\Omega \times [0,T])$ uniformly for u in a bounded subset of ${\mathbb {R}}$, and that

$$\begin{aligned} g(0,\varvec{x},t) = g(1,\varvec{x},t)=0\text { and }\frac{\partial ^2g}{\partial u^2}(u,\varvec{x},t)<0 \end{aligned}$$

(5)

for all $u>0$, $\varvec{x}\in \Omega$ and any $t>0$, so that $u\equiv 0$ is a subsolution of (1)–(3) and that any constant greater than 1 is a supersolution of (1)–(3), implying that the state space $u\in [0,1]$ is positively invariant. The notions of sub- and supersolutions are defined as in Theorem 1.19 of [3].

In [19], we considered a special case of (1)–(3) where $\Omega$ was taken as the unit interval and the reaction term g was allowed to depend only on the local population density (i.e. it was explicitly independent of space and time). In this paper, we consider the more general situation where the habitat $\Omega$ may consist of any finite number of spatial dimensions, where the reaction term g is spatially heterogeneous, and depends periodically in time. In this context, we establish a maximally sustainable pro-rata harvest rate for which a globally stable periodic attractor of the model exists, attracting all other solutions. Moreover, we establish important properties of this periodic attractor guaranteeing the existence of a pro-rata harvest rate which maximises the total harvest over a single interval of periodicity. The model we considered in [19] and the more general model in this paper are different from those in previous studies where reaction-diffusion models were employed to determine the effects of species harvesting on the prediction of persistence or extinction of a biological species. The main difference is that in this paper we assume that the habitat $\Omega$ is protected and that the population may only be harvested on the spatial boundary $\partial \Omega$. We are not aware of literature pertaining models where the population is only harvested on boundary of its habitat $\partial \Omega$.

The remainder of this paper is organised as follows. A necessary and sufficient condition is established in Sect. 2 (in terms of the sign of the principal eigenvalue) for the existence of a unique positive T-periodic attractor of solutions to (1)–(3) corresponding to any non-zero initial condition $f(\varvec{x})$. This is followed by a long-term asymptotic analysis of the solution to (1)–(3), in which the necessary and sufficient condition of Sect. 2 for the existence of a unique positive, T-periodic attractor of solutions to (1)–(3) is rendered more practical in Sect. 3 by relating it to a maximal allowable pro-rata boundary harvest rate $\kappa ^*$. Thereafter, we establish a number of fundamental properties of the aforementioned periodic attractor in Sect. 4. These properties include establishing that the periodic attractor of solutions to (1)–(3) is spatially symmetric, provided that $\Omega$ and $g(u,\varvec{x},t)$ are also spatially symmetric, and that this periodic attractor is a monotonically decreasing function of the pro-rata harvest rate $\kappa$. Finally, we establish the existence of a pro-rata harvest rate in Sect. 5 which maximises the total harvest yield over a single interval of periodicity. The paper finally closes in Sect. 6 with a brief summary of its contents.

2 Existence of a unique positive periodic solution

The function $g(u,\varvec{x},t)$ may be expressed as a Taylor expansion around the point $u\equiv 0$ as

$$\begin{aligned} g(u,\varvec{x},t) = g(0,\varvec{x},t) + \frac{\partial g}{\partial u}(0,\varvec{x},t) u + g_1(u,\varvec{x},t)u^2 \end{aligned}$$

(6)

for some negative, continuous function $g_1(u,\varvec{x},t)$ as a result of the conditions in (5), upon which the questions of the existence and stability of a positive T-Periodic solution to (1)–(3) may be settled in terms of the sign of the principal eigenvalue $\mu =\mu _1$ of the corresponding linearised periodic eigenvalue problem

$$\begin{aligned}&\frac{\partial \phi }{\partial t} - \Delta \phi -b(\varvec{x},t)\phi =\mu \phi ,\quad \varvec{x}\in \Omega \setminus \partial \Omega , \end{aligned}$$

(7)

$$\begin{aligned}&\frac{\partial \phi }{\partial \varvec{\eta }}+\kappa \phi = 0, \quad \qquad \qquad \qquad \quad \varvec{x}\in \partial \Omega , \end{aligned}$$

(8)

$$\begin{aligned}&\phi T\text {-periodic in }t, \end{aligned}$$

(9)

where

$$\begin{aligned} b(\varvec{x},t) = \frac{\partial g}{\partial u}(0,\varvec{x},t). \end{aligned}$$

(10)

It is known, for any function $g(u,\varvec{x},t)$ for which $b(\varvec{x},t)$ in (10) is T-periodic, that the principal eigenvalue $\mu =\mu _1$ of (7) is real and that its corresponding eigenfunction $\phi _1(\varvec{ x},t)$ is positive on $\Omega \times (0,\infty )$ [10, 16]. The stability properties of the trivial solution $u\equiv 0$ in terms of the principal eigenvalue of (7)–(9), is summarised as follows.

Lemma 2.1

(Stability of the trivial solution to (1)–(3)  [10, p. 67]) If the principal eigenvalue $\mu =\mu _1$ of (7)–(9) is strictly positive, then the trivial solution $u\equiv 0$ to (1)–(3) is linearly stable and hence locally asymptotically stable. If, however, this eigenvalue is strictly negative, then the aforementioned trivial solution is linearly unstable and hence globally unstable.

The following properties of sub- and supersolutions of (1)–(3) are required in order to derive a general result on the existence and asymptotic stability of a time-periodic solution to (1)–(3).

Lemma 2.2

(Sub- and supersolutions  [3, p. 80])   If  ${\overline{u}}$, ${\underline{u}} \in C^{2,1}\left( \Omega \times (0,T]\right) \cap C\left( {\overline{\Omega }}\times [0,T]\right)$ with

$$\begin{aligned}{} & {} \frac{\partial {\overline{u}}}{\partial t}-\Delta {\overline{u}} \ge g({\overline{u}},{\varvec{x}},t), \end{aligned}$$

(11)

$$\begin{aligned}{} & {} \frac{\partial {\underline{u}}}{\partial t}-\Delta {\underline{u}} \le g({\underline{u}},{\varvec{x}},t), \end{aligned}$$

(12)

${\overline{u}}({\varvec{x}},0)\ge {\underline{u}}({\varvec{x}},0)$ on $\Omega$, and either ${\overline{u}}({\varvec{x}},t)\ge {\underline{u}}({\varvec{x}},t)$ or

$$\begin{aligned} \kappa {\overline{u}}+\frac{\partial {\overline{u}}}{\partial {\varvec{\eta }}}\ge \kappa {\underline{u}}+\frac{\partial {\underline{u}}}{\partial {\varvec{\eta }}} \end{aligned}$$

(13)

on $\partial \Omega \times (0,T]$, then either ${\overline{u}}\equiv {\underline{u}}$ or ${\overline{u}}>{\underline{u}}$ uniformly on $\Omega \times (0,T]$.

If the functions ${\overline{u}}$ and ${\underline{u}}$ satisfy the hypotheses of Lemma 2.2, then ${\overline{u}}$ is a supersolution and ${\underline{u}}$ is a subsolution of (1)–(3). The following lemma guarantees the existence of at least one stable periodic solution between any subsolution and supersolution to (1)–(3).

Lemma 2.3

(Existence of a stable periodic solution  [10, p. 67]) If strict sub- and supersolutions ${\underline{u}}<{\overline{u}}$ to (1)–(3) exist, then there exists a stable T-periodic solution $u^*\in ({\underline{u}},{\overline{u}})$ to (1)–(3) for all $x\in \Omega$ and all $t>0$.

This result implies that if $\mu _1<0$, then the function ${\underline{u}}=\epsilon \phi _1$ is a strict subsolution to (1)–(3) for any positive, sufficiently small constant $\epsilon$, while any constant larger than the carrying capacity $u\equiv 1$ is a strict supersolution to (1)–(3). It therefore follows from Lemmas 2.1 and 2.3 that there exists at least one stable T-periodic solution ${u^*\in ({\underline{u}},{\overline{u}})}$ to (1)–(3) for all $\varvec{ x}\in \Omega$ and all $t>0$ if $\mu _1<0$. It similarly follows from Lemma 2.1 that the trivial solution $u\equiv 0$ is asymptotically stable if $\mu _1>0$.

We employ a standard monotone iteration generation scheme to show that there, in fact, exists a unique, stable positive T-periodic solution to (1)–(3). By starting at a suitable initial iteration $u^{(0)}$, it is possible to construct a sequence of solutions $u^{(1)},u^{(2)},u^{(3)},\dots$ successively to the linear initial-boundary value problems

$$\begin{aligned}&\text{ D.E. } \quad {\displaystyle \frac{\partial u^{(n)}}{\partial t} - u^{(n)}\Delta = g(u^{(n-1)},\varvec{x},t)}, \quad \varvec{ x}\in \Omega \setminus \partial \Omega , \quad t > 0, \end{aligned}$$

(14)

$$\begin{aligned}&\text{ B.C. }\quad {\displaystyle \frac{\partial u^{(n)}}{\partial \varvec{\eta }}+ \kappa u^{(n)} = 0}, \qquad \qquad \qquad \qquad {\varvec{x}} \in \partial \Omega , \qquad \quad t \ge 0, \end{aligned}$$

(15)

$$\begin{aligned}&\text{ I.C. } \quad u^{(n)}(\varvec{ x},0) = f(\varvec{ x})\ge 0, \qquad \qquad \qquad ~ \varvec{ x}\in \Omega \end{aligned}$$

(16)

for all $n\in {\mathbb {N}}.$ Note that since $g(u,\varvec{x},t)$ is T-periodic for all $t>0$, it follows that $g(u^{(n-1)},\varvec{x},t)$ in (14) is also T-periodic for all $n\in {\mathbb {N}}.$ Denote the sequence corresponding to $u^{(0)}={\overline{u}}$ by $\{{\overline{u}}^{(n)}\}_{n\in {\mathbb {N}}_0}$, where ${\overline{u}}$ is a supersolution of (1)–(3), and call it an upper sequence. Similarly, denote the sequence corresponding to $u^{(0)}={\underline{u}}$ by $\{{\underline{u}}^{(n)}\}_{n\in {\mathbb {N}}_0}$, where ${\underline{u}}$ is a subsolution of (1)–(3), and call it a lower sequence. It can be shown that if, for some bounded functions ${\underline{c}}(\varvec{ x},t)$ and ${\overline{c}}(\varvec{ x},t)$, the function $g(u,\varvec{x},t)$ satisfies the condition

$$\begin{aligned} -{\underline{c}}(\varvec{ x},t)(u_1-u_2)\le g(u_1,\varvec{x},t)-g(u_2,\varvec{x},t)\le {\overline{c}}(\varvec{ x},t)(u_1-u_2) \end{aligned}$$

(17)

for all ${\underline{u}}\le u_2\le u_1\le {\overline{u}}$ and all $(\varvec{ x},t)\in (\Omega \setminus \partial \Omega )\times [0,T]$, then each of the two sequences converges monotonically to a unique solution of (1)–(3). This claim follows from the next four lemmas due to Pao [16]. The first lemma establishes the well-posedness of the lower and upper sequences defined above.

Lemma 2.4

(Well-posedness of the lower and upper sequences) If the condition in (17) is satisfied, then the two sequences $\{{\underline{u}}^{(n)}\}_{n\in {\mathbb {N}}_0}$ and $\{{\overline{u}}^{(n)}\}_{n\in {\mathbb {N}}_0}$ are well-defined, and both ${\overline{u}}^{(n)}$ and ${\underline{u}}^{(n)}$ are of class $C^\alpha (\Omega \setminus \partial \Omega \times [0,T])$ for each $n\in {\mathbb {N}}_0$.

Let $\overline{{{\mathcal {S}}}}$ denote the closure of the set ${{\mathcal {S}}}$. Then the following result guarantees the monotonicity of the lower and upper sequences.

Lemma 2.5

(Monotonicity of the lower and upper sequences) The lower and upper sequences $\{{\underline{u}}^{(n)}\}_{n\in {\mathbb {N}}_0}$ and $\{{\overline{u}}^{(n)}\}_{n\in {\mathbb {N}}_0}$ satisfy the monotone property

$$\begin{aligned} {\underline{u}} \le {\underline{u}}^{(n)}\le {\underline{u}}^{(n+1)}\le {\overline{u}}^{(n+1)}\le {\overline{u}}^{(n)}\le {\overline{u}} \end{aligned}$$

(18)

for all $(\varvec{ x},t)\in \overline{\Omega \setminus \partial \Omega \times [0,T]}$ and all $n\in {\mathbb {N}}_0$.

The next lemma guarantees that the lower and upper sequences converge to minimal and maximal solutions of (1)–(3), respectively.

Lemma 2.6

(Convergence of solutions to (1)–(3)) The limits $\lim _{n\rightarrow \infty }{\overline{u}}^{(n)}={\overline{u}}^*$ and $\lim _{n\rightarrow \infty }{\underline{u}}^{(n)}={\underline{u}}^*$ exist and satisfy the relation

$$\begin{aligned} {\underline{u}} \le {\underline{u}}^{(n)}\le {\underline{u}}^{(n+1)}\le {\underline{u}}^*\le {\overline{u}}^*\le {\overline{u}}^{(n+1)}\le {\overline{u}}^{(n)}\le {\overline{u}} \end{aligned}$$

(19)

for all $(\varvec{ x},t)\in \overline{\Omega \setminus \partial \Omega \times [0,T]}$ and all $n\in {\mathbb {N}}_0$.

It is next shown that the limits in Lemma 2.6 converge towards the unique T-periodic attractor of solutions to (1)–(3). This is accomplished by invoking the following positivity lemma.

Lemma 2.7

(Positivity of solutions  [15, p. 54]) Let c be a bounded function in $\overline{\Omega \times (0,T]}$. Suppose further that $w\in C(\overline{\Omega \times (0,T]})\cap C^{2,1}(\Omega \times (0,T])$ and that w satisfies

$$\begin{aligned}&\text{ D.E. } \quad {\displaystyle \frac{\partial w}{\partial t}} -\Delta w\ge 0, \qquad \qquad {\varvec{x}}\in \Omega \setminus \partial \Omega , \end{aligned}$$

(20)

$$\begin{aligned}&\text{ B.C. }\quad {\displaystyle \frac{\partial w}{\partial {\varvec{\eta }}}+\kappa w \ge 0}, \qquad \qquad ~{\varvec{x}} \in \partial \Omega , \end{aligned}$$

(21)

$$\begin{aligned}&\text{ I.C. }\quad {\displaystyle w({\varvec{x}},0)= f({\varvec{x}})\ge 0}, \quad \,\, {\varvec{x}}\in \Omega . \end{aligned}$$

(22)

Then $w(\varvec{x},t)\ge 0$ uniformly on $\Omega \times (0,T]$. Moreover, $w(\varvec{x},t)>0$ uniformly on $\Omega \times (0,T]$, unless $w(\varvec{x},t)\equiv 0$.

Based on the results above, we show that ${\underline{u}}^*=u^*={\overline{u}}^*$ following an argument of Pao [15, p. 55].

Lemma 2.8

(Nature of the lower and upper sequence limits) If the limits ${\underline{u}}^*$ and ${\overline{u}}^*$ are solutions to (1)–(3), then ${\underline{u}}^*=u^*={\overline{u}}^*$, which is the unique solution to (1)–(3) in the interval $u\in ({\underline{u}},{\overline{u}})$.

Proof

Let $w = {\underline{u}}^*-{\overline{u}}^*\le 0$. Then w satisfies the relation

$$\begin{aligned} \frac{\partial w}{\partial t}-\Delta w = g({\underline{u}}^*,\varvec{x},t)-g({\overline{u}}^*,\varvec{x},t)\ge {\overline{c}}({\overline{u}}^*-{\underline{u}}^*)={\overline{c}}w \end{aligned}$$

(23)

on $\Omega \setminus \partial \Omega \times [0,T]$ and the boundary condition

$$\begin{aligned} \frac{\partial w}{\partial \varvec{\eta }}+ \kappa w = 0, \quad w\in \partial \Omega \end{aligned}$$

(24)

for all $t\in [0,T]$ as well as the initial condition $w(0,\varvec{ x})=0$ on $\Omega$. It therefore follows from Lemma 2.7 that $w\ge 0$ uniformaly on $\Omega \times [0,\infty )$. This implies, by the definition of w, that ${\overline{u}}^*={\underline{u}}^*$. Suppose that $u^*$ is another solution to (1)–(3) in the interval $u^*\in ({\underline{u}},{\overline{u}})$. The functions $u^*$ and ${\underline{u}}$ are therefore ordered super- and subsolutions, respectively, of (1)–(3) and since the sequence $\{{\overline{u}}^{(n)}\}_{n\in {\mathbb {N}}_0}$ corresponding to ${\overline{u}}^{(0)}=u^*$ consists of the same function for any $n\in {\mathbb {N}}_0$, it follows that $u^*\ge {\underline{u}}^*$. Similarly, by considering $u^*$ and ${\overline{u}}$ as sub- and supersolutions, respectively, it follows by the same argument that $u^*\le {\overline{u}}^*$. Consequently ${\underline{u}}^*=u^*={\overline{u}}^*$. $\square$

The following result is a direct consequence of Lemmas 2.1–2.8.

Theorem 2.9

(The existence of a unique, positive periodic attractor) The model (1)–(3) admits a positive T-periodic solution $u^*(\varvec{ x},t)$ if and only if the principal eigenvalue $\mu _1$ of (7)–(9) is strictly negative and there exists a positive supersolution ${\overline{u}}$ to (1)–(3). Moreover, this positive T-periodic solution $u^*(\varvec{ x},t)$ is unique and globally asymptotically stable. If, however, the above eigenvalue is non-negative, then the trivial solution $u(\varvec{ x},t)\equiv 0$ is globally asymptotically stable.

3 Asymptotic analysis of model solution

If the principal eigenvalue $\mu = \mu _1$ to the eigenvalue problem (7)–(9) corresponding to (1)–(3) is negative, then it follows from Theorem 2.9 that there exists a unique positive T-periodic solution $u^*(\varvec{ x},t)$ which is globally attractive among all solutions to (1)–(3) for each positive choice of $\kappa$, provided that $\mu _1$ remains negative. In this section, the latter requirement is reformulated in terms of the existence of a maximally sustainable pro-rata harvest rate $\kappa ^*$.

Let ${\tilde{b}}(t)=\min _{\varvec{x}\in \Omega }b(\varvec{x},t)$ and consider the eigenvalue problem

$$\begin{aligned}&\frac{\partial \phi }{\partial t} - \Delta \phi -{\tilde{b}}(t)\phi ={\underline{\mu }}\phi ,\quad \,\varvec{ x}\in \Omega \setminus \partial \Omega , \end{aligned}$$

(25)

$$\begin{aligned}&\frac{\partial \phi }{\partial \varvec{\eta }}+\kappa \phi = 0,\qquad \qquad \qquad \quad \varvec{ x}\in \partial \Omega , \end{aligned}$$

(26)

$$\begin{aligned}&\phi T\text {-periodic in }t, \end{aligned}$$

(27)

with corresponding principal eigenvalue $\underline{\mu }=\underline{\mu _1}$ and since $b(\mathbf{{x}},t)\ge \tilde{b(t)}$ it follows that $\mu _1\le \underline{\mu _1}$ [9]. Moreover, let $\phi (\varvec{ x},t) = \alpha (t)\varphi _1(\varvec{ x})$, where $\varphi _1(\varvec{ x})$ is the principal eigenfunction of

$$\begin{aligned}&-\Delta \varphi = \gamma \varphi ,\quad ~~\varvec{ x}\in \Omega \setminus \partial \Omega , \end{aligned}$$

(28)

$$\begin{aligned}&\frac{\partial \varphi }{\partial \varvec{\eta }}+\kappa \phi = 0,\quad \varvec{ x}\in \partial \Omega , \end{aligned}$$

(29)

with corresponding principal eigenvalue $\gamma =\gamma _1$. Combining (7) and (28) yields

$$\begin{aligned} \varphi _1\frac{{\text{ d }}\alpha }{{\text{ d }} t} + \gamma _1\alpha (t)\varphi _1 - {\tilde{b}}(t)\alpha (t)\varphi _1 = {\underline{\mu }}\alpha (t)\varphi _1, \end{aligned}$$

(30)

from which it follows that

$$\begin{aligned} \alpha (t) = \alpha (0)e^{\int _{0}^{t}{\tilde{b}}(\tau )\,{\text {d}}\tau -\gamma _1t +\mu t} \end{aligned}$$

(31)

and since $\phi$ is T-periodic in t, the requirement

$$\begin{aligned} \int _{0}^{T}{\tilde{b}}(\tau )\,{\text{ d }}\tau -\gamma _1T +{\underline{\mu }}T =0 \end{aligned}$$

(32)

must hold, or equivalently

$$\begin{aligned} {\underline{\mu }} = \frac{1}{T}\left( \gamma _1T-\int _{0}^{T}{\tilde{b}}(\tau )\,{\text{ d }}\tau \right) . \end{aligned}$$

(33)

The problem with the expression for the principal eigenvalue $\mu =\mu _1$ in (33) is that it is expressed in terms of the the principal eigenvalue $\gamma =\gamma _1$ of (28)–(29). The following result due to Cantrell and Cosner [4], however, provides an estimate for $\gamma _1$ upon consideration of the variational formulation of (28)–(29).

Theorem 3.1

(The principal eigenvalue of (28)–(29)  [4, p. 98]) If $\partial \Omega$ is piecewise continuous of class $C^{2,\alpha }$, $\Omega$ satisfies the interior cone condition¹ and $\kappa \ge 0$, then the principal eigenvalue $\gamma =\gamma _1$ of (28)–(29) is

$$\begin{aligned} \gamma _{1}=\min _{\begin{array}{c} \varphi \in W^{1,2}(\Omega )\\ \varphi \not \equiv 0 \end{array}}\left[ \frac{\int _{\Omega } |\nabla \varphi |^{2}\, {\text{ d }}\varvec{ x} + \kappa \int _{\partial \Omega } \varphi ^{2}\,{\text{ d }} \varvec{s}}{\int _{\Omega } \varphi ^{2}\,{\text{ d }}\varvec{ x}}\right] . \end{aligned}$$

(34)

Moreover, the eigenfunction $\phi _1$ is positive on $\Omega$, $\gamma _1$ depends continuously on $\kappa$ with respect to $L^\infty (\partial \Omega )$ and $\gamma _1$ is simple, in the sense that the eigenspace associated with $\gamma _1$ is one-dimensional.

Upon substituting the test function $\varphi \equiv 1$ into (34), it follows that $\gamma _1\le \kappa |\partial \Omega |/|\Omega |$, where the term $|\partial \Omega |/|\Omega |$ is the perimeter-to-area ratio of $\Omega$. Consequently, it follows from (33) that

$$\begin{aligned} \underline{\mu _1} \le \frac{1}{T}\left( \frac{\kappa |\partial \Omega |}{|\Omega |}T-\int _{0}^{T}{\tilde{b}}(\tau )\,{\text{ d }}\tau \right) . \end{aligned}$$

(35)

Therefore, since $\mu \le {\underline{\mu }}$, it follows that $\mu _1<0$ if

$$\begin{aligned} \frac{1}{T}\left( \frac{\kappa |\partial \Omega |}{|\Omega |}T-\int _{0}^{T}{\tilde{b}}(\tau )\,{\text{ d }}\tau \right) <0, \end{aligned}$$

(36)

which may be rewritten as

$$\begin{aligned} \kappa < \frac{|\Omega |}{|\partial \Omega |T}\int _{0}^{T}{\tilde{b}}(\tau )\,{\text{ d }}\tau = \kappa ^*. \end{aligned}$$

(37)

Although the upper bound on $\mu _1$ is somewhat rough as a result of the choice of test function $\varphi \equiv 1$, it is clear that there exists a largest value of $\kappa$, denoted by $\kappa ^*$, for which $\mu _1 <0$. Theorem 2.9 consequently implies the existence of a unique positive periodic solution $u^*(\varvec{ x},t)$ to (1)–(3) which is globally asymptotically stable if $f(\varvec{ x})\not \equiv 0$ for all $\kappa <\kappa ^*$. Moreover, it follows from the comparison principle that $\mu _1$ is increasing in $\kappa$ and as a result there exists some other $\kappa ^{**}\ge \kappa ^*$ such $\mu _1 > 0$. Therefore, if $\kappa \ge \kappa ^{**}$, then the solution to (1)–(3) is attracted to the zero equilibrium. The analysis above is summarised in the following result, which contains a necessary condition on the pro-rata rate of boundary harvesting of the population in order for it to persist as $t\rightarrow \infty$.

Theorem 3.2

(Existence of a unique attractor of solutions to (1)–(3)) Let $u(\varvec{ x},t)$ be the unique solution to (1)–(3) corresponding to some initial condition $f(\varvec{ x})\not \equiv 0$. If $\kappa \in [0,\kappa ^*)$, then there exists a unique positive function $u^*(\varvec{ x},t)$, defined on $\Omega \times [0,\infty )$ and T-periodic in t, such that $\lim _{t\rightarrow \infty }|u(\varvec{ x},t)-u^*(\varvec{ x},t)|=0$ uniformly for all $\varvec{ x}\in \Omega$.

The result of Theorem 3.2 is confirmed by means of a numerical example in the special case of a single spatial dimension, where $\Omega =[0,1]$.

Example 1

(Numerical example for the case where $m=1$ and $\kappa >0$) Two numerical approximations of solutions to (1)–(3) in the special case of $m=1$ spatial dimension are illustrated graphically in Fig. 1, corresponding to two different choices of $\kappa$ whilst maintaining the same reaction term and initial condition in each example. The initial condition is chosen as $f( x)=0.5(1+ \sin (4\pi x))$ and it is assumed that $g(u,x,t) =u(1-u)(0.5+x-x^2)(0.5+0.49\sin (0.5t))$, from which it follows that $T = 4\pi$ and that $b(x,t) =(0.5+x-x^2)(0.5+0.49\sin (0.5t))$. In Fig. 1a, the value of $\kappa$ is chosen as $0.18>\kappa ^{**}$ and in Fig. 1b, the value of $\kappa$ is chosen as $0.16<\kappa ^*$. In the former case, the solution decays to the zero-equilibrium solution, as expected, while in the latter, the solution indeed converges to a $4\pi$-periodic attractor, as guaranteed by Theorem 3.2. $\blacksquare$

Fig. 1

The long-term asymptotic behaviour of solutions to (1)–(3) in the special case of a single spatial dimension with pro-rata harvesting on the spatial boundaries and with a $4\pi$-periodic reaction term $g(u,t) = u(1-u)(0.5+x-x^2)(0.5+0.49\sin (0.5t))$ and the non-zero initial condition $f( x)=0.5(1+ \sin (4\pi x))$

4 Properties of the periodic attractor

In Sect. 2 is was shown that the unique, positive solution to (1)–(3) is attracted to a unique, positive T-periodic solution $u^*(\varvec{x},t)$ if $\kappa <\kappa ^*$. The T-periodic attractor $u^*(\varvec{x},t)$ of solutions to (1)–(3) satisfies the boundary value problem

$$\begin{aligned}&\text{ D.E. } \quad {\displaystyle \frac{\partial u^*}{\partial t} = \Delta u^* + g(u^*,\varvec{x},t)}, \quad \,\, {\varvec{x}}\in \Omega \setminus \partial \Omega , \quad t \ge 0, \end{aligned}$$

(38)

$$\begin{aligned}&\text{ B.C. }\quad {\displaystyle \frac{\partial u^*}{\partial \varvec{\eta }}+ \kappa u^* = 0}, \qquad \qquad \qquad \,\, {\varvec{x}} \in \partial \Omega ,\quad \qquad \,t \ge 0, \end{aligned}$$

(39)

$$\begin{aligned}&\text{ P.C. } \quad g(u,\varvec{x},t)=g(u,\varvec{x},t+kT),\quad \qquad \qquad \qquad \quad t\ge 0,\quad k\in {\mathbb {N}}. \end{aligned}$$

(40)

It is shown in this section that the unique positive T-periodic attractor of solutions to (1)–(3) is a spatially symmetric function on $\Omega$ if $\Omega$ and g exhibit spatial symmetry.

Theorem 4.1

(Symmetry of the T-periodic attractor to (1)–(3)) Suppose $g(u^*,\varvec{x},t)$ is symmetric on $\Omega$ and that $\Omega$ is a symmetric domain in ${\mathbb {R}}^{m}$, in the sense that there exists an $(m-1)$-dimensional hyperplane $\Xi$ of symmetry for $\Omega$. Let $\varvec{s}: \Omega \mapsto \Omega$ be a symmetry mapping which maps any point $\varvec{ x} \in \Omega$ to its mirror image $\varvec{s}(\varvec{ x}) \in \Omega$ orthogonally across the hyperplane of symmetry $\Xi$. Then the unique positive periodic attractor $u^{*}$ to (1)–(3) is itself symmetric on $\Omega$ in the sense that it satisfies the relationship $u^{*}(\varvec{s}(\varvec{ x}),t)=u^{*}(\varvec{ x},t)$ for any $\varvec{ x} \in \Omega$.

Proof

The function $u^{*}(\varvec{ x},t)$ satisfies (38)–(40). Moreover, the mirror image $\varvec{s}(\partial \varvec{ x})$ of any boundary point $\partial \varvec{ x} \in \partial \Omega$ is again a boundary point of $\Omega$ and satisfies the relationship

$$\begin{aligned} {\displaystyle \frac{\partial u^{*}}{\partial \varvec{\eta }}\Big |_{\varvec{ x}=\partial \varvec{ x}}=\frac{\partial u^{*}}{\partial \varvec{\eta }}\Big | _{\varvec{ x}=\varvec{s}(\partial \varvec{ x})}}\text { for all }t\in [(k-1)T,kT]\text { and any }k\in {\mathbb {N}} \end{aligned}$$

by symmetry, while the mirror image $\varvec{s}(\varvec{ x})$ of any internal point $\varvec{ x} \in \Omega \setminus \partial \Omega$ is again an interior point of $\Omega$ and similarly satisfies the relationship

$$\begin{aligned} {\displaystyle \nabla u^{*}(\varvec{s}(\varvec{ x}))=-\nabla u^{*}(\varvec{ x})}\text { for all }t\in [(k-1)T,kT]\text { and }k\in {\mathbb {N}}. \end{aligned}$$

The symmetry mapping $\varvec{s}$ implies that the temporal derivative of $u^*$ remains unchanged after the transformation, in the sense that

$$\begin{aligned} \frac{\partial u^*({\varvec{x}},t)}{\partial t} =\frac{\partial u^*(\varvec{s}({\varvec{x}}),t)}{\partial t} \end{aligned}$$

for all $\varvec{ x} \in \Omega$ and all $t\in [(k-1)T,kT]$, with $k\in {\mathbb {N}}$. It therefore follows that

$$\begin{aligned}&\frac{\partial u^*}{\partial t} - \Delta u^{*}(\varvec{s}(\varvec{ x}),t)-g\left( u^{*}(\varvec{s}(\varvec{ x}),\varvec{s}(\varvec{ x}),t),t\right) \\&\quad = \frac{\partial u^*}{\partial t} -\nabla \cdot \nabla u^{*}(\varvec{s}(\varvec{ x}),\varvec{s}(\varvec{ x}),t)-g\left( u^{*}(\varvec{s}(\varvec{ x})),t\right) \\&\quad =\frac{\partial u^*}{\partial t} -(-1)^{2} \Delta u^{*}(\varvec{ x},t) -g\left( u^{*}(\varvec{ x}),\varvec{x},t\right) \\&\quad =0 \end{aligned}$$

for all $\varvec{ x} \in \Omega$ and all $t\in [(k-1)T,kT]$, with $k\in {\mathbb {N}}$, and that

$$\begin{aligned} \kappa u^{*}(\varvec{s}(\partial \varvec{ x}),t)=\kappa u^{*}(\partial \varvec{ x},t)=-\frac{\partial u^{*}}{\partial \varvec{\eta }}\Big |_{\varvec{ x}=\partial \varvec{ x}}=-\frac{\partial u^{*}}{\partial \varvec{\eta }}\Big |_{\varvec{ x}=\varvec{s}(\partial \varvec{ x})} \end{aligned}$$

for all $\varvec{ x} \in \partial \Omega$ and all $t\in [(k-1)T,kT]$, for any $k\in {\mathbb {N}}$. As a result, $u^{*}$ also satisfies the boundary value problem

$$\begin{aligned} \text{ D.E. }&\quad {\displaystyle \frac{\partial u^*}{\partial t} = \Delta u^*(\varvec{s}(\varvec{ x}),t) + g(u^*(\varvec{s}(\varvec{ x}),t),\varvec{s}(\varvec{ x}),t)},{} & {} \quad {\varvec{x}}\in \Omega \setminus \partial \Omega ,{} & {} \quad t \ge 0,\\ \text{ B.C. }&\quad {\displaystyle \frac{\partial u^*(\varvec{s}(\varvec{ x}),t)}{\partial \varvec{\eta }}+ \kappa u^*(\varvec{s}(\varvec{ x}),t) = 0},{} & {} \quad {\varvec{x}} \in \partial \Omega ,{} & {} \quad t \ge 0,\\ \text{ P.C. }&\quad g(u,t)=g(u,t+kT),{} & {} \quad u\ge 0,{} & {} \quad t\ge 0,{} & {} \quad k\in {\mathbb {N}}. \end{aligned}$$

This shows that $u^{*}$ is invariant under the transformation $\varvec{ x} \mapsto \varvec{s}(\varvec{ x})$. $\square$

The result of Theorem 4.1 is demonstrated by means of a numerical example for the special case of a spatial domain of dimension $m=1$.

Example 2

(Numerical examples of the symmetry of $u^*$ for $m=1$) Two numerical examples are provided in Fig. 2 corresponding to two different choices of the function $g(u^*,x,t)$, whilst keeping the value of $\kappa$ the same in each example. More specifically, the spatial domain $\Omega$ is assumed to have dimension $m=1$, taken as the unit interval [0, 1]. In this case, the 0-dimensional hyperplane $x=1/2$ is a hyperplane $\Xi$ of symmetry, and the function $s : \Omega \mapsto \Omega$, defined as $s( x)=1- x$, maps any point $x\in \Omega$ to its mirror image $1- x\in \Omega$ across $\Xi$. In Fig. 2a, the T-periodic attractor $u^*( x,t)$ is illustrated graphically for the symmetric choice $g(u^*,x,t)=u^*(1-u^*)(1-x+x^2)(0.5+0.49\sin (0.5t))$, while in Fig. 2b, the periodic attractor $u^*( x,t)$ is illustrated graphically for the asymmetric choice ${g(u^*,x,t)=u^*(e^1-e^{u^*})(1-x+x^2)(0.5+0.49\sin (0.5t))}$. It is evident from these examples, that $u^*$ is symmetric on the unit interval in both cases, as guaranteed by Theorem 4.1. $\blacksquare$

Fig. 2

An illustration of the positive $4\pi$-periodic attractor for $t\in [20\pi ,24\pi ]$ of solutions to (1)–(3) in the special case of a single spatial dimension and with the $4\pi$-periodic reaction term $g(u^*,x,t) = u^*(1-u^*)(1-x+x^2)(0.5+0.49\sin (0.5t))$ in (a) and the $4\pi$-periodic reaction term $g(u^*,x,t) = u^*(e^1-e^{u^*})(1-x+x^2)(0.5+0.49\sin (0.5t))$ in (b), for $\kappa = 0.1$ in both cases

The symmetry of the positive T-periodic attractor yields insight into the spatial locations of its maximum and minima. This is summarised formally in the following theorem.

Theorem 4.2

(Maximum of the T-periodic attractor to (1)–(3)) Let $u^*(\varvec{ x},t)$ be the unique positive T-periodic attractor of solutions to (1)–(3). If $g(u^*,\varvec{x},t)\ge 0$ for all $(u^*,\varvec{x},t)\in [0,1]\times \Omega \times [0,T]$, then $u^*(\varvec{ x},t)$ attains a global minimum on boundary points of $\Omega$ orthogonal to, and furthest away from, $\Xi$ and a global maximum at midpoints of $\Xi$ within $\Omega$ for any $\kappa \in [0,\kappa ^*)$.

Proof

It follows from Theorem 2.2 that ${\overline{u}}\equiv 1$ is a supersolution to (1)–(3) so that $u^*(\varvec{ x},t)<{\overline{u}}$ for any $\kappa \in (0,\kappa ^*)$. The condition $g(u^*,\varvec{x},t)\ge 0$ for all $(u^*,\varvec{x},t)\in [0,1]\times \Omega \times [0,T]$ implies, by (1), that $\partial u^*/\partial t -\Delta u^* \ge 0$ for all $\varvec{ x}\in \Omega \setminus \partial \Omega$. Therefore, it follows from the maximum principle that $u^*$ cannot attain a non-negative minimum M on $\Omega \setminus \partial \Omega \times (0,T]$, unless $u^*\equiv M$. A constant solution is, however, incompatible with the boundary condition in (39), unless $\kappa =0$, and so, the periodic attractor $u^*(\varvec{ x},t)$ is concave on $\Omega$, uniformly for all $t\ge 0$, and achieves a minimum on mirror image boundary points across $\Xi$ (that is, on points of $\partial \Omega$) furthest away from $\Xi$ and a maximum midway along $\Xi$ within $\Omega$. $\square$

It is shown next that the positive T-periodic attractor is a monotonically decreasing function of the pro-rata harvest rate $\kappa$.

Theorem 4.3

(Monotonicity of the T-periodic attractor to (1)–(3)) Let $\{u^*_n(\varvec{ x},t),\kappa _n\}_{n\in {\mathbb {N}}}$ be a bi-sequence of unique, positive T-periodic solutions to (1)–(3) corresponding to $u_n^*(\varvec{ x},t)=u^*(\varvec{ x},t)$ when $\kappa =\kappa _n$ for $n\in {\mathbb {N}}$, with $\kappa _1,\kappa _2,\kappa _3,\dots$ monotonically increasing such that $\kappa _n\rightarrow \kappa ^*$ as $n\rightarrow \infty$. Then $u^*_n(\varvec{ x},t)$ approaches some positive T-periodic solution $\varepsilon (\varvec{x},t)$ in the sense that $u^*_n(\varvec{ x},t) \rightarrow \varepsilon (\varvec{x},t)$ as $n\rightarrow \infty$, and $u^*_{n}(\varvec{ x},t)> u^*_{n+1}(\varvec{x},t)$ uniformly on $(\varvec{ x},t)\in \Omega \times [(k-1)T,kT]$ for all $n\in {\mathbb {N}}$ and any $k\in {\mathbb {N}}$.

Proof

Let $u^*_i(\varvec{ x},t)$ and $u^*_j(\varvec{ x},t)$ be the unique positive T-periodic attractors of solutions to (38)–(40) corresponding to $\kappa =\kappa _i$ and $\kappa =\kappa _j$, respectively, for some pair $i,j\in {\mathbb {N}}$ with $j>i$. Moreover, let $c(\varvec{ x},t)$ be a bounded function on $\Omega \times [0,\infty )$ and let $w = u^*_i - u^*_j$. It follows from (1) that

$$\begin{aligned} \frac{\partial u^*_i}{\partial t} -\frac{\partial u^*_j}{\partial t}-\Delta u^*_i + \Delta u^*_j = g(u^*_i,\varvec{x},t)-g(u^*_j,\varvec{x},t) \end{aligned}$$

(41)

and hence by the mean-value theorem that

$$\begin{aligned} \frac{\partial w}{\partial t}-\Delta w +cw = cw +\frac{\partial g}{\partial u}(\zeta ,\varvec{x},t)(u^*_i-u^*_j) \end{aligned}$$

(42)

for some intermediate value $\zeta =\zeta (\varvec{ x},t)$ between $u^*_i(\varvec{ x},t)$ and $u^*_j(\varvec{ x},t)$. Therefore,

$$\begin{aligned} \frac{\partial w}{\partial t}- \Delta w +cw = 0 \end{aligned}$$

by choosing

$$\begin{aligned} c(\varvec{ x},t) = -\frac{\partial g}{\partial u}(\zeta ,t). \end{aligned}$$

(43)

Moreover, it follows from the boundary condition in (39) that

$$\begin{aligned} \frac{\partial w}{\partial \varvec{\eta }}&= -\kappa _iu^*_i+\kappa _ju^*_j \nonumber \\&= -\alpha \kappa _ju^*_i+\kappa _ju^*_j \nonumber \\&= -\kappa _j(\alpha u^*_i-u^*_j) \nonumber \\&\ge -\kappa _j( u^*_i-u^*_j) \nonumber \\&= -\kappa _jw \end{aligned}$$

(44)

for some real number $\alpha \in (0,1)$. Finally, by choosing $u^*_i(\varvec{ x},0) = u^*_j(\varvec{ x},0)$, it follows from the condition in (40) that $w(\varvec{ x},T) = w(\varvec{ x},kT) \ge 0$ for any $k\in {\mathbb {N}}$. Therefore, the function w(x, t) satisfies the initial-boundary value problem

$$\begin{aligned} \text{ D.E. }&\quad \frac{\partial w}{\partial t} - \Delta w + cw =0 ,{} & {} \quad \varvec{ x}\in \Omega \setminus \partial \Omega ,{} & {} \quad t > 0, \end{aligned}$$

(45)

$$\begin{aligned} \text{ B.C. }&\quad {\displaystyle \frac{\partial w}{\partial \varvec{\eta }} +\kappa w \ge 0},{} & {} \quad \varvec{ x}\in \partial \Omega ,\end{aligned}$$

(46)

$$\begin{aligned} \text{ I.C. }&\quad w(\varvec{ x},T) = w(\varvec{ x},kT) \ge 0 ,{} & {} \quad \varvec{ x}\in \Omega ,{} & {} \quad k\in {\mathbb {N}}, \end{aligned}$$

(47)

and so it follows that $w>0$ uniformly on $\Omega$ for all $t>0$ by Lemma 2.7. The desired result therefore follows immediately. $\square$

The result of Theorem 4.3 is demonstrated by means of a numerical example for the special case of a single spatial dimension, where $\Omega =[0,1]$.

Example 3

(Monotonicity of the T-periodic attractor to (1)–(3)) In this example, the continuum of unique positive $4\pi$-periodic solutions $u^*$ to (1)–(3) for $\kappa \in (0,\kappa ^*)$ is considered in the case where $\Omega$ has dimension $m=1$, taken as $\Omega =[0,1]$, and $g(u^*,x,t) = u^*(1-u^*)(1-x+x^2)(0.5+0.49\sin (0.5t))$. Numerical approximations of the unique positive $4\pi$-periodic solutions $u^*$ to (1)–(3) are illustrated for four different values of $\kappa$ in Fig. 3. It is evident from the figure that the entire $4\pi$-periodic solution $u^*( x,t)$ decreases with respect to increasing $\kappa$, as guaranteed by Theorem 4.3. $\blacksquare$

Fig. 3

Numerical approximations of the positive $4\pi$-periodic solution $u^*(x,t)$ to (38)–(40) corresponding to four values of $\kappa$

5 An optimal harvesting yield

Theorem 4.3 implies that the positive T-periodic attractor of solutions to (1)–(3) is uniquely determined by the pro-rata harvest rate $\kappa$ and that the population density decreases everywhere in $\Omega \times [(1-k)T,kT]$ as the harvest rate increases, for any $k\in {\mathbb {N}}$. This leads to the following question: Does there exist a pro-rata harvest rate that maximises the total harvest over a single period in the long run and, if so, what is its value? The total harvest yield in the long run over a single period of length T along the boundary of the spatial domain is given by

$$\begin{aligned} z = \kappa \int _{0}^{T}\int _{\partial \Omega } u^*(\varvec{ x},t)\,{\text{ d }}\varvec{s}\,{\text{ d }}t. \end{aligned}$$

(48)

We therefore consider the PDE-constrained optimisation problem

$$\begin{aligned} \displaystyle {\max _\kappa }&\quad z=\kappa \int _{0}^{T}\int _{\partial \Omega } u^*(\varvec{ x},t)\,{\text{ d }}\varvec{s}\,{\text{ d }}t ,{} & {} {} & {} \end{aligned}$$

(49)

$$\begin{aligned} \text {s.t.}&\quad \displaystyle { \frac{\partial u^*}{\partial t} - \Delta u^* - g(u^*,t) =0},{} & {} \quad \varvec{ x}\in \partial \Omega ,{} & {} \quad t > 0,\end{aligned}$$

(50)

$$\begin{aligned}&\quad {\displaystyle \frac{\partial u^*}{\partial \varvec{\eta }}+ \kappa u^* = 0} ,{} & {} \quad \varvec{ x}\in \partial \Omega ,{} & {} \quad t>0,\end{aligned}$$

(51)

$$\begin{aligned}&\quad g(u,t)=g(u,t+kT) ,{} & {} \quad t\ge 0,{} & {} \quad u\ge 0,{} & {} \quad k\in {\mathbb {N}},\end{aligned}$$

(52)

$$\begin{aligned}&\quad \kappa \in (0,\kappa ^*).{} & {} {} & {} {}{} & {} \end{aligned}$$

(53)

We show that there exists at least one choice $\kappa ^o\in (0,\kappa ^*)$ which solves (49)–(53).

Theorem 5.1

(The existence of an optimal solution to (49)–(53)) There exists at least one optimal solution $\kappa ^o\in (0,\kappa ^*)$ to (49)–(53).

Proof

Let $h(\kappa )=\int _{0}^{T}\int _{\partial \Omega } u^*(\varvec{ x},t)\,{\text{ d }}\varvec{s}\,{\text{ d }}t$. The function $h(\kappa ):(0,\kappa ^*)\mapsto \mathbb {R^+}$ is a strictly decreasing function of $\kappa$, as a result of Theorem 4.3, in the sense that

$$\begin{aligned} \frac{{\text{ d }}h}{\text{ d }\kappa } <0. \end{aligned}$$

(54)

Moreover, if $\kappa =0$, it follows from (38)–(40) that $u^*\equiv 1$ uniformly on $\Omega$, showing that $h(0) = |\partial \Omega |\,T$. It furthermore follows from Theorem 3.2 that $\lim _{\kappa \rightarrow \kappa ^*}u^*(\varvec{ x},t)=0$, showing that $h(\kappa ^*) = 0$. Therefore, the function $z(\kappa )=\kappa h(\kappa )$ has two real roots, one at $\kappa = 0$ and the other at $\kappa =\kappa ^*$. Moreover, it follows from (49) that

$$\begin{aligned} \frac{{\text{ d }}z}{{\text{ d }}\kappa } = h(\kappa )+\kappa \frac{{\text{ d }}h}{{\text{ d }}\kappa }, \end{aligned}$$

(55)

so that

$$\begin{aligned} \frac{{\text{ d }}z}{{\text{ d }}\kappa }\Big |_{\kappa = 0}= |\partial \Omega |T >0 \end{aligned}$$

(56)

and

$$\begin{aligned} \frac{{\text{ d }}z}{{\text{ d }}\kappa }\Big |_{\kappa = \kappa ^*}= \kappa ^*\frac{{\text{ d }}h}{{\text{ d }}\kappa } <0. \end{aligned}$$

(57)

This implies, by the monotonicity of $h(\kappa )$ as a function of $\kappa$, that there exists at least one value $\kappa ^o\in (0,\kappa ^*)$ that maximises the yield function $z=z(\kappa )$. $\square$

Theorem 5.1 demonstrates that the objective function (49) is a curve parametrised by $\kappa$, starting at the point $(z, \kappa ) = (0,0)$ and initially increasing until it reaches the point $(z, \kappa ) = (z^o,\kappa ^o)$, where $z^o=\kappa ^oh^o$ and $h^o= \int _{0}^T\int _{\partial \Omega } u^o(\varvec{x},t)\,{\text{ d }}\varvec{s}\,{\text{ d }}t$ with $u^o(\varvec{x},t)$ being the value of the positive T-periodic attractor $u^*(\varvec{x},t)$ at $(\varvec{x},t)\in \Omega \times [(k-1)T,kT]$ corresponding to $\kappa ^o$. As $\kappa$ increases beyond $\kappa ^o$, the curve decreases again. Theorem 5.1 does not, of course, guarantee the existence of a unique solution $\kappa ^o$ to (49)–(53). The fact that the objective function $z(\kappa$) is increasing when $\kappa =0$ and is decreasing when $\kappa =\kappa ^*$ nevertheless implies that there exists at least one critical point of $z(\kappa )$ on the interval $(0,\kappa ^*)$. The result of Theorem 5.1 is illustrated by means of an example in the special case of $m=1$ spatial dimension.

Example 4

(The case of optimal harvesting in the case where $m=1$) In this example, the special case of (49)–(53) is considered where $\Omega$ has dimension $m=1$, taken as $\Omega =[0,1]$, and where $g(u^*,x,t)=u^*(1-u^*)(1-x+x^2)(0.5+0.49\sin (0.5t))$. In this case, the function $h(\kappa )$ can be approximated by computing successive numerical approximations of the function $u^*(x,t)$ satisfying (50)–(52) for discrete choices of $\kappa \in [0,\kappa ^*]$. The functional form of the objective function in (49) for this specific case is illustrated in Fig. 4.

Fig. 4

A numerical approximation of the objective function (49) for the parameters specified in Example 4

It is clear from the figure that (49) is unimodal and that the optimal solution $(\kappa ^o,u^o(x,t))$ to (49)–(53) is, in fact, unique in this case. For this reason, the well-known golden section search algorithm is employed to approximate $\kappa ^o\approx 0.105$ with a corresponding yield value $z^o \approx 0.605$. The output of the algorithm for a tolerance threshold value of 0.01 on the length of the interval of uncertainty is provided in Table 1. Finally, a numerical approximation of the optimal periodic solution $u^o(x,t)$ to (38)–(40), corresponding to the choice $\kappa =\kappa ^o$, is illustrated graphically in Fig. 5.

Table 1 The output of the golden section search algorithm when approximating the solution to the special case of (49)–(53) in Example 4

Fig. 5

A numerical approximation of the optimal solution $u^o(x,t)$ to (38)–(40) corresponding to $\kappa ^o=0.105$ for the parameters specified in Example 4

6 Conclusion and discussion

In this paper, the sustainability of pro-rata harvesting of a theoretical population on the boundary of a protected habitat was investigated. Our approach was to determine a maximally sustainable pro-rata harvest rate analytically which guarantees the persistence of the population modelled as an abstract initial-boundary value problem involving a time-periodic reaction-diffusion equation and spatial emigration on the boundary of a protected habitat in the form of Robin boundary conditions. As a result of studying this abstract model, which is an over simplification of the mechanisms governing the migration and reproduction of a biological population, a number of insights were gained into the long-term effectiveness of a protection zone in which harvesting is only allowed along the boundary.

It is acknowledged that, for several reasons, the model derived in (1)–(3) is a highly idealised abstraction of the growth and movement of real biological species. Therefore, before we discuss the aforementioned insights into the long-term effectiveness of establishing a protection zone for a biological species in more detail, we point out a few of the practical limitations of the modelling assumptions underlying the derivation of our abstract model:

• 1. The first important limiting assumption is that the population density of the biological species under consideration is modelled as a continuously differentiable function in both the space and time domains. In applications, a real biological species typically consists of finite numbers of discrete entities, in which case births and deaths may be noticeable as immediate and discrete discontinuities in population numbers over time. The assumption of continuity is, however, required in order to bring to bear on the model the considerable arsenal of analytic tools available within the mature realm of the differential calculus.

• 2. The second assumption pertains to the choice of the function $g(u,\varvec{x},t)$, which governs the intrinsic growth rate of the population which is assumed to be zero only when the population density is zero, but positive for any small, positive density. This assumption is highly idealised because in reality if the population density $u(\varvec{x},t)$ approaches zero in a specific region of the habitat at some time but not in other regions of the habitat at the same point in time, species members are still expected to be able to find mates for the purpose of reproduction. If, however, the species population density is, in fact, dwindling uniformly within its habitat, the species members may not be successful in finding mates while the population is still at a positive density. While the latter phenomenon may be modelled by assuming that the intrinsic birth rate of the species becomes negative at low population densities, a situation known as the allee effect [1], the model is nevertheless simplified by assuming a birth rate that is concave as a function of population density, because this renders the model amenable to analytic methods.

• 3. Another idealised assumption is related to the diffusion process describing the ease or effectiveness of the movement of species members, which is assumed to be constant over space and time as a result of the assumption of an isotropic spatial domain $\Omega$. It is acknowledged that in most applications, real biological species members may in certain cases be able to move more easily through certain regions of a heterogeneous habitat than through other regions, or that they may be able to move more easily through a particular region of their habitat at certain times of the year than during other times due to the effects of seasonality.

• 4. By not including the first-order differential operator $\partial u(\varvec{x},t)/\partial \varvec{x}$, sometimes referred to as the drift operator, in our model we are inherently assuming an unbiased random walk as the mechanism underlying the diffusion process, as opposed to a biased walk, in the derivation of (1)–(3). The latter case would result in the incorporation of the drift operator.

• 5. Finally, we point out that the pro-rata harvesting of the species is modelled by a constant Robin boundary condition, although in reality a spatially varying boundary may be more realistic, especially in the case where portions of the boundary may be lethal.

Nevertheless, the model is useful in the sense that it lends itself to a long-term asymptotic analysis from which practical insights may be gained. The mathematical analysis commenced in Sect. 2 with the establishment of a necessary and sufficient condition (in terms of the sign of the principal eigenvalue) for the existence of a unique positive T-periodic attractor of solutions to (1)–(3) corresponding to any non-zero initial condition $f(\varvec{x})$. This was followed by a long-term asymptotic analysis of the solution to (1)–(3), in which the necessary and sufficient condition of Sect. 2 for the existence of a unique positive, T-periodic attractor of solutions to (1)–(3) was rendered more practical in Sect. 3 by relating it to a maximal allowable pro-rata boundary harvest rate $\kappa ^*$. This is an interesting result from a substantiality point of view because it shows that a no-harvest protection zone may be rendered futile if the population is over exploited on the boundary of this protection zone.

Thereafter, we established a number of fundamental properties of the aforementioned periodic attractor in Sect. 4. These properties included establishing that the periodic attractor of solutions to (1)–(3) is spatially symmetric, provided that $\Omega$ and $g(u,\varvec{x},t)$ are also spatially symmetric, and that this periodic attractor is a monotonically decreasing function of the pro-rata harvest rate $\kappa$. The latter of these two results confirms the interesting and logical insight that increasing the pro-rata harvest rate on the boundary of a protected habitat will decrease the overall population density within the protection zone. This result is important for perhaps the most useful result of this paper where we established the existence of a pro-rata harvest rate in Sect. 5 which maximises the total harvest yield over a single interval of periodicity. In this case we showed that the pro-rata boundary harvest rate which maximises the total harvest yield over a single interval of periodicity is, in fact, always less than the maximal allowable pro-rata boundary harvest rate. This implies, in this case, that less is more because if the species were to be harvested at any sustainable pro-rata boundary harvest rate which is greater than the optimal harvest rate, then the population at the boundary will become small, resulting in inefficient harvesting. This will result in a less than an optimal harvest in the long run.