## Abstract

The aim of this paper is to provide a new mathematical model for a fishery by including a stock variable for the resource. This model takes the form of an infinite delay differential equation. It is mathematically studied and a bifurcation analysis of the steady states is fulfilled. Depending on the different parameters of the problem, we show that Hopf bifurcation may occur leading to oscillating behaviours of the system. The mathematical results are finally discussed.

## 1. Introduction

Bioeconomics mainly considers the dynamics of renewable resources from the point of view of control theory (Clark 1990) focusing on the existence of a maximum sustainable yield. Most models take into account two state variables, the resource with mass in pounds or mass per unit of area and the number of firms involved in the exploitation of the resource. In the context of fisheries, the variables are the fish density and the fishing effort. We refer to earlier contributions presenting time-continuous models (Smith 1968, 1969) as well as time-discrete versions (Barbier *et al.* 2002). These mathematical models consider two equations. The first equation represents the variation with time of the resource density which is due to its natural growth and to harvesting. A second equation describes the time variation of the number of firms which is assumed to increase when the exploitation is rentable and conversely. We also refer to some fishery models (Mchich *et al.*2002, 2006), where the resource was assumed to grow logistically and to be harvested according to a Lotka–Volterra type I functional response which is also commonly referred to as a Schaefer function. In these models, it is commonly assumed that the number of firms (or the fishing effort) varies according to the difference between the benefit and the cost of the exploitation at the present time *t*. If it is rentable, more firms are entering in the resource exploitation and conversely. In this paper, we are going to assume that part of the harvested resource is storaged. For example, one can think of a fish resource that can enter into a frozen fish stock for some time before being reinjected in the market. Therefore, we shall add a third equation to a classical fishery model governing the time evolution of a fish stock variable. The aim of this paper is to investigate the effects of such a fish stock on the dynamics of the fishery. As we shall see, from the point of view of modelling, this leads to delay differential equations.

The paper is organized as follows. In §2, we present the mathematical delay differential model presented in the context of a fishery. In §3, we present the stability and bifurcation analysis of the model. Section 4 is devoted to the interpretation and to the discussion of the results and to the discussion of several perspectives for this work.

## 2. A delay differential equation model

The starting point of our modelling is the simple fishery models like those given in Barbier *et al.* (2002) and Mchich *et al.* (2002, 2006). These models consider two state variables: the density of the resource *n*(*t*) and the fishing effort *E*(*t*). Such models read as
2.1
Here the resource *n* is assumed to have a logistic growth. For simplicity, we set the carrying capacity to 1. Function *φ*(*n*,*E*) denotes the capture per unit of time. It usually takes the form *φ*(*n*,*E*)=*q**n**E*, while the parameter *p* denotes the constant price. *p**φ*(*n*,*E*) represents the instantaneous benefit. The term −*c**E* denotes the monetary cost in terms of salaries, fuel, taxes, etc. *c* is the cost per unit of fishing effort and unit of time. Thus, the second equation assumes that the fishing effort is going to increase when the fishery makes some profit and conversely. In this simple model, it is assumed that all the catch is immediately sold and contributes to the fishing effort *E*. However and more generally, one can assume that some part of the instantaneous catch is immediately sold while the rest enters into some stock in order to be manufactured for instance (frozen, smoked, and so on). After some time, the fish stock returns to the fish market to be sold. For that purpose, we generalize the previous simple model by introducing a new stock variable, denoted by *S*. Some fraction *η*∈[0,1] of the fish catch will enter the stock compartment while the other fraction 1−*η* will immediately be sold to contribute to the fishing effort *E*. Under these assumptions, we now consider the following new model:
2.2
Here the parameter *δ*>0 denotes the rate at which the fish stock returns to the fish market and contributes with price *p* to the fishing effort. For simplicity, we assume that the price is the same for the fish which is immediately sold and the fish coming from the stock. Under some suitable conditions about the initial condition, we can eliminate the variable *S* by integrating the third equation. Therefore, we can reduce the above system of ODEs leading to the following two components of delay differential equations:
2.3

In this work, we study the following delay differential equation:
2.4
where *n*_{t}(*θ*)=*n*(*t*+*θ*) and *E*_{t}(*θ*)=*E*(*t*+*θ*) denote the history functions. This problem is also supplemented by some initial conditions. In system (2.4), the kernel *w*(*θ*) takes into account the management of the stock as well as some delay to manufacture the resource.

