Hybrid dynamics of multi-species resource exploitation

Abstract

In this paper, we analyze a bio-economic model of exploitation of renewable commercial resources. To take into account the typically continuous-time modeling of biological species and, instead, of the specialized harvesting activities, which by its nature cannot change continuously, the resulting dynamic system is of the hybrid type, i.e. continuous for biological variables and discrete for the economic ones. Through a discretization of the continuous variables, the problem is then reformulated by means of a three-dimensional map. Of this map, we study the dynamic properties to understand how economic parameters influence the long-run availability of resources.

1 Introduction

The balance between the need to conserve resources and their use for economic purposes represents a challenging problem. To avoid over-exploitation that leads to the “tragedy of the commons” (see Hardin (1968)), forms of harvesting regulation are usually employed, mainly based on optimal control techniques.¹ In spite of the various forms of regulation employed, most natural resources are over-exploited, with many species being extinct or on the verge of extinction, see, for instance, McWhinnie (2009).

One of the main reasons for this failure relates to the difficulty in monitoring and enforcing optimal rules. For this reason, one could replace those forms of rigid control, difficult to impose, with milder forms of regulation. An example relating to the possible dynamics of these rules is provided in Bischi et al. (2013a) and concerns the harvesting of two species of shellfish from the Adriatic Sea (Venerupis aurea and Callista chione). These species do not interact biologically with each other and, on the basis of a proposal to regulate their exploitation, fishers must only commit themselves to harvest only one of the two resources for a given period of time. This can occur because an authority imposes that fishers have to stick to the decided strategy for a given period of time after which they can reconsider their decisions on the basis of observed profits. The dynamics of the resources, modeled by means of a continuous-time system, are therefore coupled with a discrete-time dynamics on the species to be taken, based on the profits from fishing. The model that follows is then a hybrid system that combines continuous and discrete time dynamics, see Aubin et al. (2002), Goebel et al. (2009) and Haddad et al. (2006). From a mathematical point of view, even if hybrid models are more complicated to study, they are closer to a more correct description of the reality being analyzed (see also Bischi et al. (2013b), Lamantia and Radi (2015) and Bischi et al. (2015) for related models).

In this paper, we provide a dynamic analysis of the hybrid dynamical system proposed in Bischi et al. (2013a). In addition to the reasons indicated in Bischi et al. (2013a), this type of modeling is coherent to all those situations of specialized fisheries, in which a particular species is targeted and for a certain period of time, due to technical requirements of the harvesting (types of gears) or rules, it is not possible to change the species taken, see French et al. (2019) for an analogous example. The purpose of this paper is twofold. On the one hand, we extend the model to include financial evaluation of the interests on profits in the various harvesting periods. In addition, we provide an equivalent reformulation of the model in Bischi et al. (2013a) by means of a discretization of the continuous-time dynamics of the resources in the periods in which there may be a change in the choices of the species to be taken.² These choices are then modeled by means of an evolutionary dynamic (discrete-time replicator type), which is based on the profits obtained over the previous periods (see Weibull (1995) and Hofbauer and Sigmund (1998)). The discrete-time representation of the hybrid dynamical system proposed in Bischi et al. (2013a) allows us to underline interesting aspects of the global dynamics of the model that were neglected in Bischi et al. (2013a). The first relevant aspect that we are able to underline is the type of bifurcation, that is period-doubling or flip, through which the inner equilibrium of the system loses stability when a different time scale between the decision process of quantity to harvest and the one of the species to harvest is introduced. A second relevant aspect that we underline by means of numerically generated bifurcation diagrams is the existence of chaotic dynamics. Specifically, we show a sequence of period-doubling bifurcations leading to chaotic dynamics. A third relevant aspect that we underline is the presence of coexisting attractors, that introduces path dependence and the related economic issues. Therefore, the current work adds to the existing literature by explaining for the first time the mechanisms behind the global dynamics of the hybrid bio-economic model proposed in Bischi et al. (2013a).

The structure of the paper is as follows. In Sect. 2, we summarize the essential elements of the bio-economic model considered in Bischi et al. (2013a). In Sect. 3, we reformulate that hybrid model through a three-dimensional map and propose a possible extension to take into account the compounding of risk-free interests in harvesting profits. In Sect. 4, we present the main results on the global dynamics of the model, in particular, regarding the influence of economic parameters on the dynamics of the species subject to exploitation. Section 5 concludes. All the proofs are in Appendix A.

2 The bio-economic setup

Here we briefly recap the model proposed in Bischi et al. (2013a) for the particular case of interest for this paper and present a possible extension of the model. Consider two non-interacting species subject to commercial harvesting and denote their biomass measures by $X_{1}$ and $X_{2}$. Biomass i follows logistic natural growth (see Melià and Gatto (2005)) and its time evolution can be described by an ordinary differential equation (ODE) of the type

$$\begin{aligned} {\dot{X}}_{i}\left( t\right) =X_{i}\left( t\right) \rho _{i}\left( 1-\frac{X_{i}\left( t\right) }{k_{i}}\right) -H_{i}\left( X_{i}\left( t\right) ,r_{i}\left( t\right) \right) \text {, } \quad i=1,2, \end{aligned}$$

(1)

where $r_{i}$ is the share of fishers that harvest species i while $\rho _{i}>0$ and $k_{i}>0$ are, respectively, the intrinsic rate of growth and the carrying capacity of species i. In (1), the term $H_{i}$ denotes the instantaneous harvesting of species $i=1,2$.

The model in Bischi et al. (2013a) describes the exploitation of specialized fisheries, in the sense that $N\in \ {\mathbb {N}} \setminus \{0,1 \}$ fishers are present but, given the specialized nature of the fishery, for a certain period of time each fisher can only harvest one species at once. This can be due to technical constraints linked to the adopted fishing technique (e.g. types of gears) or to the legislation (i.e. authority restricting each agent to catch only one species at a time). As a consequence, it is not possible to carry out a continuous change of the species caught but the targeted species can change only at the beginning of a fishing period.

For the purpose of this paper, let us consider here the particular case of the model in Bischi et al. (2013a) with constant selling price $a_{i}>0$ for species i. Thus, if species i is targeted, the choice on quantity $h_{i}$ to harvest is carried out by maximizing profits of the form

$$\begin{aligned} \pi _{i}\left( X_{i}\left( t\right) ,h_{i}\right) =a_{i}h_{i}-\gamma _{i}\frac{h_{i}^{2}}{X_{i}\left( t\right) }\text {, }\quad i=1,2, \end{aligned}$$

(2)

where the term

$$\begin{aligned} C_{i}\left( X_{i}\left( t\right) , h_{i}\right) =\gamma _{i}\frac{h_{i}^{2}}{X_{i}\left( t\right) }, \end{aligned}$$

(3)

denotes harvesting costs and $\gamma _{i}>0$ represents an inefficiency parameter for catching species i, see Clark (1990), Conrad and Smith (2012) and the Appendix in Bischi et al. (2013a) for details.

Profit maximization leads to a myopic Markovian harvesting rule of the type