The dynamics of system (2.2) with *φ*(*n*,*E*)=*q**n**E* is rather simple. This system can have up to three stationary states. The points (*n*,*E*)=(0,0) and (*n*,*E*)=(1,0) always exist. The first one is always unstable while the second one is stable if and only if *c*/*p**q*>1. When *c*/*p**q* crosses the critical value 1, the point (*n*,*E*)=(1,0) undergoes a transcritical bifurcation leading to the appearance of another stationary point that becomes stable when *c*/*p**q*<1.

In the model we consider here, that is, equation (2.4), the dynamics is rather different. The value *c*/*p**q*=1 is also a bifurcation point and when this parameter crosses 1, the fishery free stationary point (*n*,*E*)=(1,0) also undergoes a transcritical bifurcation. The resulting stationary point is not always stable anymore. This point can also undergo some Hopf bifurcation leading to the appearance of periodic dynamics. The next section will shed some light on this phenomenom by studying the dynamics of this system with respect to the two key parameters *η* and *δ* when the weighted function is given by the simple exponential damped rate *w*(*θ*)=*δ*e^{δθ} with .

## 3. Stability analysis and bifurcations

We consider the delay differential equation described in §2 that reads
3.1
Here the parameters *r*, *q*, *p* and *c* are all non-negative and *η*∈[0,1]. Function is a weighted function taking into account the past and such that . Moreover we assume that there exists *μ*_{0}>0 and some constant *C*>0 such that
3.2

Moreover, we have denoted by *n*_{t} and *E*_{t} the history functions. This means that are defined by *n*_{t}(*θ*)=*n*(*t*+*θ*) and *E*_{t}(*θ*)=*E*(*t*+*θ*).

This problem is supplemented with initial conditions *n*_{0} and *E*_{0} such that
3.3

Then we have the following result.

## Proposition 3.1.

*For each bounded and continuous initial conditions n*_{0} *and E*_{0} *satisfying equation (3.3), then problem (3.1) has a globally defined solution (n(t),E(t)) such that*

The proof of this result is classical (e.g. Hale & Kato 1978; Diekmann *et al.* 1995; Adimy *et al.*2004*a*,*b*, submitted and references therein) and left to the reader.

We are now interested in the asymptotic behaviour of this model. Problem (3.1) can have up to three stationary points. Since we assume that , the stationary points of equation (3.1) are given by the resolution of the system 3.4

The existence of these feasible stationary points is summarized in the following lemma.

## Lemma 3.2.

*When c/pq≥*1 *system (3.4) only has the two feasible boundary solutions*
*When c/pq<*1, *then system (3.4) has three stationary solutions, the two boundary stationary points (n*_{1}*,E*_{1}*) and (n*_{2}*,E*_{2}*) and the following interior solution:*
3.5

We now investigate the stability of these different stationary points. Owing to assumption (3.2), a suitable functional setting can be developed following the method introduced by Liu *et al*. (2008). Some details are given in appendix A where it is shown that the stability analysis of a stationary point (*n**,*E**) is related to its characteristic equation given by
3.6
where we have set

We first give some results about the stability of the boundary stationary points. We have the following result.

## Proposition 3.3.

*The trivial point (n*_{1}*,E*_{1}*) is always unstable. When c/pq≥*1, *then we have*
*In addition, if n*_{0}*≠*0 *then n(t)→*1 *when* . *When c/pq<*1, *then the boundary equlibrium (n*_{2}*,E*_{2}*) is unstable.*

## Proof.

We first consider the case of (*n*_{1},*E*_{1}). With this point, the characteristic equation reads
Thus, *λ*=*r*>0 is a solution and the result follows.

In the case of (*n*_{2},*E*_{2}), the characterisic equation reads
When *c*/*p**q*<1, we consider the map defined by
It is a decreasing map from *p**q*−*c*>0 to *p**q*(1−*η*)−*c*. Thus, the equation *φ*(*λ*)=*λ* has a solution *λ*_{0}>0 and the result follows.

Let us now consider the case where *c*/*p**q*>1. To show that *E*(*t*)→0 when , let us consider the functional defined by
Let us show that the map *t*→*V* (*E*_{t}) is decreasing. First we have
Therefore, we get

Since *n*(*t*)≤1 for every time, we obtain that
Since *p**q*−*c*<0, the map *V* is a Lyapunov functional and this proves that *E*(*t*)→0 as . The proof of the proposition is completed. ▪