$$\begin{aligned} h_{i}^{*}\left( X_{i}\left( t\right) \right) =\frac{a_{i}X_{i}\left( t\right) }{2\gamma _{i}}\text {, }\quad i=1,2, \end{aligned}$$

(4)

By inserting (4) into (2), individual instantaneous profits assume the following form

$$\begin{aligned} \pi _{i}^{*}\left( X_{i}\left( t\right) \right) =\frac{\gamma _{i}}{X_{i}\left( t\right) }\left( h_{i}^{*}\left( X_{i}\left( t\right) \right) \right) ^{2}=\frac{a_{i}^{2}X_{i}\left( t\right) }{4\gamma _{i}}\text {, } \quad i=1,2, \end{aligned}$$

(5)

As mentioned before, fishers target only one species during a certain period of time. Therefore, considering the set of N agents, in a given time interval, a fraction $r_{i}$ specializes in the harvesting of species i, so that total instantaneous harvesting in (1) can be written as follows

$$\begin{aligned} H_{i}\left( X_{i}\left( t\right) ,r_{i}\left( t\right) \right) = N r_{i}\left( t\right) h_{i}^{*}\left( X_{i}\left( t\right) \right) = N r_{i}\left( t\right) \frac{a_{i}X_{i}\left( t\right) }{2\gamma _{i}}. \end{aligned}$$

(6)

In the next section, we describe how fractions $r_{i}$ are updated over time according to evolutionary pressure.

3 Switching mechanism and bio-economic model

In this section, we describe in detail the evolutionary mechanism that determines the temporal evolution of the fraction of agents that target each species by following and presenting a possible extension of the bio-economic model in Bischi et al. (2013a). The main reason for this choice is to model a specialized fishery, characterized by harvesting that concentrates on a particular species for a time interval.

As mentioned above, a share $r_{i}$ of fishers catches only species i during an interval of time. Harvesting takes place continuously and agents select the harvesting according to (4). Thus, during a period of time, a fixed share $r_{i}$ of fishers harvest only species i characterized by logistic growth, so that the time evolution of species i in (1) can be rewritten as

$$\begin{aligned} {\dot{X}}_{i}\left( t\right) = X_{i}\left( t\right) \rho _{i}\left( 1-\theta _{i}\left( r_{i}\left( t\right) \right) -\frac{X_{i}\left( t\right) }{k_{i}}\right) \text {, } \quad i=1,2, \end{aligned}$$

(7)

where the term $\rho _{i}\theta _{i}\left( r_{i}\right)$, with $\theta _{i}\left( r_{i}\right) =\frac{Nr_{i}}{\rho _{i}}\frac{a_{i}}{2\gamma _{i}}$, is the instantaneous harvesting for unit of biomass of species i.

In Bischi et al. (2013a), the continuous-time dynamics of resources in (7) is coupled with a discrete-time dynamics for the evolution of the share of fishers devoted to harvest only species i, $i=1,2$, that is $r_{i}$. This choice is carried out by considering the average profit $\overline{\overline{\pi }}_{i}^{*}$ over a period of length $s>0$ for a representative agent targeting species i. Specifically, the value of $\overline{\overline{\pi }}_{i}^{*}$ at time $t+s$ is given by

$$\begin{aligned} \overline{\overline{\pi }}_{i}^{*}\left( X_{i}\left( t\right) ,r_{i}\left( t\right) \right) =\frac{\displaystyle \int \limits _{0}^{s}\pi _{i}^{*}\left( X_{i}\left( t+\tau \right) \right) d\tau }{s}\text {, } \quad i=1,2, \end{aligned}$$

(8)

By assuming that the magnitude of $\overline{\overline{\pi }}_{i}^{*}$ can be assessed by agents of either group, Bischi et al. (2013a) proposes a hybrid dynamical system obtained by coupling continuous-time growth and harvesting of the fish species in (7) with a discrete (or pulse) fishing strategy switching, modeled by discrete replicator dynamics (see Weibull (1995), Hofbauer and Sigmund (1998) and Cressman (2003)) of the form