We consider the stability of the interior equilibrium (*n*_{3},*E*_{3}). Its stability relies on the properties of the characteristic equation. With a general weighted function *w*, the characteristic equation is hard to study. For that purpose, we concentrate on a particular weight *w* given by
3.7
for some given *δ*>0. This assumption corresponds to the use of the stock with a constant rate *δ*>0.

Then we will prove the following result.

## Proposition 3.4.

*Let us assume that c/pq>*1 *and w(x)=δe*^{δx} *for some δ>*0. *Then the interior point (n*_{3}*,E*_{3}*) is locally asympotically stable if and only if*
3.8

## Proof.

Since *c*/*p**q*>1, the point (*n*_{3},*E*_{3}) is a feasible equilibrium. Then the stability of this point is related to the characteristic equation that reads
Since we assume that the weighted function *w* is given by equation (3.7), we obtain that
Thus, the characteristic equation is equivalent to the following polynomial equation:
The stability conditions now follow from the Routh–Hurwitz criterion and we get that (*n*_{3},*E*_{3}) is stable if and only if
This completes the proof of the proposition. ▪

## Remark.

Let us notice that when , system (3.1) with *w(θ)=δe*^{δθ} formally converges to sytem (2.2) with *φ(n,E)=qnE*. So that condition (3.8) is satisfied when *δ* is large enough.

In the same way, when *η*=0, then system (3.1) is equivalent to system (2.2) with *φ(n,E)=qnE* and we know that (*n*_{3},*E*_{3}) is locally asymptotically stable. On the other hand, when *η*=0, condition (3.8) is rewritten as
which is a condition obviously true. Therefore, proposition 3.4 recovers the well-known results for system (2.2) with *φ(n,E)=qnE.*

We now study the stability boundary (*δ*+*r**n*_{3}+*η**c*)(*n*_{3}*δ*+*n*_{3}*c**η*+*c*(1−*η*)(1−*n*_{3}))=*c**δ*(1−*n*_{3}). We will prove in the following result that this boundary stability corresponds to Hopf bifurcation for the dynamical system (e.g. Hassard *et al.* 1981; Diekmann *et al.* 1995; Magal & Ruan 2009) and references therein concerning Hopf bifurcation). To do that we will use the intensity of the delay *δ*>0 as a bifurcation parameter. We will prove the following result.

## Theorem 3.5.

*Let us assume that c/pq<*1. *Then let us assume that there exists δ*_{0}*>*0 *satisfying the boundary instability condition*
3.9
*Then δ=δ*_{0} *is a Hopf bifurcation point of system (3.1) around the equilibrium (n*_{3}*,E*_{3}*).*

## Proof.

Let us look for purely imaginary eigenvalue *λ*=*i**ω*. Using the characteristic equation leads to the equation
Thus, we get
and
Owing to assumption (3.9), we obtain that there exists a unique *ω*_{0}>0 such that
Let us denote by the characteristic equation, that is,
After simple computations, we obtain
By using the implicit function theorem, we know that there exists *ε*>0 and a *C*^{1}-map such that
3.10
Moreover, we have
3.11
To complete the proof of the theorem, it remains to check the transversality condition; this means . Since , this transversality condition is equivalent to show that
To show that, note that we have
Therefore, we obtain that
Thus, the transversality condition follows and completes the proof of the result. ▪

## Remark.

We can easily check that when (3.9) is satisfied for some *δ*_{0}>0, it is also satisfied for another value of *δ*_{1}>0. Therefore, when increasing the value *δ*>0, we first obtain a Hopf bifurcation of the coexistence stationary state (*n*_{3},*E*_{3}). When going on to increase the value of *δ*, we observe an oscillating regime and then the coexistence stationary state restabilizes. We can also notice that the condition for the Hopf bifurcation (3.9) is in particular sastified when *n*_{3} is small enough. This means that when the density of the resource is small enough, the management of the stock, described here by the parameter *δ*>0, is important and drives the dynamics of the system.

The parameter *δ* appears to be an important bifurcation parameter. However, we could also study in the same way the influence of the parameter *η*∈(0,1). It can be proved by using the same kinds of computations that when *η*∈(0,1) satisfies the boundary instability (3.9), this also leads to a Hopf bifurcation in the dynamics. The detailed computations are omitted here. Some numerical simulations will be provided in the next section in order to consider this phenomenom.

## 4. Discussion

This section is devoted to some numerical simulations for problem (3.1).

We observe from the bifurcation diagram in figure 1 that when *η* increases, this means the proportion of fished resource that directly enters the stock, then the range of instability for the delay *δ* is also increasing. This essentially means that the stock acts in favour of oscillating dynamics.

figure 1 shows that for a small proportion of storage, *η*<0.35, the fishery is always stable. Above this threshold value, instability can occur in a certain range of the parameter *δ*. Let us fix *η* above this threshold, and figure 1 shows that the fishery is stable for small as well as large *δ* values. Indeed, when *δ* is small enough, the storaged fish remains a very long time in the stock, almost never returns to the market or at a very small rate and therefore has no real effect on the fishery which remains stable. When *δ* is large enough, the fish immediately returns to the market after being storaged for a very short time and this also has no real effect on the stability of the fishery. These later cases are thus almost similar to stock absence. For intermediate *δ* values, instability occurs. This instability makes sense. Indeed, let us consider some time when the fish density is small, thus the instantaneous catch is also small, and as a consequence without any stock, the fishing effort would reduce fast because the difference between the benefit due to the actual catch and the costs would be large. However, in the case of a fish stock, fish harvested in previous times (when the resource was much more abundant) return at a constant rate *δ* to the market and this in turn supports the fishing effort. We can understand that this can induce a destabilizing effect on the fishery because in that case the fishing effort is maintained at a rather large level with respect to the actual small abundance of the fish. This in turn will provoke a more important decrease in the fish density, which is overexploited and this destabilizes the fishery. A similar and converse process can occur when the fish density is abundant.

We now take a more detailed look at the dynamics of the model by showing the trajectories of the system. The results are shown in figures 2 and 3. In figure 2, we fix the values of the parameters and only vary the delay *δ*. We obtain the expected dynamics provided by figure 1. Indeed, with a small delay (figure 2*a*), the solution stabilizes towards the positive equilibrium, when the delay increases (figure 2*b*), the trajectory oscillates, while when going on increasing the delay (figure 2*c*), the trajectory restabilizes.

In figure 3, we fix the parameter set and we vary the proportion *η*. As expected in the bifurcation diagram, we observe that increasing *η* induces some oscillations in the dynamics. We can also notice that the amplitudes of oscillation increase with respect to *η*. In figure 3*c*, we obtain about 0.75 of the amplitude for the density of fish.

We complete this numerical study of instability of the coexistence equilibrium by looking at the influence of the price variable. The results are summarized in figure 4. We observe that increasing the price leads to an increasing instability area. We can understand this because in the case of fixed fishery costs, when the price of the fish is high, the benefit of the fishery is also large. As the variation of the fishing effort depends on the difference between the benefit and the costs, a high price maintains the fishing effort at a larger level than it should be in the case of a lower price. Therefore, large prices lead to overexploitation of the resource with respect to fishery with low market prices. Overexploitation in general promotes instability and this explains why the domain of instability grows as the market price increases as shown in figure 4. More precisely, when the price increases, the resource density component of the coexistence stationary state is smaller and smaller and therefore, due to proposition 3.4 and theorem 3.5, we obtain that this acts in favour of destabilization and oscillating dynamics.

Finally, we numerically confirm in figure 5 the simple dynamics proved in proposition 3.3. Here the parameter set is chosen such that *c*/*p**q*>1. The fishing effort rapidly goes to zero while the density of fish recovers to its maximum value of 1.

We now come back to model (3.4) in order to show how to improve it. These improvements will be considered in some forthcoming works. The first improvement we could consider consists of the following system:
4.1
In this model, we add a term in the second equation in order to consider the storage cost. We can also rewrite this problem in the form of delay differential equations. This kind of problem leads to some hard questions. Indeed, the trajectories of the dynamical system do not necessarily stay positive. Therefore, we enter the question of the viability of the fishery model that has been addressed by Aubin (1991; see also Cury *et al.* 2005 and references therein). There are a lot of questions related to this topic, in particular, to describe the viability kernel. These kinds of questions will be addressed in a forthcoming work.

In this contribution, we have limited our study to the simplest case of a fishery with a stock. It appears that such a simple model already exhibits interesting dynamics. However, several improvements of this work could be considered in future studies. For example, several fish species that are commercialized could become scarce, their density being maintained at very low levels with important risks of extinction. Therefore, it would be realistic for such overexploited species to consider that the fish stock does not grow logisticaly but according to an Allee effect dynamics.

On the other hand, instead of considering a fixed market price, one could also consider that the price of the resource is not constant but varies as the result of the supply (the catch) and some demand function depending on the fish market price.

We could also consider the effects of economic feedbacks on the management of stock (the example of Japanese tuna) which is storaged (frozen) when it is abundant and unstoraged when the actual catch is small. In that case, the aim of the stock is to maintain the supply at a rather stable level in order to maintain the market price at some desired level. We intend to investigate the effects of these different factors on the fishery dynamics in future studies.