$$\begin{aligned} r\left( t+s\right) =\left\{ \begin{array}{lr} \displaystyle \frac{r\left( t\right) \overline{\overline{\pi }}_{1}^{*}\left( X_{1}\left( t\right) ,r\left( t\right) \right) }{r\left( t\right) \overline{\overline{\pi }}_{1}^{*}\left( X_{1}\left( t\right) ,r\left( t\right) \right) + \left( 1-r\left( t\right) \right) \overline{\overline{\pi }}_{2}^{*}\left( X_{2}\left( t\right) ,r\left( t\right) \right) } &{} \quad \text { if } \quad \frac{t+s}{s}=\left\lfloor \frac{t+s}{s}\right\rfloor \\ \\ r\left( \left\lfloor \frac{t+s}{s}\right\rfloor s\right) &{} \quad \text { otherwise } \end{array} \right. ,\nonumber \\ \end{aligned}$$

(9)

where $\left\lfloor x\right\rfloor$ is the largest integer not greater than x (i.e. the floor of x), $\overline{\overline{\pi }}_{i}^{*}$, $i=1,2$ are given, respectively, in (8).

Differently from Bischi et al. (2013a), here we limit the attention to describe the state of the system at the switching times ${\mathbb {S}}=\left\{ 0,s,2s,\ldots ,ns,\ldots \right\}$. Specifically, consider a time t at which fishers decide the biomass to target, that is $t\in {\mathbb {S}}$. Then, at time t the shares $r_{i}$, $i=1,2$, are determined and remain fixed during the next interval of time $\left[ t,t+s\right)$. Moreover, if at time t the biomass of species i measures $X_{i}\left( t\right)$, at time $X_{i}\left( t+s\right)$ its value can be obtained by solving the correspondent Cauchy problem with dynamics (7) and initial condition $X_{i}\left( t\right)$. Thus, assuming that specialized fishers do not change harvested species in the interval of time $\left[ t,t+s\right)$, the total biomass of species i at time $t+s$ reads

$$\begin{aligned} X_{i}\left( t+s\right) =\frac{\left( 1-\theta _{i}\left( r\left( t\right) \right) \right) k_{i}X_{i}\left( t\right) }{X_{i}\left( t\right) -\left[ X_{i}\left( t\right) -\left( 1-\theta _{i}\left( r\left( t\right) \right) \right) k_{i}\right] e^{-\left( 1-\theta _{i}\left( r\left( t\right) \right) \right) \rho _{i}s}}. \end{aligned}$$

(10)

The switching mechanism that describes how the fraction of specialized fishers changes targeted species is based on profit comparison on the relevant interval $\left[ t,t+s\right)$, with $t\in {\mathbb {S}}$. Specifically, extending the main framework in Bischi et al. (2013a), we assume that agents make an assessment of the total profits obtained by either strategy during the time interval $\left[ t,t+s\right)$, introducing a parameter $\delta \in (-1,1)$ that weighs the profits on the basis of the moment of realization. The particular case of $\delta =0$ is the setup employed in Bischi et al. (2013a). Therefore, considering instantaneous profits in (5), at the end of each non-switching period $\left[ t,t+s\right)$, with $t\in {\mathbb {S}}$, a representative agent who targeted species i assesses the evaluation of average profits $\overline{\pi }_{i}^{*}$ (compounded to time $t+s$ when $\delta <0$ or with more weight on recent profits when $\delta >0$) over the last period, by calculating

$$\begin{aligned} \overline{\pi }_{i}^{*}\left( X_{i}\left( t\right) ,r_{i}\left( t\right) \right)= & {} \frac{1}{s}\displaystyle \int \limits _{0}^{s}e^{\delta (\tau -s)}\pi _{i}^{*}\left( X_{i}\left( t+\tau \right) \right) d\tau \nonumber \\= & {} \frac{a_{i}^{2} }{4\gamma _{i}s}\displaystyle \int \limits _{0}^{s}e^{\delta (\tau -s)} X_{i}\left( t+\tau \right) d\tau \text {; } \quad i=1,2 \end{aligned}$$

(11)

where the integrand function is given by the product between the biomass of species i at time $t+\tau \in \left[ t,t+s\right)$, given in (10), times the weighting term $e^{\delta (\tau -s)}$. Clearly when $\delta <0$, $e^{\delta (\tau -s)}$ can be interpreted as a standard compounding factor, which takes into account the interests on the profits obtained by harvesting species i in the time interval, where the evaluation time is the end of each period. The case $\delta >0$ corresponds to negative interests: we consider this possibility because it has two interesting interpretations. Firstly, the phenomenon of having negative interests has been recently observed in the financial world. Secondly, it can be interpreted as a term of fading memory, which places more weight on more recent profits than on older ones, as described in Bischi et al. (2020a) (see also Bischi and Merlone (2017) and Bischi et al. (2020b) for evolutionary games with fading memory). In the particular case of $\delta =0$ we have that $\overline{\pi }_{i}^{*}$ in (11) is equal to $\overline{\overline{\pi }}_{i}^{*}$ in (8), which is the original setup in Bischi et al. (2013a). Note also that as the integrand in (11) is non-negative so are average profits $\overline{\pi }_{i}^{*}$. Expressing $X_{i}\left( t+\tau \right)$, with $\tau \in \left[ 0,s\right]$, as a function of $X_{i}\left( t\right)$ using the formula in (10) and computing the integral in (11), the (discounted/capitalized) average profits in the period $\left[ t,t+s\right)$ can be assessed as

$$\begin{aligned}&\overline{\pi }_{i}^{*}\left( X_{i}\left( t\right) ,r_{i}\left( t\right) \right) \nonumber \\&\quad = \frac{e^{-s\delta } a_{i}^{2}\left( 1-\theta _{i}\left( r_{i}\left( t\right) \right) \right) k_{i}}{4\gamma _{i}\delta s}\left[ e^{\delta s} \,_{2}F_{1}\left( 1,\Phi _{i}\left( r_{i}\left( t\right) \right) ;1+\Phi _{i}\left( r_{i}\left( t\right) \right) ;\Lambda _{i}\left( X_{i}\left( t\right) ,r_{i}\left( t\right) \right) \right) \right. \nonumber \\&\left. \qquad -_{2}F_{1}\left( 1,\Phi _{i}\left( r_{i}\left( t\right) \right) ;1+\Phi _{i}\left( r_{i}\left( t\right) \right) ;e^{-s\left( 1-\theta _{i}\left( r_{i}\left( t\right) \right) \right) \rho _{i}}\Lambda _{i}\left( X_{i}\left( t\right) ,r_{i}\left( t\right) \right) \right) \right] , \end{aligned}$$

(12)

where

$$\begin{aligned} \begin{array}{rcl} \Phi _{i}\left( r_{i}\left( t\right) \right) &{} = &{} \frac{\delta }{\rho _{i}\left( 1-\theta _{i}\left( r_{i}\left( t\right) \right) \right) }\\ \\ \Lambda _{i}\left( X_{i}\left( t\right) ,r_{i}\left( t\right) \right) &{} = &{} \frac{X_{i}\left( t\right) -k_{i}\left( 1-\theta _{i}\left( r_{i}\left( t\right) \right) \right) }{X_{i}\left( t\right) }\\ \end{array}, \end{aligned}$$

(13)

and $_{2}F_{1}\left( \cdot ,\cdot ;\cdot ;\cdot \right)$ is the Gaussian, or ordinary, hypergeometric function. In the particular case $\delta =0$ (setup in Bischi et al. (2013a)), average profits over the non-switching period $\left[ t,t+s\right)$ are given by

$$\begin{aligned}&\overline{\pi }_{i}^{*}\left( X_{i}\left( t\right) ,r_{i}\left( t\right) \right) \nonumber \\&\quad = \frac{a_{i}^{2}k_{i}}{4\gamma _{i}\rho _{i}s}\left\{ \left( 1-\theta _{i}\left( r_{i}\left( t\right) \right) \right) \rho _{i}s\right. \nonumber \\&\left. \qquad +\log \left[ \left( X_{i}\left( t\right) -\left( 1-\theta _{i}\left( r_{i}\left( t\right) \right) \right) k_{i}\right) e^{-\left( 1-\theta _{i}\left( r_{i}\left( t\right) \right) \right) \rho _{i}s}-X_{i}\left( t\right) \right] \right. \nonumber \\&\left. \qquad - \log \left( \left( \theta _{i}\left( r_{i}\left( t\right) \right) -1\right) k_{i}\right) \right\} . \end{aligned}$$

(14)

Notice that average profits over the period $\left[ t,t+s\right)$ depend on $X_{i}\left( t\right)$, the biomass of species i at time t, on $\theta _{i}\left( r_{i}\left( t\right) \right) =\frac{Nr_{i}\left( t\right) }{\rho _{i}}\frac{a_{i}}{2\gamma _{i}}$, where $r_{i}\left( t\right)$ is the share of fishers harvesting species i held constant in $\left[ t,t+s\right)$, and on the length s of the period in which the two shares of fishers are fixed.

Consider the average profits over the period $\left[ t,t+s\right)$ defined as in (11), where $r_{1}=r$ and $r_{2}=1-r$. At each time $t+s$, with $t\in {\mathbb {S}}$, the fraction r is updated accordingly. In fact, average profits are taken as a measure of fitness for the performance of each strategy, so that shares are updated according to replicator dynamics in discrete time, see Weibull (1995) and Hofbauer and Sigmund (1998), that is at each time $t\in \left\{ s,2s,3s,\ldots ,ns,\ldots \right\}$, we have:

$$\begin{aligned} r\left( t+s\right) =\frac{r\left( t\right) \overline{\pi }_{1}^{*}\left( X_{1}\left( t\right) ,r\left( t\right) \right) }{r\left( t\right) \overline{\pi }_{1}^{*}\left( X_{1}\left( t\right) ,r\left( t\right) \right) +\left( 1-r\left( t\right) \right) \overline{\pi }_{2}^{*}\left( X_{1}\left( t\right) ,r\left( t\right) \right) } . \end{aligned}$$

(15)

All in all, with a slight abuse of notation,³ the original bio-economic hybrid model is given by

$$\begin{aligned} \left( {\dot{X}}_{1}\left( t\right) ,{\dot{X}}_{2}\left( t\right) ,r\left( t\right) \right) = L\left( X_{1}\left( t\right) ,X_{2}\left( t\right) ,r\left( t\right) ,X_{1}\left( t-s\right) ,X_{2}\left( t-s\right) ,r\left( t-s\right) \right) ,\nonumber \\ \end{aligned}$$

(16)

where

$$\begin{aligned}&L\left( X_{1}\left( t\right) ,X_{2}\left( t\right) ,r\left( t\right) ,X_{1}\left( t-s\right) ,X_{2}\left( t-s\right) ,r\left( t-s\right) \right) \nonumber \\&\quad = \left\{ \begin{array}{l} X_{1}\left( t\right) \rho _{1}\left( 1-\theta _{1}\left( r\left( t\right) \right) -\frac{X_{1}\left( t\right) }{k_{1}}\right) \\ \\ X_{2}\left( t\right) \rho _{2}\left( 1-\theta _{2}\left( 1-r\left( t\right) \right) -\frac{X_{2}\left( t\right) }{k_{2}}\right) \\ \\ R\left( X_{1}\left( t-s\right) ,X_{2}\left( t-s\right) ,r\left( t-s\right) \right) \end{array} \right. \end{aligned}$$

(17)

and

$$\begin{aligned}&R\left( X_{1}\left( t\right) ,X_{2}\left( t\right) ,r\left( t\right) \right) \nonumber \\&\quad := \left\{ \begin{array}{lr} \displaystyle \frac{r\left( t\right) \overline{\pi }_{1}^{*}\left( X_{1}\left( t\right) ,r\left( t\right) \right) }{r\left( t\right) \overline{\pi }_{1}^{*}\left( X_{1}\left( t\right) ,r\left( t\right) \right) + \left( 1-r\left( t\right) \right) \overline{\pi }_{2}^{*}\left( X_{2}\left( t\right) ,r\left( t\right) \right) } &{} \quad \text { if } \quad \frac{t+s}{s}=\left\lfloor \frac{t+s}{s}\right\rfloor \\ \\ r\left( \left\lfloor \frac{t+s}{s}\right\rfloor s\right) &{} \quad \text { otherwise } \end{array} \right. \nonumber \\ \end{aligned}$$

(18)

to which the following parameter restrictions apply.

It is possible to represent this hybrid model at times $\left\{ 0,s,2s,3s,\ldots ,ns,\ldots \right\}$ through a three-dimensional map by coupling the dynamics of the resources in (10) sampled at the discrete times of switching with the discrete replicator dynamics in (15). Then, the hybrid dynamical system can be reformulated recursively as a three-dimensional map of the form:

$$\begin{aligned} \left( X_{1}\left( t+s\right) ,X_{2}\left( t+s\right) ,r\left( t+s\right) \right) = G\left( X_{1}\left( t\right) ,X_{2}\left( t\right) ,r\left( t\right) \right) , \end{aligned}$$

(19)

where the function G is defined as follows

$$\begin{aligned} G\left( X_{1}\left( t\right) ,X_{2}\left( t\right) ,r\left( t\right) \right) := \left\{ \begin{array}{l} \frac{\left( 1-\frac{Nr\left( t\right) }{\rho _{1}}\frac{a_{1}}{2\gamma _{1}}\right) k_{1}X_{1}\left( t\right) }{X_{1}\left( t\right) +\left[ \left( 1-\frac{Nr\left( t\right) }{\rho _{1}}\frac{a_{1}}{2\gamma _{1}}\right) k_{1}-X_{1}\left( t\right) \right] e^{-\left( 1-\frac{Nr\left( t\right) }{\rho _{1}}\frac{a_{1}}{2\gamma _{1}}\right) \rho _{1}s}}\\ \\ \frac{\left( 1-\frac{N\left( 1-r\left( t\right) \right) }{\rho _{2}}\frac{a_{2}}{2\gamma _{2}}\right) k_{2}X_{2}\left( t\right) }{X_{2}\left( t\right) +\left[ \left( 1-\frac{N\left( 1-r\left( t\right) \right) }{\rho _{2}}\frac{a_{2}}{2\gamma _{2}}\right) k_{2}-X_{2}\left( t\right) \right] e^{-\left( 1-\frac{N\left( 1-r\left( t\right) \right) }{\rho _{2}}\frac{a_{2}}{2\gamma _{2}}\right) \rho _{2}s}}\\ \\ \frac{r\left( t\right) \overline{\pi }_{1}^{*}\left( X_{1}\left( t\right) ,r\left( t\right) \right) }{r\left( t\right) \overline{\pi }_{1}^{*}\left( X_{1}\left( t\right) ,r\left( t\right) \right) +\left( 1-r\left( t\right) \right) \overline{\pi }_{2}^{*}\left( X_{2}\left( t\right) ,r\left( t\right) \right) } \end{array} \right. ,\nonumber \\ \end{aligned}$$

(20)

where $t\in \left\{ 0,s,2s,3s,\ldots ,ns,\ldots \right\}$ and with $\overline{\pi }_{2}^{*}$, $i=1,2$, defined as in (12) if $\delta \ne 0$ and as in (14) if $\delta =0$.

Notice that the length s of the harvesting period in which shares of fishers are fixed (non-switching period) is a parameter of the model. By taking s as the unit of time or by operating the time rescaling $z=\frac{t}{s}$, we can write the dynamical system (19) as a discrete dynamical system in standard form.

By construction of the map G, the related discrete-time dynamical system (19) represents the hybrid dynamical system L at any time $t\in {\mathbb {S}}$. This proves the following proposition:

Proposition 1

At any time $t\in {\mathbb {S}}$, the hybrid dynamical system (16)-(17) admits the discrete time representation given by map $G\left( X_{1}\left( t\right) ,X_{2}\left( t\right) ,r\left( t\right) \right)$ in (19)-(20).

Our conjecture is that any hybrid dynamical system admit a discrete time representation given a suitable construction of ${\mathbb {S}}$. At the same time, it is easy to conclude that a hybrid dynamical system cannot have a representation in terms of a continuous-time dynamical system. Indeed, the function L is not differentiable in r. Therefore, the continuous-time version of the discrete-time system (19), given by

$$\begin{aligned} \left( {\dot{X}}_{1}\left( t\right) ,{\dot{X}}_{2}\left( t\right) ,{\dot{r}}\left( t\right) \right) = F \left( X_{1}\left( t\right) ,X_{2}\left( t\right) ,r\left( t\right) \right) , \end{aligned}$$

(21)

where F is the vector of partial derivatives of G, cannot replicate the value of the hybrid model at any time $t\in {\mathbb {R}}_{+}$.

The advantage of using a discrete-time replicator dynamics to represent the value of our hybrid dynamical system at time $t\in {\mathbb {S}}$ is given by the possibility to explain the local and global dynamics of the hybrid model by means of the bifurcation theory for discrete-time dynamical systems. Therefore, in the next section we can provide a global analysis of the dynamics of the hybrid model and explain the mechanism through which equilibria lose stability. This analysis extends the results in Bischi et al. (2013a), where to provide some analytical results the dynamics of the hybrid system was approximated by means of the continuous-time dynamical system (21) to which we refer in the following as the benchmark model with continuous-time replicator dynamics.

4 Dynamic analysis

The dynamical system (19) is analyzed in the set

$$\begin{aligned}&A=\left( \left[ 0,+\infty \right) \times \left[ 0,+\infty \right) \times \left[ 0,1\right] \right) /\left( \left\{ \left. \left( 0,X_{2},0\right) \right| X_{2}\in \left[ 0,+\infty \right) \right\} \right. \nonumber \\&\left. \quad \cup \left\{ \left. \left( X_{1},0,1\right) \right| X_{1}\in \left[ 0,+\infty \right) \right\} \cup \left\{ \left. \left( 0,0,r\right) \right| r\in \left[ 0,1\right] \right\} \right) . \end{aligned}$$

(22)

The region of interest A is an invariant set for the dynamical system (19), which admits some boundary equilibria under specific conditions on the value of the parameters. A boundary equilibrium is an equilibrium of the system that lies in one of the borders of A. Possible boundary equilibria of the dynamical system (19) are

$$\begin{aligned} E_{1}^{1}=\left( 0,k_{2}\left( 1-N\frac{a_{2}}{2\gamma _{2}\rho _{2}}\right) ,0\right) \quad \text {and}\quad E_{2}^{1}=\left( k_{1}\left( 1-N\frac{a_{1}}{2\gamma _{1}\rho _{1}}\right) ,0,1\right) ,\nonumber \\ \end{aligned}$$

(23)

characterized by the extinction of the non-harvested species. These boundary equilibria $E_{i}^{1}$, with $i=1,2$, exist only under the specific condition $Na_{i}<2\gamma _{i}\rho _{i}$, with $i=1,2$, so that $E_{i}^{1}$ has positive biomass in the harvested species i.

The same condition $Na_{i}<2\gamma _{i}\rho _{i}$, with $i=1,2$, ensures the existence of another boundary equilibrium, which is characterized by harvesting of species i only, and with the other species being unexploited and at carrying capacity. Consider $i=1,2$, these further boundary equilibria are

$$\begin{aligned} E_{1}^{2}=\left( k_{1},k_{2}\left( 1-N\frac{a_{2}}{2\gamma _{2}\rho _{2} }\right) ,0\right) \quad \text { and } \quad E_{2}^{2}=\left( k_{1}\left( 1-N\frac{a_{1}}{2\gamma _{1}\rho _{1}}\right) ,k_{2},1\right) .\nonumber \\ \end{aligned}$$

(24)

The local stability of these boundary equilibria is discussed in the following proposition (the proof of which is in Appendix A).

Proposition 2

(Local stability of boundary equilibria) Consider the dynamics of system (19) in A. The following results hold:

• The boundary equilibria with extinction of the non-harvested species $E_{1}^{1}$ and $E_{2}^{1}$, see (23), are always saddle points, provided they exist in A.

• The boundary equilibria with harvesting of species 2 only $E_{1} ^{2}$, see (24), is a stable node when the ratio of selling prices is such that

  $$\begin{aligned} \frac{a_{1}}{a_{2}}<\sqrt{\frac{\gamma _{1}k_{2}\left( 2\gamma _{2}\rho _{2}-Na_{2}\right) }{2\gamma _{2}k_{1}\rho _{2}}}, \end{aligned}$$

  (25)

  and is a saddle point otherwise.

• The boundary equilibria with harvesting of species 1 only $E_{2}^{2}$, see (24), is a stable node when the ratio of selling prices is such that

  $$\begin{aligned} \frac{a_{2}}{a_{1}}<\sqrt{\frac{\gamma _{2}k_{1}\left( 2\gamma _{1}\rho _{1}-Na_{1}\right) }{2\gamma _{1}k_{2}\rho _{1}}}, \end{aligned}$$

  (26)

  and is a saddle point otherwise.

In addition to the previous kinds of equilibria, there exists an inner equilibrium with the exploitation of both species when an iso-profit condition holds, that is when average profits over each period are equal for the two species and the biomasses are in equilibrium. As stated in the following proposition (the proof of which is in Appendix A), the inner equilibrium is unique.

Proposition 3

Consider the dynamics of system (19) in the set A. A unique inner equilibrium $\left( X_{1}^{*},X_{2}^{*},r^{*}\right)$ exists, where

$$\begin{aligned} X_{i}^{*}= & {} \frac{a_{j}^{2}k_{1}k_{2}\gamma _{i}\left( 2a_{2}\gamma _{1}\rho _{1}+2a_{1}\gamma _{2}\rho _{2}-a_{1}a_{2}N\right) }{2\left( a_{2}^{3}k_{2}\gamma _{1}^{2}\rho _{1}+a_{1}^{3}k_{1}{\gamma }_{2}^{2}\rho _{2}\right) }\text {, }\quad i=1,2;i\ne j\nonumber \\ \nonumber \\ r^{*}= & {} \frac{\gamma _{1}\rho _{1}\left( a_{2}^{3} k_{2}N\gamma _{1}-2a_{2}^{2}k_{2}\gamma _{1}\gamma _{2}\rho _{2}+2a_{1}^{2}k_{1}\gamma _{2}^{2}\rho _{2}\right) }{N\left( a_{2}^{3} k_{2}\gamma _{1}^{2}\rho _{1}+a_{1}^{3}k_{1}\gamma _{2}^{2}\rho _{2}\right) }. \end{aligned}$$

(27)

The conditions for existence and stability of equilibria of the discrete-time dynamical system (19) apply also to the hybrid model (16) since the two systems are equivalent in ${\mathbb {S}}$ as specified in Proposition 1. Therefore, the following result follows straightforward from Propositions 1 and 3.

Corollary 1

The hybrid dynamical system (16) has a unique inner equilibrium $\left( X_{1}^{*},X_{2}^{*},r^{*}\right)$ defined as in (27).

The result in Corollary 1 proves the conjecture in Bischi et al. (2013a) about the uniqueness of the inner equilibrium for the hybrid dynamical system (16). Moreover, the unique inner equilibrium in (27) coincides with the unique inner equilibrium of the benchmark model with continuous replicator analyzed in Bischi et al. (2013a), to which we refer for details on the conditions on the parameters that guarantee their existence. In particular, we here recall that, according to Bischi et al. (2013a), in the benchmark model with continuous-time replicator dynamics whenever equilibrium (27) exists it is locally asymptotically stable.

We take this as a starting point for our dynamic analysis. From the equation of total harvesting for species i, see (6), it is immediate to observe that the strain on the species increases as its price increases. Given a starting scenario, we then consider the effect of a price increment of one species, say species 1, all other parameters fixed. The bifurcation diagrams in Fig. 1 depicts this case with $a_{1}\in \left[ 0,300\right]$ and other parameters set as in the baseline example in Bischi et al. (2013a), that are

$$\begin{aligned} \rho _{1}=80;k_{1}=50;\gamma _{1}=9;a_{2}=60;\rho _{2}=140;k_{2}=80;\gamma _{2}=9;N=10; \end{aligned}$$

(28)

and with a length of the interval $s=1.5$ and $\delta =0$. When harvesting species 1 is not remunerative enough, that is its price is well below the price for species 2, species 1 is unexploited and at around its carrying capacity $k_{1}$. In this case, the system converges to corner equilibrium $E_{1}^{2}$ in which $r=0$ that is all fishers target species 2, which converges at its modified carrying capacity.

A transcritical bifurcation occurs at $a_{1}\approx 66.24$, so that a stability exchange occurs between the corner equilibrium $E_{1}^{2}$ and the inner equilibrium $\left( X_{1}^{*},X_{2}^{*},r^{*}\right)$ in (27) that becomes stable. For higher prices of species 1, its exploitation becomes highly remunerative and, as its price increases, its biomass continuously decreases with price increments. At the same time, as more fishers target species 1, the harvesting pressure on species 2 is reduced and the biomass of species 2 increases. Here there is an interesting biological interaction effect between the two species that is due to specialized harvesting. In fact, as prices and landing of species 1 increase, less pressure is placed on species 2 which can grow more freely. As a consequence of the stock externality in the cost function (3), harvesting species 2 becomes cheaper to be withdrawn. In other words, the abundance of species 2 here entails a lower cost of its harvesting which brings higher profits for fishing species 2 despite its relatively low price. The effect is that, for a price of about $a_{1}\approx 127.6$, the fraction of agents who find it convenient to fish species 1 reaches a maximum and then it begins to decrease as the price increases. As a consequence, also the biomass of species 2 has a maximum point although biomass of species 1 continues to decrease as its price increases, see the bifurcation diagrams in Fig. 1(a)-(b), (the bifurcation diagram for the state variable r is shown in Fig. 4). At a price $a_{1}\approx 193.467$, the inner equilibrium $\left( X_{1}^{*},X_{2}^{*},r^{*}\right)$ in (27) loses stability through a flip bifurcation: the selling price of species $a_{1}$ is relatively high but, given the scarcity of species 1, so is its cost of harvesting. The agents’ attempt to harvest the resource of greater economic value, combined with its scarcity and the consequent low cost of the withdrawal of the species 2, gives rise to endogenous instability due to the tendency of fishers to move towards the more profitable species. The instability, originating from a cascade of flip bifurcations, leads the system to chaotic dynamics for species 1 and on species 2 through the dynamics of the fraction of specialized fishers r. Note that with the onset of equilibrium instability, and in particular with the first cycle of period four, cyclical or chaotic dynamics imply not only greater fluctuations in the resource biomass but also extremely low biomass values. In the model, extinction is avoided because biomass can only assume strictly positive values due to the logistic dynamics in (7), but in real cases, the resource (species 1 in this example) can become extinct because of small exogenous shocks.

Fig. 1

One-dimensional bifurcation diagrams with respect to the parameter $a_{1}$ ranging in the interval $\left[ 0,300\right]$. In Panel (a), resp. (b), the long run dynamics of the state variable $X_{1}\left( t\right)$, resp. $X_{2}\left( t\right)$, is represented as a function of the bifurcation parameter $a_{1}$. For $a_{1}=0$, the long run dynamics is represented by the equilibrium $E_{1}^{2}$. For $a_{1}\approx 66.24$ the long run dynamics is represented by the inner equilibrium $\left( X_{1}^{*},X_{2}^{*},r^{*}\right)$ in (27). At $a_{1}\approx 193.467$, the inner equilibrium undergoes a flip bifurcation and for $a_{1}>193.467$ a sequence of period-doubling bifurcations leads to chaotic dynamics. Parameters as in (28) and with a length of the interval $s=1.5$ and $\delta =0$

Let us now take such a chaotic scenario for an exploration of the influence of the length of the non-switching period s on the overall dynamics. For this purpose, we maintain all parameters as in the previous case, see (28), with $a_{1}=300$ and again without risk-free interest accumulation ($\delta =0$). In Bischi et al. (2013a), it is shown that under continuous-time adjustments in biomass and in the replicator dynamics (benchmark model with continuous-time replicator dynamics), the inner equilibrium, whenever it exists, is locally asymptotically stable. With the parameters of the example, the inner equilibrium in (27) is given by

$$\begin{aligned} \left( X_{1}^{*},X_{2}^{*},r^{*}\right) =(2.78345,69.5863,0.45328). \end{aligned}$$

(29)

From (27), the length s of the non-switching period has no influence on equilibrium values but, it strongly influences the dynamic behavior of the system in the hybrid case, as shown below. Indeed, differently from the benchmark model with continuous replicator dynamics considered in Bischi et al. (2013a), under hybrid dynamics with strong pressure on species 1 the inner equilibrium may be unstable and chaotic dynamics arise. Specifically, with species 1 that oscillates in the interval $\left( 0,k_{1}\right)$ with periodic windows of low periodicity. However, the increment of the length of the non-switching interval can have a stabilizing effect, which is evident for switching times greater than $s\approx 1.52$, as the generic trajectory converges to a cycle of period 3. Observe also another interesting phenomenon of coexistence of attractors, which is evident in the bifurcation diagram of Fig. 2(a)-(b) for switching times in the interval $\left( 1.75,1.82\right)$, where blue and red trajectories are obtained with initial conditions $\left( X_{1}(0),X_{2}(0),r(0)\right)$ given by $\left( 30,40,0.2\right)$ and $\left( 30,40,0.4\right)$ respectively. In this example, with the same initial biomass for the two species, a difference in the initial share of harvesters determines the long-run dynamics convergence to a cycle of period 3 (blue points) or to a complex attractor or to a different cycle (red points), which is a clear phenomenon of path-dependence. Note that on the cycle of period three, the biomass of species 1 approaches extremely low levels and therefore the danger of extinction of the resource is greater than on the chaotic attractor, in which the biomass never decreases below a threshold value (and also the oscillations of the resource are lower along the chaotic attractor than on the 3-cycle).

Fig. 2

One-dimensional bifurcation diagrams with respect to the parameter s ranging in the interval $\left[ 0,2\right]$. In Panel (a), resp. (b), the long run dynamics of the state variable $X_{1}\left( t\right)$, resp. $X_{2}\left( t\right)$, is represented as a function of the bifurcation parameter s. For s greater than around 1.75 there are two attractors. On attractor is depicted in blue, the other attractor is depicted in red. Parameters: $a_{1}=300$, $\delta =0$ and the remaining values of the parameters as in (28) (color figure online)

Figure 3 shows the kind of attractors obtained by varying the length s of the non-switching interval. For sufficiently low values, the generic trajectory converges to a limit cycle, with high variability for species 1 but very low variability for species 2, see Fig. 3(a) obtained with all parameters as in Fig. 2, and $s=0.001$. An increment in s leads to the evident formation of a chaotic attractor, see Fig. 3(b)-(c) obtained with $s=0.5$ and $s=1$ respectively. Compared to the case with instantaneous adjustments in the replicator dynamics, described in Bischi et al. (2013a), the presence of a finite and small amplitude of the non-switching period leads to cyclical oscillations that mainly impact the species of high commercial value. However, longer periods of specialization on a resource lead to fluctuations even in the species of lower commercial value. Figure 3(d) obtained for $s=1.75$ shows the coexistence of attractors which was evident in Fig. 2: in this case, a low period attractor (red points, obtained with initial condition $\left( 40,75,0.5\right)$) coexists with a three-piece chaotic attractor (blue points, obtained with initial condition $\left( 30,65,0.2\right)$).

Fig. 3

Phase spaces and long run trajectories of the dynamical system (19). Parameters as in 2, that is $a_{3}=300$ and the remaining parameter values as in (28). Panel (a), $s=0.001$. Panel (b), $s=0.5$. Panel (c), $s=1$. Panel (d), $s=1.75$

Now we consider the effect of weighting differently past profits. As already observed before, from (27) the instantaneous risk-free interest rate $\delta$, as well as the length of the non-switching period s, have no influence on equilibrium values but they do influence its stability. This effect was shown for s in the bifurcation diagram of Fig. 2. Figure 4 re-proposes the example of Fig. 1 when interests are taken into account. Here all parameters are as in Fig. 1 but $\delta =0.3$, that is more weight is placed on recent profits along the non-switching period. Here for comparison purposes, we show how the share r of harvesters of species 1 varies with the selling price of species 1 in the case without interest calculations (blue points) and with interest calculations (red points). When convergence to equilibrium occurs, the dynamics are identical. However, instability occurs for lower price levels when interests are neglected and also the bifurcation points subsequent to the point of the first flip bifurcation are obtained for lower values of prices when interests are disregarded. In any case, the type of complexity observed is similar, given that the road to chaos still opens up with a cascade of flip bifurcations.

Fig. 4

One-dimensional bifurcation diagrams with respect to the parameter $a_{1}$ ranging in the interval $\left[ 0,300\right]$. The bifurcation shows the long run dynamics of the state variable $X_{2}\left( t\right)$ that is represented as a function of the bifurcation parameter $a_{1}$. The blue trajectory is obtained setting $\delta =0$, the red trajectory is obtained setting $\delta = 0.3$. The remaining parameters are as in (28) and $s=1.5$, that is as in Fig. 1 (color figure online)

The last numerical example we propose tries to briefly explore the effect of having positive interests on profits ($\delta <0$) or negative interests/fading memory ($\delta >0$) for the dynamics of the system. Figure 5 presents a bifurcation diagram for our baseline example with $\delta \in \left( -1,1\right)$ and all other parameters as in (28) with $a_{1}=205$, so that high landing pressure is imposed on species 1, and $s=8$, so that a sufficiently long interval in which interests are calculated is imposed. Here it is evident that the effect of memory is not univocal: for sufficiently low or high positive values of $\delta$, the dynamics converge to a cycle (respectively of period 3 and of period 2), while for intermediate values of $\delta$, mainly chaotic dynamics occurs with periodicity windows. Here we have shown a typical scenario in the case of the instability of the inner equilibrium, otherwise the presence of memory would not have led to any difference in the dynamics of the system. Note that higher levels of memory (higher positive values of $\delta$) are associated with a lower variability of both resources subject to withdrawal. This is due to the fact that memory on profits places more weight on the recent performances and this creates an attitude that reduces harvesting pressure. The effect that we observe from the bifurcation diagrams in Fig. 5 is a reduction of fluctuations in the level of biomass. Opposite is the effect with a negative $\delta$, that is when interests on past profits are considered: the higher the interest rate on profits (lower negative values of $\delta$), the greater the incentive for agents to withdraw resources earlier, thus reducing their ability to regenerate.

Fig. 5

One-dimensional bifurcation diagrams with respect to the parameter $\delta$ ranging in the interval $\left[ -1,1\right]$. In Panel (a), resp. (b), the long run dynamics of the state variable $X_{1}\left( t\right)$, resp. $X_{2}\left( t\right)$, is represented as a function of the bifurcation parameter $\delta$. Reducing $\delta$ we observe a sequence of period-doubling bifurcations. Parameters as in (28) with $s=8$ and $a_{1}=300$

5 Conclusions

In this paper, we have analyzed a dynamic model that represents the dynamics of two renewable natural resources subject to harvesting. The model is hybrid since it seeks to consistently represent biological dynamics, obedient to continuous-time dynamics, with specialized fishing dynamics, which by its nature can only be reviewed at discrete temporal instants. The model takes its cue from Bischi et al. (2013a) by introducing a weighting function (that can be interpreted as a risk-free interest rate or a rate of fading memory) in the profit accumulation and discretizing the continuous dynamics of logistic growth by reformulating the model through a three-dimensional map of which the main dynamic properties are presented here. In particular, it is possible to determine the bifurcation values of the inner equilibrium and the paths that lead to greater dynamic complexity. As is natural, greater pressure in terms of a greater commercial value of a resource leads to greater exploitation which, therefore, leads to chaotic fluctuations of all resources involved in the specialized harvesting. From this point of view, even indirectly, the increase in the price of a resource can lead to fluctuations in another resource whose price remains constant if the latter is included within the same specialized fishery. Alongside the phenomena of coexistence and path-dependency, we find in this reformulation that the increase in the width of the interval in which harvesting remains focused on a species does not have a unique impact on the conservation of the resource, but its increase could be beneficial. Furthermore, we find that investing profits at a risk-free rate can have a less conservative and more fluctuating effect on resources, whereas the opposite occurs when fading memory is assumed. The model in question is a possible starting point for further investigations on hybrid modeling in the exploitation of resources and provides a possible methodological contribution for the reformulation of some existing models in a more coherent way with the underlying bio-economic reality.

Appendix A

Proof of Proposition 2

Let us prove the stability conditions for equilibria $E_{1}^{1}$ and $E_{1}^{2}$, as for equilibria $E_{2}^{1}$ and $E_{2}^{2}$ it is only necessary to swap the indices. Equilibrium $E_{1}^{1}=\left( 0,k_{2}\left( 1-N\frac{a_{2}}{2\gamma _{2}\rho _{2}}\right) ,0\right)$ is feasible when biomass of species 2 is nonnegative, that is when $\frac{Na_{2}}{2\gamma _{2}\rho _{2}}<1$. In $E_{1}^{1}$ the Jacobian matrix of the system (19) is

$$\begin{aligned} J\left( E_{1}^{1}\right) =\left( \begin{array} [c]{ccc} \displaystyle e^{s\rho _{1}} &{} 0 &{} 0\\ \\ 0 &{} e^{-\rho _{2}s+\frac{Na_{2}s}{2\gamma _{2}}} &{} \frac{Na_{2}k_{2}\left( 1-e^{\frac{Na_{2}s}{2\gamma _{2}}-\rho _{2}s}\right) }{2\gamma _{2}\rho _{2}}\\ \\ \displaystyle 0 &{} 0 &{} \frac{\overline{\pi }_{1}^{*}\left( 0,0\right) }{\overline{\pi }_{2}^{*}\left( k_{2}\left( 1-N\frac{a_{2}}{2\gamma _{2}\rho _{2}}\right) ,0\right) } \end{array} \right) , \end{aligned}$$

(30)

where $\overline{\pi }_{1}^{*}$ and $\overline{\pi }_{2}^{*}$ are defined as in (11). Being the Jacobian matrix upper triangular, the eigenvalues are $\lambda _{1}^{E_{1}^{1}} = e^{s\rho _{1}}$, which is always greater than one since s and $\rho _{1}$ are positive parameters, $\lambda _{2}^{E_{1}^{1}} = e^{-\rho _{2}s+\frac{Na_{2}s}{2\gamma _{2}}}$, which is always in $\left( 0,1\right)$ as $E_{1}^{1}\in A$ implies $Na_{2}<2\gamma _{2}\rho _{2}$, and $\lambda _{3}^{E_{1}^{1}} = 0$, since $\overline{\pi }_{1}^{*}=0$ any time $X_{1}=0$. This proves that equilibrium $E_{1}^{1}$ whenever exists is always a saddle point. Concerning equilibrium $E_{1}^{2}=\left( k_{1},k_{2}\left( 1-N\frac{a_{2}}{2\gamma _{2}\rho _{2}}\right) ,0\right)$ is meaningful when $E_{1}^{1}$ exists, that is when $\frac{Na_{2}}{2\gamma _{2}\rho _{2}}<1$. At this equilibrium the Jacobian matrix of system (19) is

$$\begin{aligned} J\left( E_{1}^{2}\right) =\left( \begin{array} [c]{ccc} \displaystyle e^{-s\rho _{1}} &{} 0 &{} \frac{a_{1}k_{1}N\left( e^{-s\rho _{1}}-1\right) }{2\gamma _{1}\rho _{1}}\\ \\ \displaystyle 0 &{} e^{-\rho _{2}s+\frac{Na_{2}s}{2\gamma _{2}}} &{} \frac{Na_{2}k_{2}\left( 1-e^{\frac{Na_{2}s}{2\gamma _{2}}-\rho _{2}s}\right) }{2\gamma _{2}\rho _{2}}\\ \\ \displaystyle 0 &{} 0 &{} \frac{\overline{\pi }_{1}^{*}\left( k_{1},0\right) }{\overline{\pi }_{2}^{*}\left( k_{2}\left( 1-N\frac{a_{2}}{2\gamma _{2}\rho _{2}}\right) ,0\right) } \end{array} \right) , \end{aligned}$$

(31)

where $\overline{\pi }_{1}^{*}$ and $\overline{\pi }_{2}^{*}$ are defined as in (11). Being the Jacobian upper triangular, the eigenvalues are $\lambda _{1}^{E_{1}^{2}} = e^{s\rho _{1}}$, which belongs to $\left( 0,1\right)$ since s and $\rho _{1}$ are positive parameters, $\lambda _{2}^{E_{1}^{2}} = e^{-\rho _{2}s+\frac{Na_{2}s}{2\gamma _{2}}}$, which is always in $\left( 0,1\right)$ as $E_{1}^{2}\in A$ implies $Na_{2}<2\gamma _{2}\rho _{2}$, and $\lambda _{3}^{E_{1}^{2}} = \frac{\overline{\pi }_{1}^{*}\left( k_{1},0\right) }{\overline{\pi }_{2}^{*}\left( k_{2}\left( 1-N\frac{a_{2}}{2\gamma _{2}\rho _{2}}\right) ,0\right) }$, where

$$\begin{aligned} \overline{\pi }_{1}^{*}\left( k_{1},0\right) = \frac{a_{1}^{2}k_{1}\left( e^{s\delta }-1\right) }{4\gamma _{1}\delta }, \end{aligned}$$

(32)

and

$$\begin{aligned} \overline{\pi }_{2}^{*}\left( k_{2}\left( 1-N\frac{a_{2}}{2\gamma _{2}\rho _{2}}\right) ,0\right) = \frac{a_{2}^{2}k_{2}\left( 2\gamma _{2}\rho _{2}-a_{2}N\right) \left( e^{s\delta }-1\right) }{8\gamma _{2}^{2}\rho _{2}\delta }>0, \end{aligned}$$

(33)

which is positive any time $E_{1}^{2}\in A$. Therefore, $\lambda _{3}^{E_{1}^{2}}$ is in modulus less than one and $E_{1}^{2}$ is a stable node provided that condition (25) holds. Clearly, when in (25) the reverse inequality holds, equilibrium $E_{1}^{2}$ is a saddle. This completes the proof. $\square$

Proof of Proposition 3

From the equilibrium condition of the biomass, it is immediate to observe from (10) that the species i must be at the modified carrying capacity

$$\begin{aligned} X_{i}^{*}=k_{i}\left( 1-\theta _{i}\left( r_{i}^{*}\right) \right) =k_{i}\left( 1-\frac{Nr_{i}^{*}}{\rho _{i}}\frac{a_{i}}{2\gamma _{i}}\right) \text {, } \quad i=1,2, \end{aligned}$$

(34)

where $r_{i}^{*}$ is the equilibrium share of harvesters targeting species i. Employing the equilibrium biomass in the calculation of average profits in (11), it is then immediate to obtain the corresponding equilibrium profits for species i, that are

$$\begin{aligned} \overline{\overline{\pi }}_{i}^{*}\left( X_{i}^{*},r^{*}\right) =\frac{a_{1}^{2}\left( 1-\frac{Nr_{i}^{*}}{\rho _{i}}\frac{a_{i}}{2\gamma _{i}}\right) k_{i}}{4\gamma _{i}} , \end{aligned}$$

(35)

when $\delta =0$, and

$$\begin{aligned} \overline{\pi }_{i}^{*}\left( X_{i}^{*},r^{*}\right) =\frac{a_{i}^{2}\left( 1-\frac{Nr_{i}^{*}}{\rho _{i}}\frac{a_{i}}{2\gamma _{i}}\right) k_{i}\left( e^{\delta s}-1\right) }{4\gamma _{i}\delta s}=\overline{\overline{\pi }}_{i}^{*}\left( X_{i}^{*},r^{*}\right) \frac{e^{\delta s}-1}{\delta s}, \end{aligned}$$

(36)

when $\delta >0$. By (15), an equilibrium share $r_{i}^{*}$ of harvesters for species i is obtained by imposing the iso-profit condition $\overline{\pi }_{1}^{*}\left( X_{1}^{*},r^{*}\right) =\overline{\pi }_{2}^{*}\left( X_{2}^{*},r^{*}\right)$ which shows the uniqueness of the equilibrium (27) for the dynamical system (19). This completes the proof. $\square$