## Appendix A

In this appendix, we provide a mathematical functional framework in order to justify the compuations given in §3 to fulfil the bifurcation analysis. Here we use the methodology introduced by Liu *et al.* (2008). The aim of this appendix is to consider a linear retarded equation, to rewrite it as an abstract Cauchy problem and finally to show the applicability of the results of Magal & Ruan (2009) concerning the bifurcation analysis by checking that the essential growth rate for some operator is negative. To do so, we introduce a suitable functional framework. Let *γ*>0 be given such that 0<*γ*<*μ*_{0}, where *μ*_{0} is defined in equation (3.2). Then we consider the Banach space *C*_{b} of continuous and bounded maps from into and we introduce the Banach space
which is endowed with the norm
Then the corresponding linear retarded equation considered in this paper reads as
A1
where *B* is some *N*×*N* real matrix while operator has the form
wherein *D* is some given *N*×*N* real matrix. In order to tackle this retarded equation we shall follow the method introduced by Liu *et al.* (2008). Let us rewrite equation (A1) as an abstract Cauchy problem with a non-densely defined operator. For that purpose, we define a map *u*≡*u*(*t*,*θ*) on by *u*(*t*,*θ*)=*x*(*t*+*θ*). Next one can check that function *u* satisfies the following partial differential equations:
A2
This formulation allows us to study equation (A1) as a non-densely defined Cauchy problem in the Banach space
endowed with the usual product norm. We define the operator by
wherein
We also define *L*:*X*_{0}→*X* with and
Then by setting (as in Liu *et al.* 2008)
and identifying the two functions *v*(*t*) and *u*(*t*), one obtains an equivalent formulation for equation (A2) as an abstract Cauchy problem
A3

Then we shall show the following result that allows us to justify the computation for bifurcation given in this paper (see Magal & Ruan 2009).

### Theorem A.1.

*The linear operator* *is a Hille–Yosida operator. Moreover, the part (A*+*L*)_{0} *in X*_{0} *is the infinitesimal generator of a strongly continuous semigroup {T*_{(A+L)0}(*t*)}_{t≥0} *on X*_{0}. *Moreover, the essential growth rate ω*_{ess}((*A*+*L*)_{0}) *satisfies*

The proof of this result requires several steps. Here we only show the main steps of the proof and skip the details. We refer to Liu *et al.* (2008) for more details (we also refer to Ducrot *et al.* (in press, submitted) for details, definitions and similar methodology). The first step is given in the following lemma.

### Lemma A.2.

*The linear operator* *is a Hille–Yosida operator.*

Indeed, one can easily check that and The result follows from simple computations.

We now consider the part *A*_{0} of *A* in *X*_{0} that is defined by
and
From this definition, we obtain that *A*_{0} is defined by
From lemma A.2, we know that *A*_{0} is the infinitesimal generator of a *C*_{0}-semigroup on *X*_{0}. Moreover, one can derive an explicit expression for the semigroup {*T*_{A0}(*t*)}_{t≥0}
where *T**φ*=*φ*(0) while
Let us now notice that the semigroup {*T*_{A0}(*t*)}_{t≥0} is quasi-compact since *T* is a finite-rank operator and thus compact. Next for each *t*≥0, we have
Thus, we get that *ω*_{ess}(*A*_{0})≤−*γ*. Now since operator *L*:*X*_{0}→*X* is bounded (since *γ*<*μ*_{0}) and compact (because of finite rank) one can conclude that is a Hille–Yosida operator and
Thus, by using the result of Ducrot *et al.* (2008), one gets that
and the result is proved.

## Footnotes

One contribution of 17 to a Theme Issue ‘From biological and clinical experiments to mathematical models’.

- © 2009 The Royal Society